Changeset 737960d in git
 Timestamp:
 Mar 1, 2021, 12:24:11 PM (3 years ago)
 Branches:
 (u'fiekerDuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
 Children:
 c95beaa99dc4305cb7dbe128571600a256952295
 Parents:
 31f12ca584f1b01f47ae377d9783a9818653141c
 gitauthor:
 Hans Schoenemann <hannes@mathematik.unikl.de>20210301 12:24:11+01:00
 gitcommitter:
 Hans Schoenemann <hannes@mathematik.unikl.de>20210301 12:24:41+01:00
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/grobcov.lib
r31f12ca r737960d 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="version grobcov.lib 4.2.0.0 Derc_2020 "; // $Id$ 3 // version N10; June 2020; 2 version="version grobcov.lib 4.2.0 February_2021 "; // $Id$; 3 // version N12; February 2021; 4 4 5 info=" 5 6 LIBRARY: grobcov.lib … … 120 121 hypothesis and thesis in ADGT. 121 122 122 This version was finished on 26/5/2020, 123 123 The last version N11 (2021) has improved the routines for locus 124 and allows to determine a parametric locus. 125 126 This version was finished on 1/2/2021, 124 127 125 128 NOTATIONS: Before calling any routine of the library grobcov, … … 176 179 177 180 locus(G); Special routine for determining the geometrical 178 locus of points verifying given conditions. Given 179 a parametric ideal J in Q[x1,..,xn][u1,..,um] 180 with parameters (x1,..,xn) and variables 181 (u1,..,um), representing the system determining 182 the ndimensional locus with tracer point (x1,..,xn) 183 verifying certain properties, 184 one can apply locus to the system F, for obtaining the locus. 181 locus of points verifying given conditions. To use it, the 182 ring R=(0,a1,..,ap,x1,..xn),(u1,..um,v1..vn),lp; 183 must be declared, where 184 (a1,..ap) are parameters (optative), 185 (x1,..xn) are the variabes of the tracer point, 186 (u1,..,um) are auxiliary variables, 187 (v1,..,vn) are the mover variables. 188 Then the input to locus must be the 189 parametric ideal F defined in R. 185 190 186 191 locus provides all the components of the locus and … … 188 193 \"Normal\", \"Special\", \"Accumulation\", 189 194 \"Degenerate\". 190 The mover are the uvariables.195 The mover variables are the last n variables. 191 196 The user can ventually restrict them to a subset of them 192 197 for geometrical reasons but this can change the true 193 198 taxonomy. 199 locus also allows to determine a parmetric locus 200 depending on p parameters a1,..ap using then 201 the option \"numpar\",p. 194 202 195 203 locusdg(G); Is a special routine that determines the … … 203 211 Q[x1,..,xn][t1,..,tm] depending on an ideal of 204 212 constraints C in Q[t1,..,tm]. It computes the 205 locus of the envelop, and dete rmines the206 different components as well as its taxonomy:213 locus of the envelop, and detemines the 214 different components as well as their taxonomies: 207 215 \"Normal\", \"Special\", \"Accumulation\", 208 216 \"Degenerate\". (See help for locus). … … 257 265 If H1=1 then H1 is not considered, and analogously for T1. 258 266 259 ConsLevels(A); Given a list of locally c losed sets, constructs the267 ConsLevels(A); Given a list of locally colsed sets, constructs the 260 268 canonical representation of the levels of A an its complement. 261 269 … … 313 321 // Wibmer's Lemma: 1992016 314 322 // Final data October 2016 315 // Updated locus (messages) // Updated locus (messages)323 // Updated locus (messages) 316 324 // Final data Mars 2017 317 325 // Release N4: (public) … … 323 331 // Release N10: May 2020. 324 332 // Updated locus. New determination of the taxonomies 333 // Release N11: February 2021. 334 // Improved the routines for locus. Accept parametric locus as option. 325 335 326 336 //*************Auxiliary routines************** … … 966 976 mu=mu*nu; 967 977 rho=qlc[1]*qlm; 978 //"T_nu="; nu; "mu="; mu; "rho="; rho; 968 979 p=nu*prho*F[i]; 969 980 r=nu*r; 970 981 for (j=1;j<=size(F);j++){q[j]=nu*q[j];} 971 982 q[i]=q[i]+rho; 983 //"T_q[i]="; q[i]; 972 984 divoc=1; 973 985 } … … 979 991 p=plcp*lpp; 980 992 } 993 //"T_r="; r; "p="; p; 981 994 } 982 995 list res=r,q,mu; … … 984 997 } 985 998 example 986 { 987 "EXAMPLE:"; echo = 2; 999 { "RXAMPLE:";echo = 2; 1000 // Division of a polynom by an ideal 1001 988 1002 if(defined(R)){kill R;} 989 1003 ring R=(0,a,b,c),(x,y),dp; 990 1004 short=0; 1005 991 1006 // Divisor 992 1007 poly f=(abac)*xy+(ab)*x+(5c); 1008 993 1009 // Dividends 994 ideal F=ax+b,cy+a; 1010 ideal F=ax+b, 1011 cy+a; 1012 995 1013 // (Remainder, quotients, factor) 996 1014 def r=pdivi(f,F); 997 1015 r; 1016 998 1017 // Verifying the division 999 1018 r[3]*f(r[2][1]*F[1]+r[2][2]*F[2]+r[1]); … … 1190 1209 }; 1191 1210 example 1192 { 1193 "EXAMPLE:"; echo = 2; 1194 if(defined(R)){kill R;} 1195 ring R=(0,a,b,c),(x,y),dp; 1196 short=0; 1197 // parametric polynom 1198 poly f=(b^21)*x^3*y+(c^21)*x*y^2+(c^2*bb)*x+(abc)*y; 1199 ideal p=c1; 1200 ideal q=ab; 1201 // normaform of f on V(p) \ V(q) 1202 pnormalf(f,p,q); 1211 { "EXAMPLE:"; echo = 2; 1212 1213 if(defined(R)){kill R;} 1214 ring R=(0,a,b,c),(x,y),dp; 1215 short=0; 1216 1217 // parametric polynom 1218 poly f=(b^21)*x^3*y+(c^21)*x*y^2+(c^2*bb)*x+(abc)*y; 1219 // ideals defining V(p)\V(q) 1220 ideal p=c1; 1221 ideal q=ab; 1222 1223 // pnormaform of f on V(p) \ V(q) 1224 pnormalf(f,p,q); 1203 1225 } 1204 1226 … … 1488 1510 ideals with P included in Q, representing the set 1489 1511 V(P) \ V(Q) = V(N) \ V(M) 1490 KEYWORDS: locally closed set; canon ical form1512 KEYWORDS: locally closed set; canoncial form 1491 1513 EXAMPLE: Crep; shows an example" 1492 1514 { … … 1518 1540 } 1519 1541 example 1520 { 1521 "EXAMPLE:"; echo = 2; 1542 { "EXAMPLE:"; echo = 2; 1522 1543 short=0; 1523 1544 if(defined(R)){kill R;} … … 1525 1546 ideal p=a*b; 1526 1547 ideal q=a,b2; 1548 1527 1549 // Crepresentation of V(p) \ V(q) 1528 1550 Crep(p,q); … … 1588 1610 Output: [Comp_1, .. , Comp_s ] where 1589 1611 Comp_i=[p_i,[p_i1,..,p_is_i]] 1590 KEYWORDS: locally closed set; canon ical form1612 KEYWORDS: locally closed set; canoncial form 1591 1613 EXAMPLE: Prep; shows an example" 1592 1614 { … … 1606 1628 } 1607 1629 example 1608 { 1609 "EXAMPLE:"; echo = 2; 1630 { "EXAMPLE:"; echo = 2; 1610 1631 short=0; 1611 1632 if(defined(R)){kill R;} … … 1613 1634 ideal p=a*b;; 1614 1635 ideal q=a,b1; 1636 1615 1637 // Prepresentation of V(p) \ V(q) 1616 1638 Prep(p,q); … … 1656 1678 } 1657 1679 } 1680 //"T_before="; prep; 1658 1681 if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));} 1682 //"T_Prep="; prep; 1683 //def Lout=CompleteA(prep,prep); 1684 //"T_Lout="; Lout; 1659 1685 return(prep); 1660 1686 } … … 1681 1707 P included in Q, representing the 1682 1708 set V(P) \ V(Q) 1683 KEYWORDS: locally closed set; canon ical form1709 KEYWORDS: locally closed set; canoncial form 1684 1710 EXAMPLE: PtoCrep; shows an example" 1685 1711 { … … 1701 1727 example 1702 1728 { 1703 "EXAMPLE:"; echo = 2; 1704 if(defined(R)){kill R;} 1705 ring R=0,(a,b,c),lp; 1706 short=0; 1707 ideal p=a*(a^2+b^2+c^225); 1708 ideal q=a*(a3),b4; 1709 // Crepresentaion of V(p) \ V(q) 1710 def Cr=Crep(p,q); 1711 Cr; 1712 // Crepresentation of V(p) \ V(q) 1713 def L=Prep(p,q); 1714 L; 1715 PtoCrep(L); 1729 echo = 2; 1730 //EXAMPLE: 1731 1732 if(defined(R)){kill R;} 1733 ring R=0,(a,b,c),lp; 1734 short=0; 1735 1736 ideal p=a*(a^2+b^2+c^225); 1737 ideal q=a*(a3),b4; 1738 1739 // Crepresentaion of V(p) \ V(q) 1740 def Cr=Crep(p,q); 1741 Cr; 1742 1743 // Prepresentation of V(p) \ V(q) 1744 def L=Prep(p,q); 1745 L; 1746 1747 PtoCrep(L); 1716 1748 } 1717 1749 … … 1734 1766 { 1735 1767 option(returnSB); 1768 //"T_Lp[i]="; Lp[i]; 1736 1769 N=Lp[i][1]; 1737 1770 Lb=Lp[i][2]; 1771 //"T_Lb="; Lb; 1738 1772 ida=intersect(ida,N); 1739 1773 for(j=1;j<=size(Lb);j++) … … 1947 1981 def BH=imap(RR,B2); 1948 1982 def LH=imap(RR,LL); 1983 //"T_BH="; BH; 1984 //"T_LH="; LH; 1949 1985 for (i=1;i<=size(BH);i++) 1950 1986 { … … 1958 1994 def RPH=DH+PH; 1959 1995 def GSH=KSW(BH,LH); 1996 //"T_GSH="; GSH; 1960 1997 //setglobalrings(); 1961 1998 // DEHOMOGENIZING THE RESULT … … 2011 2048 example 2012 2049 { 2013 "EXAMPLE:"; echo = 2; 2014 // Casas conjecture for degree 4: 2015 if(defined(R)){kill R;} 2016 ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp; 2017 short=0; 2018 ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0), 2019 x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1), 2020 x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0), 2021 x2^2+(2*a3)*x2+(a2), 2022 x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0), 2023 x3+(a3); 2024 cgsdr(F); 2050 echo = 2; 2051 // EXAMPLE: 2052 // Casas conjecture for degree 4: 2053 2054 // CasasAlvero conjecture states that on a field of characteristic 0, 2055 // if a polynomial of degree n in x has a common root whith each of its 2056 // n1 derivatives (not assumed to be the same), then it is of the form 2057 // P(x) = k(x + a)^n, i.e. the common roots must all be the same. 2058 2059 if(defined(R)){kill R;} 2060 ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp; 2061 short=0; 2062 2063 ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0), 2064 x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1), 2065 x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0), 2066 x2^2+(2*a3)*x2+(a2), 2067 x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0), 2068 x3+(a3); 2069 2070 cgsdr(F); 2025 2071 } 2026 2072 … … 2091 2137 P[i]=L[Q[i][1]][Q[i][2]]; 2092 2138 } 2139 //"T_P="; P; 2093 2140 // P is the list of the maximal components of the union 2094 2141 // with the corresponding initial holes. … … 2128 2175 int i; int j; int k; list H; list C; list T; 2129 2176 list L0; list P0; list P; list Q0; list Q; 2177 //"T_L="; L; 2130 2178 for (i=1;i<=size(L);i++) 2131 2179 { … … 2136 2184 } 2137 2185 } 2186 //"T_P0="; P0; 2138 2187 Q0=selectminideals(P0); 2188 //"T_Q0="; Q0; 2139 2189 for (i=1;i<=size(Q0);i++) 2140 2190 { … … 2142 2192 P[i]=L[Q0[i][1]];// [Q[i][2]]; 2143 2193 } 2194 //"T_P="; P; 2144 2195 // P is the list of the maximal components of the union 2145 2196 // with the corresponding initial holes. … … 2249 2300 static proc addpartfine(list H, list C0) 2250 2301 { 2302 //"T_H="; H; 2251 2303 int i; int j; int k; int te; intvec notQ; int l; list sel; 2252 2304 intvec jtesC; … … 2335 2387 static proc needcombine(list LCU,ideal N) 2336 2388 { 2389 //"Deciding if combine is needed";; 2337 2390 ideal BB; 2338 2391 int tes=0; int m=1; int j; int k; poly sp; … … 2509 2562 LCU[k][1]=B; 2510 2563 } 2564 //"Deciding if combine is needed"; 2511 2565 crep=PtoCrep(prep); 2512 2566 if(size(LCU)>1){tes=1;} … … 2789 2843 example 2790 2844 { 2791 "EXAMPLE:"; echo = 2; 2845 echo = 2; 2846 // EXAMPLE 1: 2792 2847 // Casas conjecture for degree 4: 2793 if(defined(R)){kill R;} 2794 ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp; 2795 short=0; 2796 ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0), 2797 x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1), 2798 x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0), 2799 x2^2+(2*a3)*x2+(a2), 2800 x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0), 2801 x3+(a3); 2802 grobcov(F); 2803 2804 // EXAMPLE: 2848 2849 // CasasAlvero conjecture states that on a field of characteristic 0, 2850 // if a polynomial of degree n in x has a common root whith each of its 2851 // n1 derivatives (not assumed to be the same), then it is of the form 2852 // P(x) = k(x + a)^n, i.e. the common roots must all be the same. 2853 2854 if(defined(R)){kill R;} 2855 ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp; 2856 short=0; 2857 2858 ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0), 2859 x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1), 2860 x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0), 2861 x2^2+(2*a3)*x2+(a2), 2862 x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0), 2863 x3+(a3); 2864 2865 grobcov(F); 2866 2867 // EXAMPLE 2 2805 2868 // M. Rychlik robot; 2806 2869 // Complexity and Applications of Parametric Algorithms of 2807 // 2870 // Computational Algebraic Geometry.; 2808 2871 // In: Dynamics of Algorithms, R. de la Llave, L. Petzold and J. Lorenz eds.; 2809 2872 // IMA Volumes in Mathematics and its Applications, 2810 // SpringerVerlag 118: 129 (2000).; 2811 // (18. Mathematical robotics: Problem 4, twoarm robot)." 2812 if (defined(R)){kill R;} 2813 ring R=(0,a,b,l2,l3),(c3,s3,c1,s1), dp; 2814 short=0; 2815 ideal S12=al3*c3l2*c1,bl3*s3l2*s1,c1^2+s1^21,c3^2+s3^21; 2816 S12; 2817 grobcov(S12); 2873 // SpringerVerlag 118: 129 (2000).; 2874 // (18. Mathematical robotics: Problem 4, twoarm robot). 2875 2876 if (defined(R)){kill R;} 2877 ring R=(0,a,b,l2,l3),(c3,s3,c1,s1), dp; 2878 short=0; 2879 2880 ideal S12=al3*c3l2*c1,bl3*s3l2*s1,c1^2+s1^21,c3^2+s3^21; 2881 S12; 2882 2883 grobcov(S12); 2818 2884 } 2819 2885 … … 2916 2982 example 2917 2983 { 2918 echo = 2; "EXAMPLE:"; 2919 if(defined(R)){kill R;} 2920 ring R=(0,a1,a2),(x),lp; 2921 short=0; 2922 poly f=(a1^24*a1+a2^2a2)*x+(a1^416*a1+a2^34*a2); 2923 ideal p=a1*a2; 2924 ideal q=a2^2a2,a1*a2,a1^24*a1; 2925 extendpoly(f,p,q);" "; 2926 // EXAMPLE; 2984 echo = 2; 2985 // EXAMPLE 1 2986 2987 if(defined(R)){kill R;} 2988 ring R=(0,a1,a2),(x),lp; 2989 short=0; 2990 2991 poly f=(a1^24*a1+a2^2a2)*x+(a1^416*a1+a2^34*a2); 2992 ideal p=a1*a2; 2993 ideal q=a2^2a2,a1*a2,a1^24*a1; 2994 2995 extendpoly(f,p,q); 2996 2997 // EXAMPLE 2 2998 2927 2999 if (defined(R)){kill R;} 2928 3000 ring R=(0,a0,b0,c0,a1,b1,c1,a2,b2,c2),(x), dp; 2929 3001 short=0; 3002 2930 3003 poly f=(b1*a2*c2c1*a2*b2)*x+(a1*c2^2+b1*b2*c2+c1*a2*c2c1*b2^2); 2931 3004 ideal p= … … 3090 3163 example 3091 3164 { 3092 "EXAMPLE:"; echo = 2; 3165 echo = 2; 3166 // EXAMPLE 3167 3093 3168 if(defined(R)){kill R;} 3094 3169 ring R=(0,a0,b0,c0,a1,b1,c1),(x), dp; 3095 3170 short=0; 3171 3096 3172 ideal S=a0*x^2+b0*x+c0, 3097 3173 a1*x^2+b1*x+c1; 3174 3098 3175 def GCS=grobcov(S,"rep",2); 3099 3176 // grobcov(S) with both P and C representations 3100 3177 GCS; 3101 3178 3102 3179 def FGC=extendGC(GCS,"rep",1); 3103 3180 // Full representation 3104 3181 FGC; 3105 3182 } … … 3243 3320 } 3244 3321 combpolys=reform(combpolys,numdens); 3322 //"T_combpolys="; combpolys; 3245 3323 //setring(RR); 3246 3324 //def combpolysT=imap(P,combpolys); … … 3376 3454 // where the w_l=(n_l1,..,n_ls) are intvec of length size(L), where 3377 3455 // n_lt fixes which element of (v_t1,..,v_tk_t) is to be 3378 // cho sen to form the tth (Q,P) for the lth element of the sheaf3456 // choosen to form the tth (Q,P) for the lth element of the sheaf 3379 3457 // representing the Iregular function. 3380 3458 // The selection is done to obtian the minimal number of elements … … 3956 4034 { 3957 4035 CG=KSWtocgsdr(CG); 4036 //"T_CG="; CG; 3958 4037 if( size(CG)>0) 3959 4038 { … … 4349 4428 { 4350 4429 t=T[i]; 4430 //"T_t"; t; 4351 4431 P=0; Q=1; 4352 4432 for(j=1;j<=n;j++) … … 4362 4442 } 4363 4443 Cb=Crep0(P,Q); 4444 //"T_Cb="; Cb; 4364 4445 if(size(Cb)!=0) 4365 4446 { … … 4450 4531 { 4451 4532 K=FirstLevel(B); 4533 //"T_K="; K; 4452 4534 L[size(L)+1]=K[1]; 4453 4535 B=K[2]; … … 4460 4542 example 4461 4543 { 4462 "EXAMPLE:"; echo = 2; 4544 echo = 2; 4545 // EXAMPLE: 4546 4463 4547 if(defined(R)){kill R;} 4464 4548 ring R=0,(x,y,z),lp; 4465 4549 short=0; 4550 4466 4551 ideal P1=(x^2+y^2+z^21); 4467 4552 ideal Q1=z,x^2+y^21; … … 4479 4564 def LL=ConsLevels(L); 4480 4565 LL; 4566 4481 4567 def LLL=Levels(LL); 4482 4568 LLL; … … 4517 4603 example 4518 4604 { 4519 "EXAMPLE:"; echo = 2; 4605 echo = 2; 4606 // EXAMPLE: 4607 4520 4608 if(defined(R)){kill R;} 4521 4609 ring R=0,(x,y,z),lp; 4522 4610 short=0; 4611 4523 4612 ideal P1=(x^2+y^2+z^21); 4524 4613 ideal Q1=z,x^2+y^21; … … 4536 4625 def LL=ConsLevels(L); 4537 4626 LL; 4627 4538 4628 def LLL=Levels(LL); 4539 4629 LLL; … … 4570 4660 if (size(B) mod 2==1){B[size(B)+1]=ideal(1);} 4571 4661 if (size(A) mod 2==1){A[size(A)+1]=ideal(1);} 4662 //"T_A=";A; 4663 // "T_B="; B; 4572 4664 n=size(A) div 2; 4573 4665 m=(size(B) div 2)1; … … 4588 4680 //string("T_j=",j); 4589 4681 ABdw=intersectpar(list(A[2*i],B[2*j+1])); 4682 //"T_ABdw="; ABdw; 4590 4683 ABup=A[2*i1]; 4684 //"T_ABup1="; ABup; 4591 4685 if(j>0) 4592 4686 { … … 4596 4690 } 4597 4691 } 4692 //"T_ABup2="; ABup; 4598 4693 ABupC=Crep(ABup,ideal(1)); 4694 //"T_ABupC="; ABupC; 4599 4695 ABup=ABupC[1]; 4696 //"T_ABup="; ABup; 4600 4697 if(ABup==1){t=0;} 4601 4698 //if(equalideals(ABup,ideal(1))){t=0;} … … 4603 4700 { 4604 4701 M=Crep(ABup,ABdw); 4702 //"T_M="; M; 4605 4703 //if(not(equalideals(M[1],ideal(1)))) {L[size(L)+1]=M;} 4606 4704 if(not(size(M)==0)) {L[size(L)+1]=M;} 4607 4705 } 4706 //"L="; L; 4608 4707 j++; 4609 4708 } … … 4614 4713 example 4615 4714 { 4616 "EXAMPLE:"; echo = 2; 4715 echo = 2; 4716 // EXAMPLE: 4717 4617 4718 if(defined(R)){kill R;} 4618 4719 ring R=(0,x,y,z,t),(x1,y1),lp; … … 4635 4736 def LL=DifConsLCSets(L1,L2); 4636 4737 LL; 4738 4637 4739 def LLL=ConsLevels(LL); 4638 4740 LLL; 4741 4639 4742 def LLLL=Levels(LLL); 4640 4743 LLLL; … … 4644 4747 4645 4748 //******************** Begin locus and envelop ****************************** 4749 4750 // Routines for defining different rings acting in the basic ring RR=Q[a,x][u,v], in lp order, where 4751 // a= parameters of the locus problem 4752 // x= tracer variables 4753 // u= auxiiary variables 4754 // v= mover variables 4755 4756 // Transforms the ringlist of Q[x_1,..,x_j] into the ringlist of Q[x_1,..,x_{n1},x_{m+1},..,x_j] 4757 // I.e., deletes the varibles x_n to x_m 4758 // To be used with the same order for all variables 4759 static proc Ldelnm(list LQx,int n,int m) 4760 { 4761 int i; 4762 int npara= m n+1; 4763 def RR=basering; 4764 def LR=LQx; 4765 int nt=size(LR[2]); 4766 def L1=LR[2]; 4767 for(i=n;i<= m;i++) {L1=delete(L1,n);} 4768 LR[2]=L1; 4769 intvec v; 4770 for(i=1;i<=ntnpara;i++){v[i]=1;} 4771 LR[3][1][2]=v; 4772 return(LR); 4773 } 4774 4775 // Transforms the ringlist of Q[a][x] into the ringlist of Q[a,x] 4776 // To be used with the same lp order 4777 proc La_xToLax(list La_x) 4778 { 4779 if(typeof(La_x[1])==typeof(0)){return(La_x);} 4780 list Lax=La_x[1]; 4781 if(Lax[1]=0){return(La_x);} 4782 list Va=Lax[2]; 4783 int na=size(Lax[2]); 4784 //"na=";na; 4785 list Vx=La_x[2]; 4786 list Vax=Va+Vx; 4787 int nx=size(Vx); 4788 //"nx="; nx; 4789 intvec vv; 4790 int i; 4791 for(i=1;i<=na+nx;i++){vv[i]=1;} 4792 Lax[2]=Vax; 4793 Lax[3][1][2]=vv; 4794 return(Lax); 4795 } 4796 4797 // Transforms the ringlist of Q[a,x] into the ringlist of Q[a][x] 4798 // To be used with the same lp order 4799 proc LaxToLa_x(list Lax,int nx) 4800 { 4801 //"T_Lax=",Lax; 4802 int nax=size(Lax[2]); 4803 int na=naxnx; 4804 if(na==0){return(Lax);} 4805 else 4806 { 4807 //string("T_ na=",na,", nx=",nx); 4808 list La_x; 4809 list Vax=Lax[2]; 4810 list Va; 4811 list Vx; 4812 int i; 4813 for(i=1;i<=na;i++){Va[size(Va)+1]=Vax[i];} 4814 intvec vva; 4815 for(i=1;i<=na;i++){vva[i]=1;} 4816 intvec vvx; 4817 for(i=1;i<=nx;i++){Vx[size(Vx)+1]=Vax[na+i];} 4818 for(i=1;i<=nx;i++){vvx[i]=1;} 4819 La_x[1]=Lax; 4820 La_x[1][2]=Va; 4821 La_x[1][3][1][2]=vva; 4822 La_x[2]=Vx; 4823 list lax3; 4824 lax3=Lax[3]; 4825 lax3[1][2]=vvx; 4826 La_x[3]=lax3; 4827 La_x[4]=Lax[4]; 4828 return(La_x); 4829 } 4830 } 4831 4832 // // Transforms the ringlist of Q[a,x] into the ringlist of Q[a][x] 4833 // // To be used with the same lp order 4834 // proc LaxToLa_x(list Lax,int nx) 4835 // { 4836 // //"T_Lax=",Lax; 4837 // int nax=size(Lax[2]); 4838 // int na=naxnx; 4839 // if(na==0){return(Lax);} 4840 // else 4841 // { 4842 // //string("T_ na=",na,", nx=",nx); 4843 // list La_x; 4844 // list Vax=Lax[2]; 4845 // list Va; 4846 // list Vx; 4847 // int i; 4848 // for(i=1;i<=na;i++){Va[size(Va)+1]=Vax[i];} 4849 // intvec vva; 4850 // for(i=1;i<=na;i++){vva[i]=1;} 4851 // intvec vvx; 4852 // for(i=1;i<=nx;i++){Vx[size(Vx)+1]=Vax[na+i];} 4853 // for(i=1;i<=nx;i++){vvx[i]=1;} 4854 // La_x[1]=Lax; 4855 // La_x[1][2]=Va; 4856 // La_x[1][3][1][2]=vva; 4857 // La_x[2]=Vx; 4858 // list lax3; 4859 // lax3=Lax[3]; 4860 // lax3[1][2]=vvx; 4861 // La_x[3]=lax3; 4862 // La_x[4]=Lax[4]; 4863 // return(La_x); 4864 // } 4865 // } 4866 4867 // proc Lax_uvToLa_v(list Lax_uv,int na, int nv) 4868 // { 4869 // int i; 4870 // def Lax=Lax_uv[1]; 4871 // int nax=size(Lax[2]); 4872 // int nuv=size(Lax_uv[2]); 4873 // def La=Lax; 4874 // list La2; 4875 // if(na===0){ 4876 // list Luv=Lax_uv; 4877 // Luv[1]=0; 4878 // list Lv= 4879 // ; 4880 // 4881 // } 4882 // for(i=1;i<=na;i++){} 4883 // } 4884 4885 // Transforms the set of ringlists of Q[a] and Q[x] into the ringlist of Q[a][x] 4886 // To be used with the same lp order 4887 proc LaLxToLa_x(list La,list Lx) 4888 { 4889 if(size(La)==0){return(Lx);} 4890 def L1=La; 4891 def L2=Lx; 4892 list L; 4893 L[1]=L1; 4894 L[2]=L2[2]; 4895 L[3]=L2[3]; 4896 L[4]=L2[4]; 4897 return(L); 4898 } 4899 4900 // Transforms the ringlist of Q[a,x] into the ring of Q[a] 4901 // proc LaxToLa(list Lax, int na) 4902 // { 4903 // int ntot=size(Lax[2]); 4904 // list La=Lax; 4905 // int i; 4906 // list V; 4907 // for(i=1;i<=ntotna;i++){V[i]=Lax[2][na+i];} 4908 // La[2]=V; 4909 // intvec vv; 4910 // for(i=1;i<=ntotna;i++){vv[i]=1;} 4911 // La[3][1][2]=vv; 4912 // return(La); 4913 // } 4914 4915 4916 // Transforms the set of ringlists of Q[a] and Q[x] into the ringlist of Q[a,x] 4917 // To be used with the same lp order 4918 proc LaLxToLax(list La,list Lx) 4919 { 4920 list L=LaLxToLa_x(La,Lx); 4921 list Lax=La_xToLax(L); 4922 return(Lax); 4923 } 4924 4925 // Transforms the ringlist of Q[a,x] into the ringlist of Q[a] 4926 // To be used with the same lp order 4927 proc LaxToLa(list Lax,int na) 4928 { 4929 list La; 4930 if(na==0){return(La);} 4931 else 4932 { 4933 La=Lax; 4934 list La2; 4935 for(i=1;i<=na;i++){La2[i]=Lax[2][i];} 4936 La[2]=La2; 4937 intvec va; 4938 for(i=1;i<=na;i++){va[i]=1;} 4939 La[3][1][2]=va; 4940 return(La); 4941 } 4942 } 4943 4944 // Transforms the ringlist of Q[a,x][u,v] into the ringlist of Q[a][v] 4945 // To be used with the same lp order 4946 proc Lax_uvToLa_v(list Lax_uv,int na, int nv) 4947 { 4948 //string("T_na=",na,"; nv=",nv); 4949 int i; 4950 list Lax=Lax_uv[1]; 4951 int nax=size(Lax[2]); 4952 int nuv=size(Lax_uv[2]); 4953 list Lv=Lax_uv; 4954 int nx=naxna; 4955 int nu=nuvnv; 4956 Lv[1]=0; 4957 list Lv2; 4958 intvec vv; 4959 for(i=1;i<=nv;i++){Lv2[i]=Lv[2][nu+i];} 4960 for(i=1;i<=nv;i++){vv[nu+i]=1;} 4961 Lv[2]=Lv2; 4962 Lv[3][1][2]=vv; 4963 if(na==0){return(Lv);} 4964 else 4965 { 4966 list La=Lax; 4967 list La2; 4968 intvec va; 4969 for(i=1;i<=na;i++){La2[i]=Lax[2][i];} 4970 for(i=1;i<=na; i++){va[i]=1;} 4971 La[2]=La2; 4972 La[3][1][2]=va; 4973 //"T_La="; La; 4974 //"T_Lv="; Lv; 4975 list La_v=LaLxToLa_x(La,Lv); 4976 return(La_v); 4977 } 4978 } 4979 4980 // Transforms the ringlist of Q[a,x] [u,v]into the ringlist of Q[x,u,a,v] 4981 // To be used with the same lp order 4982 proc Lax_uvToLxuav(list Lax_uv, int na, int nv) 4983 { 4984 //string("T_na=",na,"; nv=",nv); 4985 int i; 4986 //"T_Lax_uv="; Lax_uv; 4987 int nax=size(Lax_uv[1][2]); 4988 int nuv=size(Lax_uv[2]); 4989 int nx=naxna; 4990 int nu=nuvnv; 4991 //string("T_nax=",nax,"; nuv=",nuv,"; nx=",nx,"; nu=",nu); 4992 list Lxuav=Lax_uv[1]; 4993 list L2; 4994 for(i=1;i<=nx;i++){L2[i]=Lax_uv[1][2][na+i];} 4995 for(i=1;i<=nu;i++){L2[nx+i]=Lax_uv[2][i];} 4996 for(i=1;i<=na;i++){L2[nx+nu+i]=Lax_uv[1][2][i];} 4997 for(i=1;i<=nv;i++){L2[nx+nu+na+i]=Lax_uv[2][nu+i];} 4998 Lxuav[2]=L2; 4999 intvec vv; 5000 for(i=1;i<=nax+nuv;i++){vv[i]=1;} 5001 Lxuav[3][1][2]=vv; 5002 return(Lxuav); 5003 } 5004 5005 4646 5006 4647 5007 // indepparameters … … 4719 5079 attrib(NP,"IsSB",1); 4720 5080 int d=dim(std(NP)); 5081 //"T_d="; d; 4721 5082 if(d==0){te=0;} 4722 5083 setring(RR); … … 4724 5085 } 4725 5086 4726 // DimPar 4727 // Auxiliary routine to NorSing determining the dimension of a parametric ideal 4728 // Does not use @P and define it directly because is assumes that 4729 // variables and parameters have been inverted 4730 static proc DimPar(ideal E) 5087 // DimPar(E,nax,nx): 5088 // Auxilliary routine of locus2 determining the dimension of a component of the locus in 5089 // the ring Q[a][x] 5090 static proc DimPar(ideal E,nax,nx) 4731 5091 { 5092 //" ";"T_E in DimPar="; E; 4732 5093 def RRH=basering; 4733 5094 def RHx=ringlist(RRH); 4734 5095 def P=ring(RHx[1]); 4735 setring(P); 5096 list Lax=ringlist(P); 5097 //"Lax="; Lax; 5098 //int nax=size(Lax[2]); 5099 int na=naxnx; 5100 //string("T_na=",na,"; nx=",nx); 5101 list La_x=LaxToLa_x(Lax,nx); 5102 //"T_La_x="; La_x; 5103 def Qa_x=ring(La_x); 5104 setring(Qa_x); 5105 //setring(P); 4736 5106 def E2=std(imap(RRH,E)); 5107 //"T_E2 in DimPar="; E2; 4737 5108 attrib(E2,"IsSB",1); 4738 5109 def d=dim(E2); 5110 //string("T_d in DimPar=", d);" "; 4739 5111 setring RRH; 4740 5112 return(d); 4741 5113 } 4742 5114 4743 // DimComp 4744 // Auxiliary routine to locus2 determining the dimension of a parametric ideal 4745 // it is identical to DimPar but adds information about the character of the component 4746 static proc dimComp(ideal PA) 4747 { 4748 def RR=basering; 4749 list Rx=ringlist(RR); 4750 int n=size(Rx[1][2]); 4751 def P=ring(Rx[1]); 4752 setring(P); 4753 list Lout; 4754 def CP=imap(RR,PA); 4755 attrib(CP,"IsSB",1); 4756 int d=dim(std(CP)); 4757 if(d==n1){Lout[1]=d;Lout[2]="Degenerate"; } 4758 else {Lout[1]=d; Lout[2]="Accumulation";} 4759 setring RR; 4760 return(Lout); 5115 // DimComp 5116 // Auxilliary routine of locus2 determining the dimension of a parametric ideal 5117 // it is identical to DimPar but adds infromation about the character of the component 5118 static proc DimComp(ideal PA, int nax,int nx) 5119 { 5120 // def RR=basering; 5121 // list Rx=ringlist(RR); 5122 // int nax=size(Rx[1][2]); 5123 // int na=naxnx; 5124 // def P=ring(Rx[1]); 5125 // setring(P); 5126 // list Lout; 5127 // def CP=imap(RR,PA); 5128 // attrib(CP,"IsSB",1); 5129 // int d=dim(std(CP)); 5130 5131 list Lout; 5132 int d=DimPar(PA,nax,nx); 5133 if(d==nax1){Lout[1]=d;Lout[2]="Degenerate"; } 5134 else {Lout[1]=d; Lout[2]="Accumulation";} 5135 //"T_Lout="; Lout; 5136 setring RR; 5137 return(Lout); 4761 5138 } 4762 5139 … … 4782 5159 return(L3); 4783 5160 } 4784 4785 // NorSing4786 // Input:4787 // ideal B: the basis of a component of the grobcov4788 // ideal E: the top of the component (assumed to be of dimension > 0 (single equation)4789 // ideal N: the holes of the component4790 // Output:4791 // int d: the dimension of B on the variables (antiimage).4792 // if d>0 then the component is 'Normal'4793 // if d==0 then the component is 'Singular'4794 static proc NorSing(ideal B, ideal E, ideal N, list #)4795 {4796 int i; int j; int Fenv=0; int env; int dd;4797 list DD=#;4798 def RR=basering;4799 ideal vmov;4800 int version=0;4801 int nv=nvars(RR);4802 int moverdim=nv; //npars(RR);4803 if(nv<4){version=1;}4804 int d;4805 poly F;4806 for(i=1;i<=(size(DD) div 2);i++)4807 {4808 if(DD[2*i1]=="moverdim"){moverdim=DD[2*i];}4809 if(DD[2*i1]=="vmov"){vmov=DD[2*i];}4810 if(DD[2*i1]=="movdim"){moverdim=DD[2*i];}4811 if(DD[2*i1]=="version"){version=DD[2*i];}4812 if(DD[2*i1]=="family"){F=DD[2*i];}4813 }4814 if(F!=0){Fenv=1;}4815 list L0;4816 if(dimP0(E)==0){L0=2,"Normal";} // 2 es fals pero ha de ser >0 encara que sigui 04817 else4818 {4819 if(version==0)4820 {4821 // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part4822 // independent of parameters giving the variables with dimension 04823 dd=indepparameters(B);4824 if (dd==1){d=0; L0=d,string(B);} // ,determineF(B,F,E)4825 else{d=1; L0=2,"Normal";}4826 }4827 else4828 {4829 def RH=ringlist(RR);4830 def H=RH;4831 H[1]=0;4832 H[2]=RH[1][2]+RH[2];4833 int n=size(H[2]);4834 intvec ll;4835 for(i=1;i<=n;i++)4836 {4837 ll[i]=1;4838 }4839 H[3][1][1]="lp";4840 H[3][1][2]=ll;4841 def RRH=ring(H);4842 setring(RRH);4843 ideal BH=imap(RR,B);4844 ideal EH=imap(RR,E);4845 ideal NH=imap(RR,N);4846 if(Fenv==1){poly FH=imap(RR,F);}4847 for(i=1;i<=size(EH);i++){BH[size(BH)+1]=EH[i];}4848 BH=std(BH); // MOLT COSTOS!!!4849 ideal G;4850 ideal r; poly q;4851 for(i=1;i<=size(BH);i++)4852 {4853 r=factorize(BH[i],1);4854 q=1;4855 for(j=1;j<=size(r);j++)4856 {4857 if((pdivi(r[j],NH)[1] != 0) or (equalideals(ideal(NH),ideal(1))))4858 {4859 q=q*r[j];4860 }4861 }4862 if(q!=1){G[size(G)+1]=q;}4863 }4864 setring RR;4865 def GG=imap(RRH,G);4866 ideal GGG;4867 if(defined(L0)){kill L0; list L0;}4868 for(i=1;i<=size(GG);i++)4869 {4870 if(indepparameters(GG[i])){GGG[size(GGG)+1]=GG[i];}4871 }4872 GGG=std(GGG);4873 ideal GLM;4874 for(i=1;i<=size(GGG);i++)4875 {4876 GLM[i]=leadmonom(GGG[i]);4877 }4878 attrib(GLM,"IsSB",1);4879 d=dim(std(GLM));4880 string antiimi=string(GGG);4881 L0=d,antiimi;4882 if(d==0) // d<dim(V(E) \ V(N))?4883 {4884 " ";string("Antiimage of Special component = ", GGG);4885 }4886 else4887 {4888 L0[2]="Normal";4889 }4890 }4891 }4892 return(L0);4893 }4894 5161 4895 5162 static proc determineF(ideal A,poly F,ideal E) … … 4911 5178 def RRH=ring(H); 4912 5179 5180 //" ";string("Antiimage of Special component = ", GGG); 4913 5181 4914 5182 setring(RRH); … … 4919 5187 ideal M=std(AA+FH); 4920 5188 def rh=reduce(EH,M,5); 5189 //"T_AA="; AA; "T_FH="; FH; "T_EH="; EH; "T_rh="; rh; 4921 5190 if(rh==0){env=1;} else{env=0;} 4922 5191 setring RR; 4923 5192 //L0[3]=env; 5193 //"T_env="; env; 4924 5194 return(env); 4925 5195 } 4926 5196 4927 5197 4928 // locus0(G): Private routine used by locus (the public routine), that 4929 // builds the different components. 4930 // input: The output G of the grobcov (in generic representation, which is the default option) 4931 // Options: The algorithm allows the following options as pair of arguments: 4932 // "moverdim", d : by default moverdim is m=number of variables, 4933 // but it can be set to less values. 4934 // "version", v : There are two versions of the algorithm. ('version',1) is 4935 // a full algorithm that always distinguishes correctly between 'Normal' 4936 // and 'Special' components, whereas ('version',0) can decalre a component 4937 // as 'Normal' being really 'Special', but is more effective. By default ('version',1) 4938 // is used when the number of variables is less than 4 and 0 if not. 4939 // The user can force to use one or other version, but it is not recommended. 4940 // "system", ideal F: if the initial system is passed as an argument. 4941 // This is actually not used. 4942 // "comments", c: by default it is 0, but it can be set to 1. 4943 // Usually locus problems have mover coordinates, variables and tracer coordinates. 4944 // The mover coordinates are to be placed as the last variables, and by default, 4945 // its number is dim @P, but as option it can be set to other values. 4946 // output: 4947 // list, the canonical Prepresentation of the Normal and NonNormal locus: 4948 // The Normal locus has two kind of components: Normal and Special. 4949 // The Nonnormal locus has two kind of components: Accumulation and Degenerate. 4950 // This routine is complemented by locus that calls it in order to eliminate problems 4951 // with degenerate points of the mover. 4952 // The output components are given as 4953 // ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k) 4954 // The components are given in canonical Prepresentation of the subset. 4955 // If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level 4956 // gives the depth of the component. 4957 static proc locus0(list GG, list #) 4958 { 4959 int te=0; 4960 int t1=1; int t2=1; int i; 4961 def RR=basering; 4962 def Rx=ringlist(RR); 4963 def P=ring(Rx[1]); 4964 def RX=Rx; 4965 RX[1]=0; 4966 def D=ring(RX); 4967 def RP=D+P; 4968 // if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; } 4969 // setglobalrings(); 4970 //Options 4971 list DD=#; 4972 ideal vmov; 4973 int dimpar=size(Rx[1][2]); 4974 int dimvar=size(Rx[2]); 4975 int nv=nvars(RR); 4976 int moverdim=nv; //size(ringlist(RR)[1][2]); 4977 for(i=1;i<=nv;i++) 4978 { 4979 //vmov[size(vmov)+1]=var(i+nvmoverdim); 4980 vmov[size(vmov)+1]=var(i); 4981 } 4982 //string("T_moverdim=",moverdim,"; dimvar=",dimvar,"; dimpar=",dimpar,"; vmov=",vmov); 4983 int version=0; 4984 //if(nv<4){version=1;} 4985 int comment=0; 4986 ideal Fm; 4987 poly F; 4988 for(i=1;i<=(size(DD) div 2);i++) 4989 { 4990 if(DD[2*i1]=="moverdim"){moverdim=DD[2*i];} 4991 if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 4992 if(DD[2*i1]=="version"){version=DD[2*i];} 4993 if(DD[2*i1]=="system"){Fm=DD[2*i];} 4994 if(DD[2*i1]=="comment"){comment=DD[2*i];} 4995 if(DD[2*i1]=="family"){F=DD[2*i];} 4996 } 4997 list HHH; 4998 if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1) 4999 {return(HHH);} 5000 list G1; list G2; 5001 def G=GG; 5002 list Q1; list Q2; 5003 int d; int j; int k; 5004 t1=1; 5005 for(i=1;i<=size(G);i++) 5006 { 5007 attrib(G[i][1],"IsSB",1); 5008 d=dim(std(G[i][1])); 5009 if(d==0){G1[size(G1)+1]=G[i];} 5010 else 5011 { 5012 if(d>0){G2[size(G2)+1]=G[i];} 5013 } 5014 } 5015 if(size(G1)==0){t1=0;} // normal point components 5016 if(size(G2)==0){t2=0;} // nonnormal point components 5017 setring(RP); 5018 if(t1) 5019 { 5020 list G1RP=imap(RR,G1); 5021 } 5022 else {list G1RP;} 5023 list P1RP; 5024 ideal B; 5025 for(i=1;i<=size(G1RP);i++) 5026 { 5027 kill B; 5028 ideal B; 5029 for(k=1;k<=size(G1RP[i][3]);k++) 5030 { 5031 attrib(G1RP[i][3][k][1],"IsSB",1); 5032 G1RP[i][3][k][1]=std(G1RP[i][3][k][1]); 5033 for(j=1;j<=size(G1RP[i][2]);j++) 5034 { 5035 B[j]=reduce(G1RP[i][2][j],G1RP[i][3][k][1]); 5036 } 5037 B=std(B); 5038 P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B); 5039 } 5040 } //In P1RP the basis has been reduced wrt the top 5041 setring(RR); 5042 ideal h; 5043 list P1; 5044 if(t1) 5045 { 5046 P1=imap(RP,P1RP); 5047 for(i=1;i<=size(P1);i++) 5048 { 5049 for(j=1;j<=size(P1[i][3]);j++) 5050 { 5051 h=factorize(P1[i][3][j],1); 5052 P1[i][3][j]=h[1]; 5053 for(k=2;k<=size(h);k++) 5054 { 5055 P1[i][3][j]=P1[i][3][j]*h[k]; 5056 } 5057 } 5058 } 5059 } //In P1 factors in the basis independent of the parameters have been obtained 5060 ideal BB; int dd; list NS; 5061 // "T_AQUI_3"; 5062 for(i=1;i<=size(P1);i++) 5063 { 5064 NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD); 5065 dd=NS[1]; 5066 if(dd==0){P1[i][3]=NS;} //"Special"; //dd==0 5067 else{P1[i][3]="Normal";} 5068 } 5069 list P2; 5070 for(i=1;i<=size(G2);i++) 5071 { 5072 for(k=1;k<=size(G2[i][3]);k++) 5073 { 5074 P2[size(P2)+1]=list(G2[i][3][k][1],G2[i][3][k][2]); 5075 } 5076 } 5077 list l; 5078 for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1; 5079 for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2; 5080 5081 // if(defined(@P)){te=1; kill @P; kill @RP; kill @R;} 5082 // setglobalrings(); 5083 5084 setring(P); 5085 ideal J; 5086 if(t1==1) 5087 { 5088 def C1=imap(RR,P1); 5089 def L1=AddLocus(C1); 5090 } 5091 else{list C1; list L1; kill P1; list P1;} 5092 if(t2==1) 5093 { 5094 def C2=imap(RR,P2); 5095 def L2=AddLocus(C2); 5096 } 5097 else{list L2; list C2; kill P2; list P2;} 5098 for(i=1;i<=size(L2);i++) 5099 { 5100 J=std(L2[i][2]); 5101 d=dim(J); 5102 if(d+1==dimpar) 5103 { 5104 L2[i][4]=string("Degenerate",L2[i][4]); 5105 } 5106 else{L2[i][4]=string("Accumulation",L2[i][4]);} 5107 } 5108 5109 list LN; 5110 if(t1==1) 5111 { 5112 for(i=1;i<=size(L1);i++){LN[size(LN)+1]=L1[i];} 5113 } 5114 if(t2==1) 5115 { 5116 for(i=1;i<=size(L2);i++){LN[size(LN)+1]=L2[i];} 5117 } 5118 int tLN=1; 5119 if(size(LN)==0){tLN=0;} 5120 setring(RR); 5121 if(tLN) 5122 { 5123 def L=imap(P,LN); 5124 for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}} 5125 //list LL; 5126 for(i=1;i<=size(L);i++) 5127 { 5128 if(typeof(L[i][4])=="list") {L[i][4][1]="Special";} 5129 l[1]=L[i][2]; 5130 l[2]=L[i][3]; 5131 l[3]=L[i][4]; 5132 l[4]=L[i][5]; 5133 L[i]=l; 5134 } 5135 } 5136 else{list L;} 5137 // if(te==1){kill @P; kill @RP; kill @R;} 5138 return(L); 5139 } 5140 5141 // locus2(G): Private routine used by locus (the public routine), that 5142 // builds the different components. 5143 // input: The output G of the grobcov (in generic representation, which is the default option) 5144 // The ideal F defining the locus problem (G is the grobcov of F and has been 5198 5199 // locus2(G,F,moverdim,vmov,na): 5200 // Private routine used by locus (the public routine), that 5201 // builds the different component, and inputs for locus2 5202 // input: G= grobcov(S), already computed inside locus 5203 // F= the ideal defining the locus problem (G is the grobcov of F and has been 5204 // moverdim=number of mover variables 5205 // vmov= the ideal of the mover variables 5145 5206 // already determined by locus. 5146 // Options: The algorithm allows the following options as pair of arguments: 5147 // "moverdim", d : by default moverdim is the number of variables but it can 5148 // be set to minor values. 5149 // "version", v : There are two versions of the algorithm. ('version',1) is 5150 // a full algorithm that always distinguishes correctly between 'Normal' 5151 // and 'Special' components, whereas ('version',0) can decalre a component 5152 // as 'Normal' being really 'Special', but is more effective. By default ('version',1) 5153 // is used when the number of variables is less than 4 and 0 if not. 5154 // The user can force to use one or other version, but it is not recommended. 5155 // "system", ideal F: if the initial systrem is passed as an argument. This is actually 5156 // not used. 5157 // "comments", c: by default it is 0, but it can be set to 1. 5158 // Usually locus problems have mover coordinates, variables and tracer coordinates. 5159 // The mover coordinates are to be placed as the last variables, and by default, 5160 // its number is equal to the number of parameters of tracer variables, but as option 5161 // it can be set to other values. 5207 // na= number of parameteres of the locus problem (usually=0); 5208 // The arguments are determined by locus, and passed to locus2. 5162 5209 // output: 5163 5210 // list, the canonical Prepresentation of the Normal and NonNormal locus: … … 5171 5218 // If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level 5172 5219 // gives the depth of the component. 5173 // static proc locus2(list G, ideal F, list #) 5174 static proc locus2(list G, ideal F, list #) 5175 { 5176 int st=timer; 5177 list Snor; list Snonor; 5178 int d; int i; int j; int mt=0; 5179 def RR=basering; 5180 def Rx=ringlist(RR); 5181 def RP=ring(Rx[1]); 5182 int tax=1; 5183 list DD=#; 5184 list GG=G; 5185 int n=size(Rx[1][2]); 5186 int moverdim=n; 5187 for(i=1;i<=(size(DD) div 2);i++) 5188 { 5189 if(DD[2*i1]=="moverdim"){moverdim=DD[2*i]; mt=1; } 5190 if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 5191 // if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 5192 // if(DD[2*i1]=="version"){version=DD[2*i];} 5193 // if(DD[2*i1]=="comment"){comment=DD[2*i];} 5194 // if(DD[2*i1]=="grobcov"){GG=DD[2*i];} 5195 // if(DD[2*i1]=="antiimage"){tax=DD[2*i];} 5196 } 5197 if(mt==0){DD[size(DD)+1]="moverdim"; DD[size(DD)+1]=moverdim;} 5198 for(i=1;i<=size(GG);i++) 5199 { 5200 attrib(GG[i][1],"IsSB",1); 5201 GG[i][1]=std(GG[i][1]); 5202 d=dim(GG[i][1]); 5203 if(d==0) 5220 static proc locus2(list G, ideal F, int moverdim, ideal vmov, int na) 5221 { 5222 int st=timer; 5223 list Snor; list Snonor; 5224 int d; int i; int j; //int mt=0; 5225 def RR=basering; 5226 def Rx=ringlist(RR); 5227 def RP=ring(Rx[1]); 5228 def LP=ringlist(RP); 5229 int nax=size(LP[2]); 5230 int nx=naxna; 5231 int nv=moverdim; 5232 int tax=1; 5233 list GG=G; 5234 int n=size(Rx[1][2]); 5235 for(i=1;i<=size(GG);i++) 5204 5236 { 5205 for(j=1;j<=size(GG[i][3]);j++) 5206 { 5207 Snor[size(Snor)+1]=GG[i][3][j]; 5208 } 5209 } 5210 else 5211 { 5212 if(d>0) 5237 attrib(GG[i][1],"IsSB",1); 5238 GG[i][1]=std(GG[i][1]); 5239 d=dim(GG[i][1]); 5240 if(d==0) 5213 5241 { 5214 5242 for(j=1;j<=size(GG[i][3]);j++) 5215 5243 { 5216 Snonor[size(Snonor)+1]=GG[i][3][j]; 5217 } 5244 Snor[size(Snor)+1]=GG[i][3][j]; 5245 } 5246 } 5247 else 5248 { 5249 if(d>0) 5250 { 5251 for(j=1;j<=size(GG[i][3]);j++) 5252 { 5253 Snonor[size(Snonor)+1]=GG[i][3][j]; 5254 } 5255 } 5218 5256 } 5219 5257 } 5220 } 5221 int tnor=size(Snor); int tnonor=size(Snonor); 5222 setring RP; 5223 list SnorP; 5224 list SnonorP; 5225 if(tnor) 5226 { 5227 SnorP=imap(RR,Snor); 5228 st=timer; 5229 SnorP=LCUnionN(SnorP); 5230 } 5231 if(tnonor) 5232 { 5233 SnonorP=imap(RR,Snonor); 5234 SnonorP=LCUnionN(SnonorP); 5235 } 5236 setring RR; 5237 ideal C; list N; list BAC; list AI; 5238 list NSC; list DAC; 5239 list L; 5240 ideal B; 5241 int k; 5242 int j0; int k0; int te; 5243 poly kkk=1; 5244 ideal AI0; 5245 int dimP; 5246 5247 if(tnor) 5248 { 5249 Snor=imap(RP,SnorP); 5250 for(i=1;i<=size(Snor);i++) 5258 //"T_Snor="; Snor; 5259 //"T_Snonor="; Snonor; 5260 int tnor=size(Snor); int tnonor=size(Snonor); 5261 setring RP; 5262 list SnorP; 5263 list SnonorP; 5264 if(tnor) 5251 5265 { 5252 C=Snor[i][1]; 5253 N=Snor[i][2]; 5254 dimP=DimPar(C); 5255 AI=NS(F,G,C,N,DD); 5256 AI0=AI[2]; 5257 if(size(AI)==3){AI[2]="normal"; AI[3]=ideal(0);} 5258 else 5266 SnorP=imap(RR,Snor); 5267 st=timer; 5268 SnorP=LCUnionN(SnorP); 5269 } 5270 if(tnonor) 5271 { 5272 SnonorP=imap(RR,Snonor); 5273 SnonorP=LCUnionN(SnonorP); 5274 } 5275 //"T_SnorP after LCUnion="; SnorP; 5276 // "T_SnonorP after LCUnion="; SnonorP; 5277 setring RR; 5278 ideal C; list N; list BAC; list AI; 5279 list NSC; list DAC; 5280 list L; 5281 ideal B; 5282 int k; 5283 int j0; int k0; int te; 5284 poly kkk=1; 5285 ideal AI0; 5286 int dimP; 5287 5288 if(tnor) 5289 { 5290 Snor=imap(RP,SnorP); 5291 for(i=1;i<=size(Snor);i++) 5259 5292 { 5260 if(AI[1]>=dimP) 5261 { 5262 AI[2]="Normal"; 5263 if(tax>0){AI[3]=AI0;} 5264 } 5265 else 5266 { 5267 if(AI[1]<n1){AI[2]="Special"; AI[3]=AI0;} 5268 } 5293 C=Snor[i][1]; 5294 N=Snor[i][2]; 5295 dimP=DimPar(C,nax,nx); 5296 //"T_G="; G; 5297 AI=NS(F,G,C,N,moverdim,na,vmov,dimP); 5298 Snor[i][size(Snor[i])+1]=AI; 5269 5299 } 5270 Snor[i][size(Snor[i])+1]=AI; 5300 for(i=1;i<=size(Snor);i++) 5301 { 5302 L[size(L)+1]=Snor[i]; 5303 } 5304 } 5305 ideal AINN; 5306 if(tnonor) 5307 { 5308 Snonor=imap(RP,SnonorP); 5309 //"T_Snonor="; Snonor; 5310 //"T_G="; G; 5311 for(i=1;i<=size(Snonor);i++) 5312 { 5313 DAC=DimComp(Snonor[i][1],nax,nx); 5314 Snonor[i][size(Snonor[i])+1]=DAC; 5315 } 5316 for(i=1;i<=size(Snonor);i++) 5317 { 5318 L[size(L)+1]=Snonor[i]; 5319 } 5271 5320 } 5272 for(i=1;i<=size(Snor);i++) 5273 { 5274 L[size(L)+1]=Snor[i]; 5275 } 5276 } 5277 ideal AINN; 5278 if(tnonor) 5279 { 5280 Snonor=imap(RP,SnonorP); 5281 for(i=1;i<=size(Snonor);i++) 5282 { 5283 DAC=dimComp(Snonor[i][1]); 5284 Snonor[i][size(Snonor[i])+1]=DAC; 5285 // if(tax==1) 5286 // { 5287 // AINN=AISnonor(Snonor[i][1],GG,DD); 5288 // Snonor[i][3][3]=AINN; 5289 // // AQUI ES ON AFEGIM LA ANTIMATGE DE LES NONOR 5290 // } 5291 } 5292 for(i=1;i<=size(Snonor);i++) 5293 { 5294 L[size(L)+1]=Snonor[i]; 5295 } 5296 } 5297 return(L); 5298 } 5299 5300 // Auxiliary algorithm of locus2. 5321 return(L); 5322 } 5323 5324 // Auxilliary algorithm of locus2. 5301 5325 // The algorithm searches the basis corresponding to C, in the grobcov. 5302 5326 // It reduces the basis modulo the component. … … 5304 5328 // For each hole of the component 5305 5329 // it searches the segment where the hole is included 5306 // and select form its basis the polynomials5330 // and selects the polynomials from its basis 5307 5331 // only dependent on the variables. 5308 5332 // These polynomials are nonnull in an open set of … … 5320 5344 // taxonomy is not determined, and (n,'normal', 0) is returned as third 5321 5345 // element of the component. (n is the dimension of the space). 5322 static proc NS(ideal F,list G, ideal C, list N, list #)5346 static proc NS(ideal F,list G, ideal C, list N, int nv, int na,ideal vmov,int dimC) 5323 5347 { 5324 5348 // Initializing and defining rings 5325 int i; int j; int k; int te; int j0;int k0; int m; 5326 list DD=#; 5327 def RR=basering; 5349 int i; int j; int k; int te; int j0;int k0; int m; 5350 def RR=basering; 5351 def Lax_uv=ringlist(RR); 5352 Lax_uv[3][1][1]="lp"; 5353 int nax=size(Lax_uv[1][2]); 5354 int nuv=size(Lax_uv[2]); 5355 int nx=naxna; 5356 int nu=nuvnv; 5357 //"Lax_uv="; Lax_uv; 5358 def Lax=Lax_uv[1]; 5359 def Qax=ring(Lax); // ring Q[a,x] 5360 //"T_Lax="; Lax; 5361 def La_x=LaxToLa_x(Lax,nx); // ring Q[a][x] 5362 //"T_La_x="; La_x; 5363 def La_v=Lax_uvToLa_v(Lax_uv,na,nv); // ring Q[a][v] 5364 //"T_La_v="; La_v; 5365 def Lxuav=Lax_uvToLxuav(Lax_uv,na,nv); 5366 //"T_Lxuav="; Lxuav; 5367 5368 5369 // old rings 5328 5370 def Rx=ringlist(RR); 5329 5371 def Lx=Rx; 5330 def P=ring(Rx[1]); // ring Q[bx]5372 def P=ring(Rx[1]); // ring Q[a,x]] 5331 5373 Lx[1]=0; 5332 def D=ring(Lx); // ring Q[bu]5333 def PR0=P+D; 5374 def D=ring(Lx); // ring Q[u,v] 5375 def PR0=P+D; // ring Q[a,x,u,v] 5334 5376 def PRx=ringlist(PR0); 5335 5377 PRx[3][2][1]="lp"; 5336 def PR=ring(PRx); // ring Q[bx,bu] in lex order 5337 int n=size(Rx[1][2]); // number of parameters 5338 def nv=nvars(RR); // number of variables 5339 int moverdim=n; // number of mover variables 5340 ideal vmov; // mover variables (bw) 5341 for(i=1;i<=(size(DD) div 2);i++) 5342 { 5343 if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 5344 if(DD[2*i1]=="moverdim"){moverdim=DD[2*i];} 5345 // if(DD[2*i1]=="version"){version=DD[2*i];} 5346 // if(DD[2*i1]=="comment"){comment=DD[2*i];} 5347 // if(DD[2*i1]=="grobcov"){GG=DD[2*i];} 5348 // if(DD[2*i1]=="antiimage"){tax=DD[2*i];} 5349 } 5350 5351 //string("T_moverdim= ", moverdim); 5352 // vmov = the last vmov variables = mover variables 5353 for(i=1;i<=moverdim;i++) 5354 { 5355 vmov[size(vmov)+1]=var(i+nvmoverdim); 5356 } 5357 //string("T_nv=",nv," moverdim=",moverdim); 5358 // for(i=1;i<=nv;i++) 5359 // { 5360 // vmov[size(vmov)+1]=var(i); 5361 // } 5362 5378 // "T_PRx="; PRx; 5379 def PR=ring(PRx); // ring Q[a,x,u,v] in lex order 5380 // end of old rings 5381 5382 for(i=1;i<=nv;i++) 5383 { 5384 vmov[size(vmov)+1]=var(i+nuvnv); 5385 } 5386 //string("T_nv=",nv," moverdim=",nv); 5363 5387 int ddeg; int dp; 5364 5388 list LK; … … 5366 5390 for(i=1;i<=nv;i++){bu[i]=var(i);} 5367 5391 ideal mv; 5368 // ULL5369 // ideal of mover varibles5370 // for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);}5371 5392 for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);} 5372 5393 … … 5385 5406 if(te==1){"ERROR";} 5386 5407 def B=G[j0][2]; // Aixo aniria be per les nonor 5387 5388 // Searching the elements in Q[v_m] on basis different from B that are nul there 5408 //"T_B=G[j0][2]="; B; 5409 //string("T_k0=",k0," G[",j0,"]="); G[j0]; 5410 5411 // Searching the elements in Q[v_m] on basis differents from B that are nul there 5389 5412 // and cannot become 0 on the antiimage of B. They are placed on NoNul 5390 5413 list NoNul; … … 5440 5463 } 5441 5464 5442 5443 5444 5445 5446 5465 // Adding F to B 5466 for(i=1;i<=size(F);i++) 5467 { 5468 B[size(B)+1]=F[i]; 5469 } 5447 5470 5448 5471 def E=NoNul; … … 5457 5480 setring(RR); 5458 5481 // if(n+nv>10 or ddeg>=10){LK=n,ideal(0),"normal"; return(LK);} 5459 if(ddeg>=10){LK=n ,ideal(0),"normal"; return(LK);} // 8 instead of 10 ?5482 if(ddeg>=10){LK=nv,ideal(0),"normal"; return(LK);} // 8 instead of 10 ? 5460 5483 5461 5484 // Reducing basis B modulo C in the ring PR lex(x,u) 5462 setring(PR); 5463 def buPR=imap(RR,bu); 5464 def EE=imap(RR,E); 5465 ideal BR; 5466 def BB=imap(RR,B); 5485 // setring(PR); 5486 5487 def Qxuav=ring(Lxuav); 5488 setring(Qxuav); 5489 def BR=imap(RR,B); 5490 ideal vamov; 5491 for(i=nx+nu+1;i<=nax+nuv;i++){vamov[size(vamov)+1]=var(i);} 5492 5493 //BR=std(BR); 5467 5494 def CC=imap(RR,C); 5495 for(i=1;i<=size(CC);i++){BR[size(BR)+1]=CC[i];} 5496 BR=std(BR); 5497 // for(i=1;i<=size(CC);i++){BR[size(BR)+1]=CC[i];} 5468 5498 attrib(CC,"IsSB",1); 5469 CC=std(CC); 5470 for(i=1;i<=size(BB);i++) 5471 { 5472 BR[i]=reduce(BB[i],CC); 5473 } 5474 5475 // Computing the GB(B) in the ring PR lex(x,u) 5476 BR=std(BR); 5477 ideal BRP; 5478 5479 // Eliminating the polynomials in B containing variables v_t 5480 for(i=1;i<=size(BR);i++) 5481 { 5482 if(subset(variables(BR[i]),buPR )){BRP[size(BRP)+1]=BR[i];} 5483 } 5484 BR=BRP; 5485 5486 // Factoring and eliminating the factors in N 5487 list H; 5488 ideal BP; 5489 poly f; 5490 int tj; 5491 //ideal vv=imap(RR,bu); 5492 for(i=1;i<=size(BR);i++) 5493 { 5494 f=1; 5495 H=factorize(BR[i]); 5496 if((subset(variables(H[1][2]),buPR)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];} 5497 else 5498 { 5499 for(j=1;j<=size(H[1]);j++) 5500 { 5501 tj=subset(H[1][j],EE); 5502 if((subset(variables(H[1][j]),buPR)) and (tj==0)){f=f*H[1][j];} 5503 } 5504 if(leadcoef(f)!=f){BP[size(BP)+1]=f;} 5505 } 5506 } 5507 5508 // Determining the GB basis of B in Q[u] 5509 BP=std(BP); 5510 // Determining the dimension of the antiimage 5511 // in the ring D 5512 setring D; 5513 def KK=imap(PR,BP); 5514 KK=std(KK); 5515 ideal KKM; 5516 ideal vmovM=imap(RR,vmov); 5517 5518 // selecting the element in Q[w] 5519 for(i=1;i<=size(KK);i++) 5520 { 5521 if(subset(variables(KK[i]),vmovM)) 5522 { 5523 KKM[size(KKM)+1]=KK[i]; 5524 } 5525 } 5526 5527 // Determining the dimensions 5528 list AIM=dimM(KKM,nv,moverdim); 5529 //int d2=AIM[1]; 5530 def AAIM=AIM[2]; 5499 ideal AIM; 5500 // "T_BR="; BR; 5501 for(i=1;i<=size(BR);i++){if(subset(variables(BR[i]),vamov)){AIM[size(AIM)+1]=BR[i];}} 5502 //"T_AIM="; AIM; 5503 5504 list La_v0=imap(RR,La_v); 5505 def Qa_v=ring(La_v0); 5506 setring Qa_v; 5507 def AIMa_v=imap(Qxuav,AIM); 5508 AIMa_v=std(AIMa_v); 5509 int dimAIM=dim(AIMa_v); 5510 //"T_AIMa_v="; AIMa_v; 5511 //string("T_dimAIM=",dimAIM); 5531 5512 setring(RR); 5532 //def BBR=imap(D,KK); 5533 def LKK=imap(D,AIM); 5534 //LKK=d2,string(BBR); 5535 return(LKK); 5536 } 5537 5538 // Auxiliary algorithm of locus2. 5539 // The algorithm searches the basis corresponding to C, in the grobcov. 5540 // It reduces the basis modulo the component. 5541 // The result is the reduced basis BR. 5542 // for each hole of the component 5543 // it searches the segment where the hole is included 5544 // and select form its basis the polynomials 5545 // only dependent on the variables. 5546 // These polynomials are nonnull in an open set of 5547 // the component, and are included in the list NoNul of nonnull factors 5548 // input: C: The top of the the segment of a nonnormal locus 5549 // G the grobcov of F 5550 // output: B: the antiimage of the input segment 5551 static proc AISnonor(ideal C,list G, list #) 5552 { 5553 // Initializing and defining rings 5554 int i; int j; int k; int te; int j0;int k0; int m; 5555 list DD=#; 5556 def RR=basering; 5557 def Rx=ringlist(RR); 5558 def Lx=Rx; 5559 def P=ring(Rx[1]); // ring Q[bx] 5560 Lx[1]=0; 5561 def D=ring(Lx); // ring Q[bu] 5562 def PR0=P+D; 5563 def PRx=ringlist(PR0); 5564 PRx[3][2][1]="lp"; 5565 def PR=ring(PRx); // ring Q[bx,bu] in lex order 5566 int n=size(Rx[1][2]); // number of parameters 5567 def nv=nvars(RR); // number of variables 5568 int moverdim=n; // number of parameters 5569 ideal vmov; // mover variables (bw) 5570 5571 for(i=1;i<=(size(DD) div 2);i++) 5572 { 5573 if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 5574 if(DD[2*i1]=="moverdim"){moverdim=DD[2*i];} 5575 // if(DD[2*i1]=="version"){version=DD[2*i];} 5576 // if(DD[2*i1]=="comment"){comment=DD[2*i];} 5577 // if(DD[2*i1]=="grobcov"){GG=DD[2*i];} 5578 // if(DD[2*i1]=="antiimage"){tax=DD[2*i];} 5579 } 5580 for(i=1;i<=moverdim;i++) 5581 { 5582 vmov[size(vmov)+1]=var(i+nvmoverdim); 5583 } 5584 5585 int ddeg; int dp; 5586 list LK; 5587 ideal bu; // ideal of all variables (u) 5588 for(i=1;i<=nv;i++){bu[i]=var(i);} 5589 ideal mv; // ideal of mover varibles 5590 for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);} 5591 5592 // Searching the basis associated to Nonor 5593 j=2; te=1; 5594 while((te) and (j<=size(G))) 5595 { 5596 k=1; 5597 while((te) and (k<=size(G[j][3]))) 5598 { 5599 if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;} 5600 k++; 5601 } 5602 j++; 5603 } 5604 if(te==1){"ERROR";} 5605 def B=G[j0][2]; 5606 return(B); 5607 } 5608 5609 // NOT USED 5610 // Auxiliary algorithm of locus2. 5611 // The algorithm searches the basis corresponding to C, in the grobcov. 5612 // It reduces the basis modulo the component. 5613 // The result is the reduced basis BR. 5614 // for each hole of the component 5615 // it searches the segment where the hole is included 5616 // and select form its basis the polynomials 5617 // only dependent on the variables. 5618 // These polynomials are nonnull in an open set of 5619 // the component, and are included in the list NoNul of nonnull factors 5620 // input: F: the ideal of the locus problem 5621 // G the grobcov of F 5622 // C the top of a component of normal points 5623 // N the holes of the component 5624 // output: (d,tax,a) 5625 // where d is the dimension of the antiimage 5626 // a is the antiimage of the component and 5627 // tax is the taxonomy \"Normal\" if d is equal to the dimension of C 5628 // and \"Special\" if it is smaller. 5629 static proc NS5(ideal F,list G, ideal C, list N) 5630 { 5631 // Initializing and defining rings 5632 int i; int j; int k; int te; int j0;int k0; int m; 5633 def RR=basering; 5634 def Rx=ringlist(RR); 5635 def Lx=Rx; 5636 def P=ring(Rx[1]); 5637 Lx[1]=0; 5638 def D=ring(Lx); 5639 def RP=D+P; 5640 def PR0=P+D; 5641 def PRx=ringlist(PR0); 5642 PRx[3][2][1]="lp"; 5643 def PR=ring(PRx); 5644 int n=size(Rx[1][2]); // dimension of the tracer space 5645 def nv=nvars(RR); // number of variables 5646 ideal v; for(i=1;i<=nv;i++){v[i]=var(i);} // variables 5647 int nm; // number of mover variables 5648 ideal vmov; // mover variables 5649 if(nv>=n){nm=n;} 5650 else{nm=nv;} 5651 for(i=1;i<=nm;i++) 5652 { 5653 vmov[size(vmov)+1]=var(i+nvnm); 5654 } 5655 5656 // Searching the basis associated to C 5657 j=2; te=1; 5658 while((te) and (j<=size(G))) 5659 { 5660 k=1; 5661 while((te) and (k<=size(G[j][3]))) 5662 { 5663 if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;} 5664 k++; 5665 } 5666 j++; 5667 } 5668 if(te==1){"ERROR";} 5669 def B=G[j0][2]; // basis associated to C 5670 5671 // Searching the elements in Q[vmov] on bases different from B that are nul there 5672 // and cannot become 0 on the antiimage of C. They are placed on NoNul 5673 list NoNul; // elements in Q[vmov] on bases different from B 5674 ideal BNoNul; 5675 ideal covertop; // basis of the segment whose top is a hole of our segment 5676 int te1; 5677 for(i=1;i<=size(N);i++) 5678 { 5679 j=2; te=1; 5680 while(te and j<=size(G)) 5681 { 5682 if(j!=j0) 5683 { 5684 k=1; 5685 while(te and k<=size(G[j][3])) 5686 { 5687 covertop=G[j][3][k][1]; 5688 if(equalideals(covertop,N[i])) 5689 { 5690 te=0; te1=1; BNoNul=G[j][2]; 5691 } 5692 else 5693 { 5694 if(redPbasis(covertop,N[i])) 5695 { 5696 te=0; te1=1; m=1; 5697 while( te1 and m<=size(G[j][3][k][2]) ) 5698 { 5699 if(equalideals(G[j][3][k][2][m] ,N[i] )==1){te1=0;} 5700 m++; 5701 } 5702 } 5703 if(te1==1){ BNoNul=G[j][2];} 5704 } 5705 k++; 5706 } 5707 5708 if((te==0) and (te1==1)) 5709 { 5710 // Selecting the elements independent of the parameters, 5711 // They will be non null on the segment 5712 for(m=1;m<=size(BNoNul);m++) 5713 { 5714 if(indepparameterspoly(BNoNul[m])) 5715 { 5716 NoNul[size(NoNul)+1]=BNoNul[m]; 5717 } 5718 } 5719 } 5720 } 5721 j++; 5722 } 5723 } 5724 def NN=NoNul; 5725 poly kkk=1; 5726 if(size(NN)==0){NN[1]=kkk;} 5727 5728 // Union of F and B 5729 for(i=1;i<=size(F);i++) 5730 { 5731 B[size(B)+1]=F[i]; 5732 } 5733 5734 // Reducing basis B modulo C in the ring PR 5735 setring(PR); 5736 def vv=imap(RR,v); 5737 ideal BBR; 5738 def BB=imap(RR,B); 5739 def CC=imap(RR,C); 5740 attrib(CC,"IsSB",1); 5741 CC=std(CC); 5742 for(i=1;i<=size(BB);i++) 5743 { 5744 BBR[i]=reduce(BB[i],CC); 5745 } 5746 5747 // Selecting the elements of the ideal in Q[v] 5748 // setring RP; 5749 ideal Bv; 5750 // def Bv=imap(PR,BBR); 5751 //ideal vvv=imap(RR,v); 5752 // Bv=std(Bv); 5753 ideal BR; 5754 for(i=1;i<=size(BBR);i++) 5755 { 5756 if(subset(variables(BBR[i]),vv)){BR[size(BR)+1]=BBR[i];} 5757 } 5758 // setring PR; 5759 // def BR=imap(RP,BRv); 5760 def NNN=imap(RR,NN); 5761 5762 // Factoring and eliminating the factors in N 5763 list H; 5764 ideal BP; 5765 poly f; 5766 int tj; 5767 // ideal vv=imap(RR,v); 5768 for(i=1;i<=size(BR);i++) 5769 { 5770 f=1; 5771 H=factorize(BR[i]); 5772 if((subset(variables(H[1][2]),vv)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];} 5773 else 5774 { 5775 for(j=1;j<=size(H[1]);j++) 5776 { 5777 tj=subset(H[1][j],NNN); 5778 if((subset(variables(H[1][j]),vv)) and (tj==0)){f=f*H[1][j];} 5779 } 5780 if(leadcoef(f)!=f){BP[size(BP)+1]=f;} 5781 } 5782 } 5783 5784 // Obtaining Am, the antiimage of C on vmov 5785 setring D; 5786 ideal A=imap(PR,BP); 5787 ideal vm=imap(RR,vmov); 5788 ideal Am; 5789 for(i=1;i<=size(A);i++) 5790 { 5791 if (subset(variables(A[i]),vm)){Am[size(Am)+1]=A[i];} 5792 } 5793 5794 list AIM=dimM(Am,nv,nm); 5795 setring(RR); 5796 def LAM=imap(D,AIM); 5797 return(LAM); 5798 } 5799 5800 static proc dimM(ideal KKM, int nv, int nm) 5513 def AIMRR=imap(Qa_v,AIMa_v); 5514 string TaxComp; 5515 if(dimAIM==dimC){TaxComp="Normal";} 5516 else{TaxComp="Special";} 5517 list NSA=dimAIM,TaxComp,AIMRR; 5518 return(NSA); 5519 } 5520 5521 static proc DimM(ideal KKM, int na, int nv) 5801 5522 { 5802 5523 def RR=basering; … … 5804 5525 int i; 5805 5526 def Rx=ringlist(RR); 5806 for(i=1;i<=nm;i++) 5527 5528 for(i=1;i<=nv;i++) 5807 5529 { 5808 5530 L[i]=Rx[2][nvnm+i]; … … 5813 5535 Rx[3][1][2]=iv; 5814 5536 def DM=ring(Rx); 5537 //"Rx="; Rx; 5815 5538 setring(DM); 5816 5539 ideal KKMD=imap(RR,KKM); … … 5821 5544 def KAIM=imap(DM,KKMD); 5822 5545 list LAIM=d,KAIM; 5546 // "T_LAIM="; LAIM; 5823 5547 return(LAIM); 5824 5548 } … … 5858 5582 ideal vpar; 5859 5583 ideal vvar; 5584 //"T_n="; n; 5585 //"T_nv="; nv; 5860 5586 for(i=1;i<=n;i++){vpar[size(vpar)+1]=par(i);} 5861 5587 for(i=1;i<=nv;i++){vvar[size(vvar)+1]=var(i);} 5588 //string("T_vpar = ", vpar," vvar = ",vvar); 5862 5589 def P=ring(Rx[1]); 5863 5590 Rx[1]=0; … … 5870 5597 //string("T_vvpar = ",vvpar); 5871 5598 ideal B=std(FF); 5599 //"T_B="; B; 5872 5600 ideal Bel; 5601 //"T_vvpar="; vvpar; 5873 5602 for(i=1;i<=size(B);i++) 5874 5603 { 5875 5604 if(subset(variables(B[i]),vvpar)) {Bel[size(Bel)+1]=B[i];} 5876 5605 } 5606 //"T_Bel="; Bel; 5877 5607 list H; 5878 5608 list FH; … … 5904 5634 // locus(F): Special routine for determining the locus of points 5905 5635 // of geometrical constructions. 5906 // input: The ideal of the locus equations 5636 // input: The ideal of the locus equations defined in the 5637 // ring Q[a1,..,ap,x1,..xn][u1,..um,v1,..vn] 5907 5638 // output: 5908 5639 // The output components are given as … … 5911 5642 // The third element tax is: 5912 5643 // for normal point components, tax=(d,taxonomy,antiimage) 5913 // being d=dimension of the antiimage on the mover var aibles,5644 // being d=dimension of the antiimage on the mover variables, 5914 5645 // taxonomy='Normal' or 'Special', and 5915 5646 // antiimage=ideal of the antiimage over the mover variables 5916 // which are taken to be the last n variables.5647 // which by default are taken to be the last n variables. 5917 5648 // for nonnormal point components, tax =(d,taxonomy) 5918 5649 // being d=dimension of the component and … … 5922 5653 // Normal component: 5923 5654 //  each point in the component has 0dimensional antiimage. 5924 //  the antiimage in the mover coordinates is equal to the dimension of the component 5655 //  the antiimage in the mover coordinates is equal to the dimension of the component. 5925 5656 // Special component: 5926 5657 //  each point in the component has 0dimensional antiimage. 5927 //  the antiimage in the mover coordinates is smaller than the dimension of the component5658 //  the antiimage on the mover variables is smaller than the dimension of the component. 5928 5659 // The nonnormal locus has two kind of components: Accumulation and Degenerate. 5929 5660 // Accumulation points: … … 5936 5667 // taxonomy is not determined, and (n,'normal', 0) is returned as third 5937 5668 // element of the component. (n is the dimension of the space). 5938 5939 5669 proc locus(ideal F, list #) 5940 "USAGE: locus(ideal F [,options])5670 "USAGE: locus(ideal F [,options]) 5941 5671 Special routine for determining the locus of points of 5942 5672 a geometrical construction. 5943 INPUT: The input ideal must be the set equations 5944 defining the locus. Calling sequence: 5945 locus(F[,options]); 5946 The input ring must be a parametrical ideal 5947 in Q[x1,..,xn][u1,..,um], 5948 (x=tracer variables, u=mover variables). 5949 The tracer variables are x1,..xn = dim(xspace). 5950 By default, the mover variables are the last n uvariables. 5673 INPUT: The input ideal must be the ideal of the set equations 5674 defining the locus, defined in the ring 5675 ring Q(0,a1,..,ap,x1,..xn)(u1,..um,v1,..vn),lp; 5676 Calling sequence: 5677 locus(F [,options]); 5678 a=fixed parameters,x=tracer variables, u=auxiliary variables, v=mover variables. 5679 The parameters a are optative. If they are used, then the option \"numpar\=,np 5680 must be declared, being np the number of fixed parameters. 5681 The tracer variables are x1,..xn, where n is the dimension of the space. 5682 By default, the mover variables are the last n variables. 5951 5683 Its number can be forced by the user to the last 5952 5684 k variables by adding the option \"moverdim\",k. 5953 5685 Nevertheless, this option is recommended only 5954 5686 to experiment, and can provide incorrect taxonomies. 5687 The remaining variables are auxiliary. 5955 5688 OPTIONS: An option is a pair of arguments: string, integer. 5956 5689 To modify the default options, pairs of arguments … … 5959 5692 pair of arguments: 5960 5693 5961 \"moverdim\", k to restrict the movervariables of the 5962 antiimage to the last k variables. 5963 By defaulat k is equal to the last n uvariables, 5694 \"numpar\", np in order to consider the first np parameters of the ring 5695 to be fixed parameters of the locus, being the tracer variables 5696 the remaining parameters. 5697 To be used for a paramteric locus. (New in release N12). 5698 5699 \"moverdim\", k to force the movervariables to be the last 5700 k variables. This determines the antiimage and its dimension. 5701 By defaulat k is equal to the last n variables, 5964 5702 We can experiment with a different value, 5965 5703 but this can produce an error in the character 5966 5704 \"Normal\" or \"Special\" of a locus component. 5967 5705 5968 \"grobcov\", L, where Lis the list of a previous computed grobcov(F).5706 \"grobcov\", G, where G is the list of a previous computed grobcov(F). 5969 5707 It is to be used when we modify externally the grobcov, 5970 5708 for example to obtain the real grobcov. 5971 5709 5972 \"comments\", c: by default it is 0, but it can be set 5973 to 1. 5710 \"comments\", c: by default it is 0, but it can be set to 1. 5974 5711 RETURN: The output is a list of the components: 5975 5712 ((p1,(p11,..p1s_1),tax_1), .., (pk,(pk1,..pks_k),tax_k) … … 5998 5735  each point in the component has 0dimensional 5999 5736 antiimage. 6000  the antiimage in the mover coordinates is smaller6001 than the dimension of the component5737  the antiimage in the mover coordinates has dimension 5738 smaller than the dimension of the component 6002 5739 The nonnormal locus has two kind of components: 6003 5740 Accumulation and Degenerate. … … 6013 5750 then the taxonomy is not determined, and (n,'normal', 0) 6014 5751 is returned as third element of the component. (n is the 6015 dimension of the space).5752 dimension of the tracer space). 6016 5753 6017 5754 Given a parametric ideal F representing the system F … … 6031 5768 EXAMPLE: locus; shows an example" 6032 5769 { 6033 int tes=0; int i; int m; int mm; int n;5770 int tes=0; int i; int m; int mm; // int n; 6034 5771 def RR=basering; 6035 5772 list GG; 6036 5773 //Options 6037 5774 list DD=#; 6038 ideal vmov; 6039 int nv=nvars(RR); // number of variables 6040 // int moverdim=nv; //size(ringlist(RR)[1][2]); // number of parameters 6041 int moverdim=size(ringlist(RR)[1][2]); // number of parameters 6042 if(moverdim>nv){moverdim=nv;} 6043 // for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);} 6044 int version=2; 6045 // if(nv<4){version=2;} 5775 int nax=size(ringlist(RR)[1][2]); // number of parameters + tracer variables 5776 int nuv=nvars(RR); // number of variables 5777 int na=0; int nx=nax; 5778 int moverdim=nx; // number of tracer variables 5779 if(moverdim>nuv){moverdim=nuv;} 5780 // int version=2; 6046 5781 int comment=0; 6047 5782 int tax=1; … … 6049 5784 for(i=1;i<=(size(DD) div 2);i++) 6050 5785 { 6051 if(DD[2*i1]=="moverdim"){moverdim=DD[2*i];} 6052 if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 6053 // if(DD[2*i1]=="version"){version=DD[2*i];} 5786 if(DD[2*i1]=="numpar"){na=DD[2*i]; nx=naxna; moverdim=nx;} 6054 5787 if(DD[2*i1]=="comment"){comment=DD[2*i];} 6055 5788 if(DD[2*i1]=="grobcov"){GG=DD[2*i];} 6056 5789 } 6057 // if(version != 2){string("versions 0 and 1 are not available"); version=2;} 6058 for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nvmoverdim);} 5790 for(i=1;i<=(size(DD) div 2);i++) 5791 { 5792 if(DD[2*i1]=="moverdim"){moverdim=DD[2*i];} 5793 } 5794 int nv=moverdim; 5795 if(moverdim>nuv){moverdim=nuv;} 5796 5797 ideal vmov; 5798 //string("T_nuv=",nuv,"; moverdim=",moverdim); 5799 for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nuvmoverdim);} 6059 5800 if(size(GG)==0){GG=grobcov(F);} 6060 5801 int j; int k; int te; … … 6063 5804 list nGP; 6064 5805 if (equalideals(B0,ideal(1)) ) 6065 {return(locus2(GG,F,DD));} 6066 // if(version==2){return(locus2(GG,F,DD));} 6067 // else{return(locus0(GG,DD));} 6068 // } 5806 {return(locus2(GG,F,moverdim,vmov,na));} 6069 5807 else 6070 5808 { 6071 n=nvars(RR);6072 5809 ideal vB; 6073 5810 ideal N; … … 6078 5815 attrib(N,"IsSB",1); 6079 5816 N=std(N); 6080 if((size(N))>=2)5817 if((size(N))>=2) 6081 5818 { 6082 5819 //def dN=dim(N); … … 6120 5857 nGP=GG; 6121 5858 " ";string("Unavoidable ",moverdim,"dimensional locus"); 6122 //string("Try option 'vmov',ideal(of mover variables) to avoid some point of the mover");6123 //" ";"Elements of the basis of the generic segment in mover variables="; N;" ";6124 5859 list L; return(L); 6125 5860 } 6126 5861 } 6127 if(comment>0){"Input for locus2 GB="; nGP; "input for locus F="; F;} 6128 if(version==2) 6129 { 6130 def LL=locus2(nGP,F,DD); 6131 } 6132 else{def LL=locus0(nGP,DD);} 6133 // kill @RP; kill @P; kill @R; 5862 5863 // if(comment>0){"Input for locus2 GB="; nGP; "input for locus F="; F;} 5864 // if(version==2) 5865 // { 5866 // "T_nGP enter for locus2="; nGP; 5867 // def LL=locus2(nGP,F,moverdim,vmov,na); 5868 // } 5869 // else{ def LL=locus0(nGP,moverdim,vmov); } 5870 5871 def LL=locus2(nGP,F,moverdim,vmov,na); 5872 6134 5873 return(LL); 6135 5874 } 6136 5875 example 6137 { 6138 "EXAMPLE:"; echo = 2; 5876 { "EXAMPLE:"; echo = 2; 5877 5878 // EXAMPLE 1 5879 5880 // Conchoid, Pascal's Limacon. 5881 5882 // 1. Given a circle: x1^2+y1^24 5883 // 2. and a mover point M(x1,y1) on it 5884 // 3. Consider the fix point P(0,2) on the circle 5885 // 4. Consider the line l passing through M and P 5886 // 5. The tracer T(x,y) are the points on l at fixed distance 1 to M. 5887 6139 5888 if(defined(R)){kill R;} 6140 5889 ring R=(0,x,y),(x1,y1),dp; … … 6146 5895 locus(S96); 6147 5896 6148 // EXAMPLE 5897 // EXAMPLE 2 5898 6149 5899 // Consider two parallel planes z1=1 and z1=1, and two orthogonal parabolas on them. 6150 5900 // Determine the locus generated by the lines that rely the two parabolas … … 6152 5902 6153 5903 if(defined(R)){kill R;} 6154 ring R=(0,x,y,z),( lam,x2,y2,z2,x1,y1,z1), lp;5904 ring R=(0,x,y,z),(x2,y2,z2,z1,y1,x1,lam), lp; 6155 5905 short=0; 6156 5906 6157 5907 ideal L=z1+1, 6158 x1^2y1, 6159 z21, 6160 y2^2x2, 6161 4*x1*y21, 6162 xx1lam*(x2x1), 6163 yy1lam*(y2y1), 6164 zz1lam*(z2z1); 6165 6166 locus(L); // uses "moverdim",7 6167 6168 // Now observe the effect on the antiimage of option moverdim: 6169 // It does not take into account that lam is free. 6170 locus(L,"moverdim",3); 6171 6172 // Observe that with the option "mpverdim",3 the taxonomy becomes Special because considering only x1,y1,z1 as mover variables does not take into account that lam is free in the global antiimage. 5908 x1^2y1, 5909 z21, 5910 y2^2x2, 5911 4*x1*y21, 5912 xx1lam*(x2x1), 5913 yy1lam*(y2y1), 5914 zz1lam*(z2z1); 5915 5916 locus(L); // uses "moverdim",3 5917 // Observe the choose of the mover variables: the last 3 variables y1,x1,lam 5918 // If we choose x1,y1,z1 instead, the taxonomy becomes "Special" because 5919 // z1=1 is fix and do not really correspond to the mover variables. 5920 5921 // EXAMPLE 3 of parametric locus: 5922 5923 // Determining the equation of a general ellipse; 5924 // Uncentered elipse; 5925 5926 // Parameters (a,b,a0,b0,p): 5927 // a=large semiaxis, b=small semiaxis, 5928 // (a0,b0) = center of the ellipse, 5929 // (a1,b1) and (2*a0a1,2*b0b1) the focus, 5930 // p the slope of the line of the aaxis of the ellipse. 5931 5932 // Determine the equation of the ellipse. 5933 5934 // We must use the option "numpar",5 in order to consider 5935 // the first 5 parameters as free parameters for the locus 5936 5937 // Auxiliary variabes: 5938 // d1=distance from focus (a1,b1) to the mover point M(x1,y1), 5939 // d2=distance from focus (a2,b2) to the mover point M(x1,y1), 5940 // f=focus distance= distance from (a0,b0) to (a1,b1). 5941 5942 // Mover point (x1,y1) = tracer point (x,y). 5943 5944 if(defined(R1)){kill R1;} 5945 ring R1=(0,a,b,a0,b0,p,x,y),(d1,d2,f,a1,b1,x1,y1),lp; 5946 5947 ideal F3=b1b0p*(a1a0), 5948 //b2b0+p*(a1a0), 5949 //a1+a22*a0, 5950 //b1+b22*b0, 5951 f^2(a1a0)^2(b1b0)^2, 5952 f^2a^2b^2, 5953 (x1a1)^2+(y1b1)^2d1^2, 5954 (x12*a0+a1)^2+(y12*b0+b1)^2d2^2, 5955 d1+d22*a, 5956 xx1, 5957 yy1; 5958 5959 def G3=grobcov(F3); 5960 5961 def Loc3=locus(F3,"grobcov",G3,"numpar",5); Loc3; 5962 5963 // General ellipse: 5964 5965 def C=Loc3[1][1][1]; 5966 C; 5967 5968 // Centered ellipse of semiaxes (a,b): 5969 5970 def C0=subst(C,a0,0,b0,0,p,0); 5971 C0; 6173 5972 } 6174 5973 … … 6211 6010 } 6212 6011 example 6213 { 6214 "EXAMPLE:"; echo = 2; 6012 { "EXAMPLE:"; echo = 2; 6215 6013 if(defined(R)){kill R;}; 6216 6014 ring R=(0,a,b),(x,y),dp; 6217 6015 short=0; 6016 6218 6017 // Concoid 6219 6018 ideal S96=x^2+y^24,(b2)*xa*y+2*a,(ax)^2+(by)^21; … … 6300 6099 } 6301 6100 example 6302 { 6303 "EXAMPLE:"; echo = 2; 6101 { "EXAMPLE:"; echo = 2; 6304 6102 if(defined(R)){kill R;} 6305 6103 ring R=(0,x,y),(x1,y1),dp; 6306 6104 short=0; 6105 6307 6106 ideal S=x1^2+y1^24,(y2)*x1x*y1+2*x,(xx1)^2+(yy1)^21; 6308 6107 def L=locus(S); … … 6329 6128 the parameteres of the hypersurfaces, that are 6330 6129 considered as variables of the parametric ring. 6130 In the actual version, parametric envelope are accepted. 6131 To include fixed parameters a1,..ap, to the problem, one must 6132 declare them as the first parameters of the ring. if the 6133 the number of free parameters is p, the option \"numpar\",p 6134 is required. 6331 6135 Calling sequence: 6332 ring R=(0, x_1,..,x_n),(u_1,..,u_m),lp;6333 poly F=F( x_1,..,x_n,u_1,..,u_m);6334 ideal C=g_1( u_1,..u_m),..,g_s(u_1,..u_m);6136 ring R=(0,a1,..,ap,x_1,..,x_n),(u_1,..,u_m),lp; 6137 poly F=F(a1,..ap,x_1,..,x_n,u_1,..,u_m); 6138 ideal C=g_1(a1,..,ap,u_1,..u_m),..,g_s(a1,..ap,u_1,..u_m); 6335 6139 envelop(F,C[,options]); where s<m. 6336 x1,..,xn are the tr tacer variables.6140 x1,..,xn are the tracer variables. 6337 6141 u_1,..,u_m are the auxiliary variables. 6142 a1,..,ap are the fixed parameters if they exist 6143 If the problem is a parametric envelope, and a's exist, 6144 then the option \"numpar\",p m must be given. 6338 6145 By default the las n variables are the mover variables. 6146 See the EXAMPLE of parametric envelop by calling 6147 example envelop, 6339 6148 RETURN: The output is a list of the components [C_1, .. , C_n] 6340 6149 of the locus. Each component is given by … … 6365 6174 \"moverdim\", k: by default it is equal to n, the number of 6366 6175 xtracer variables. 6176 \"numpar\",p when fixed parameters are included 6367 6177 NOTE: grobcov and locus are called internally. 6368 The basering R, must be of the form Q[ x][u]6369 (x= parameters, u=variables).6178 The basering R, must be of the form Q[a,x][u] 6179 (x=variables, u=auxiliary variables), (a fixed parameters). 6370 6180 This routine uses the generalized definition of envelop 6371 6181 introduced in the book … … 6376 6186 { 6377 6187 def RR=basering; 6378 int tes=0; int i; int j; int k; int m; 6379 int d; 6380 int dp; 6381 ideal BBB; 6382 // new 6383 ideal CC=C; 6384 CC[size(CC)+1]=F; 6385 int movreal=size(variables(CC)); 6386 //Options 6188 list LRR=ringlist(RR); 6189 int nax=size(LRR[1][2]); 6190 int nuv=size(LRR[2]); 6191 6387 6192 list DD=#; 6388 6389 list Rx=ringlist(RR); 6390 ideal vmov; 6391 int moverdim=movreal; //size(Rx[1][2]); 6392 int nv=nvars(RR); 6393 if(moverdim>nv){moverdim=nv;} 6394 for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nvmoverdim);} 6395 //for(i=1;i<=nv;i++) {vmov[size(vmov)+1]=var(i);} 6396 int version=2; 6397 int comment=0; 6398 ideal Fm; 6399 int tax=1; 6400 6401 int familyinfo; 6402 for(i=1;i<=(size(DD) div 2);i++) 6403 { 6404 if(DD[2*i1]=="antiimage"){tax=DD[2*i];} 6405 if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 6406 if(DD[2*i1]=="version"){version=2;DD[2*i]=2;} 6407 if(DD[2*i1]=="comment"){comment=DD[2*i];} 6408 if(DD[2*i1]=="familyinfo"){familyinfo=0;DD[2*i]=0;} 6409 }; 6410 int ng=size(C); 6411 ideal S=F; 6412 for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];} 6413 int s=nvng; 6414 if(s>0) 6415 { 6416 matrix M[ng+1][ng+1]; 6417 def cc=comb(nv,ng+1); 6193 int na=0; 6194 int nx=nax; 6195 int nu=0; 6196 int nv=nuv; 6197 int i; int j; int k; 6198 //string("T_ nax=",nax,"; nx=",nx,"; nuv=",nuv,"; nv=",nv); 6199 int tnumpar=0; 6200 // int tnumvar=0; 6201 //"T_DD="; DD; 6202 for(i=1;i<=size(DD) div 2;i++) 6203 { 6204 if(DD[2*i1]=="numpar"){na=DD[2*i];tnumpar=1;} 6205 // if(DD[2*i1]=="numvar"){nv=DD[2*i];tnumvar=1;} 6206 } 6207 if(tnumpar==0){DD[size(DD)+1]="numpar"; DD[size(DD)+1]=na;} 6208 // if(tnumvar==0){DD[size(DD)+1]="numvar"; DD[size(DD)+1]=nv;} 6209 nx=naxna; 6210 nu=nuvnv; 6211 //string("T_ nax=",nax,"; nx=",nx,"; nuv=",nuv,"; nv=",nv); 6212 ideal Vnv; 6213 ideal Vnonv; 6214 for(i=1;i<=nu;i++){Vnonv[size(Vnonv)+1]=var(i);} 6215 //"T_Vnonv="; Vnonv; 6216 for(i=nu+1;i<=nuv;i++){Vnv[size(Vnv)+1]=var(i);} 6217 //"T_Vnv="; Vnv; 6218 ideal Cnor; 6219 ideal Cr=F; 6220 for(i=1;i<=size(C);i++) 6221 { 6222 if(subset(variables(C[i]),Vnonv)){Cnor[size(Cnor)+1]=C[i];} 6223 else{Cr[size(Cr)+1]=C[i];} 6224 } 6225 int nr=size(Cr); 6226 //string("T_nr=", nr,"; nv=",nv); 6227 if(nr>0) 6228 { 6229 matrix M[nr][nr]; 6230 def cc=comb(nv,nr); 6231 //"T_cc="; cc; 6232 //string("T_nv=",nv," nr=",nr); 6418 6233 poly J; 6419 6234 for(k=1;k<=size(cc);k++) 6420 6235 { 6421 for( j=1;j<=ng+1;j++)6236 for(i=1;i<=nr;i++) 6422 6237 { 6423 M[1,j]=diff(F,var(cc[k][j])); 6424 } 6425 for(i=1;i<=ng;i++) 6426 { 6427 for(j=1;j<=ng+1;j++) 6238 for(j=1;j<=nr;j++) 6428 6239 { 6429 M[i +1,j]=diff(C[i],var(cc[k][j]));6240 M[i,j]=diff(Cr[i],var(cc[k][j])); 6430 6241 } 6431 6242 } 6432 6243 J=det(M); 6433 S[size(S)+1]=J; 6434 } 6435 } 6436 if(comment>0){"System S before grobcov ="; S;} 6437 //def G=grobcov(S,DD); 6438 list HHH; 6439 DD[size(DD)+1]="moverdim"; // new option 6440 DD[size(DD)+1]=moverdim; 6441 // string("moverdim=",moverdim); 6244 Cr[size(Cr)+1]=J; 6245 } 6246 } 6247 ideal S=Cnor; 6248 for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];} 6249 for(i=1;i<=size(Cr);i++){S[size(S)+1]=Cr[i];} 6250 //"T_S="; S; 6442 6251 def L=locus(S,DD); 6443 // "L="; L;6444 list GL;6445 ideal fam; ideal env;6446 6447 def P=ring(Rx[1]);6448 Rx[1]=0;6449 def D=ring(Rx);6450 def RP=P+D;6451 list LL;6452 list NormalComp;6453 ideal Gi;6454 ideal BBBB;6455 poly B0;6456 6252 return(L); 6457 6253 } 6458 6254 example 6459 { 6460 "EXAMPLE:"; echo = 2; 6461 if(defined(R)){kill R;} 6462 ring R=(0,x,y),(r,s,y1,x1),lp; 6463 poly F=(xx1)^2+(yy1)^2r; 6464 ideal g=(x12*(s+r))^2+(y1s)^2s; 6465 6466 def E=envelop(F,g); 6467 E; 6468 6469 E[1][1]; 6470 6471 string("Taxonomy=", E[1][3][2]); 6472 6473 6474 // EXAMPLE 6255 { "EAXMPLE:"; echo=2; 6256 6257 // EXAMPLE 1 6475 6258 // Steiner Deltoid 6476 6259 // 1. Consider the circle x1^2+y1^21=0, and a mover point M(x1,y1) on it. 6477 6260 // 2. Consider the triangle A(0,1), B(1,0), C(1,0). 6478 6261 // 3. Consider lines passing through M perpendicular to two sides of ABC triangle. 6479 // 4. Obtain the envelop of the lines above. 6262 // 4. Determine the envelope of the lines above. 6263 6480 6264 if(defined(R)){kill R;} 6481 6265 ring R=(0,x,y),(x1,y1,x2,y2),lp; 6482 6266 short=0; 6267 6483 6268 ideal C=(x1)^2+(y1)^21, 6484 6269 x2+y21, … … 6486 6271 matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1; 6487 6272 poly F=det(M); 6488 // Curves Family F 6273 6274 // The lines of family F are 6489 6275 F; 6490 // Conditions C= 6276 6277 // The conditions C are 6491 6278 C; 6279 6492 6280 envelop(F,C); 6281 6282 // EXAMPLE 2 6283 // Parametric envelope 6284 6285 // Let c be the circle centered at the origin O(0,0) and having radius 1. 6286 // M(x1,y1) be a mover point gliding on c. 6287 // Let A(a0,b0) be a parametric fixed point: 6288 // Consider the set of lines parallel to the line AO passing thoug M. 6289 6290 // Determine the envelope of these lines 6291 6292 // We let the fixed point A coordinates as free parameters of the envelope. 6293 // We have to declare the existence of two parameters when 6294 // defining the ring in which we call envelop, 6295 // and set a0,b0 as the first variables of the parametric ring 6296 // The ring is thus 6297 6298 if(defined(R1)){kill R1;} 6299 ring R1=(0,a0,b0,x,y),(x1,y1),lp; 6300 short=0; 6301 6302 // The lines are F1 6303 poly F1=b0*(xx1)a0*(yy1); 6304 6305 // and the mover is on the circle c 6306 ideal C1=x1^2+y1^21; 6307 // The call is thus 6308 6309 def E1=envelop(F1,C1,"numpar",2); 6310 E1; 6311 6312 // The interesting first component EC1 is 6313 def EC1=E1[1][1][1]; 6314 EC1; 6315 6316 // that is equivalent to (a0*yb0*x)^2a0^2b0^2. 6317 // As expected it consists of the two lines 6318 // a0*yb0*x  sqrt(a0^2+b0^2), 6319 // a0*yb0*x + sqrt(a0^2+b0^2), 6320 // parallel to the line OM passing at the 6321 // points of the circle in the line perpendicular to OA. 6322 6323 // EXAMPLE 3 6324 // Parametric envelope 6325 6326 // Let c be the circle centered at the origin O(a1,b1) and having radiusr, 6327 // where a1,b1,r are fixed parameters 6328 // M(x1,y1) be a mover point gliding on c. 6329 // Let A(a0,b0) be a parametric fixed point: 6330 // Consider the set of lines parallel to the line AO passing thoug M. 6331 6332 // Determine the envelope of these lines 6333 6334 // We let the fixed point A,point M and r as free parameters of the envelope. 6335 // We have to declare the existence of 5 parameters when 6336 // defining the ring in which we call envelop, 6337 // and set a0,b0,a1,b1,r as the first variables of the parametric ring 6338 // The ring is thus 6339 6340 if(defined(R1)){kill R1;} 6341 ring R1=(0,a0,b0,a1,b1,r,x,y),(x1,y1),lp; 6342 short=0; 6343 6344 // The lines are F1 6345 poly F1=b0*(xx1)a0*(yy1); 6346 6347 // and the mover is on the circle c 6348 ideal C1=(x1a1)^2+(y1b1)^2r^2; 6349 // The call is thus 6350 6351 def E1=envelop(F1,C1,"numpar",5); 6352 E1; 6353 6354 // The interesting first component EC1 is 6355 def EC1=E1[1][1][1]; 6356 EC1; 6357 6358 // which corresponds to the product of two lines 6359 // parallel to the line AM and intercepting the circle 6360 // on the intersection of the line perpendicuar 6361 // to line AM passing through A 6493 6362 } 6494 6363 … … 6538 6407 ideal BBB; 6539 6408 //Options 6540 list DD=#;6409 // list DD=#; 6541 6410 ideal vmov; 6542 6411 int nv=nvars(RR); 6543 6412 for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);} 6544 int numpars=size(ringlist(RR)[1][2]);6545 int version=0;6546 if(nv<4){version=1;}6413 // int numpars=size(ringlist(RR)[1][2]); 6414 // int version=0; 6415 // if(nv<4){version=1;} 6547 6416 int comment=0; 6548 int familyinfo ;6417 int familyinfo=0; 6549 6418 ideal Fm; 6550 for(i=1;i<=(size(DD) div 2);i++)6551 {6552 if(DD[2*i1]=="vmov"){vmov=DD[2*i];}6553 // if(DD[2*i1]=="version"){version=DD[2*i];}6554 if(DD[2*i1]=="comment"){comment=DD[2*i];}6555 if(DD[2*i1]=="familyinfo"){familyinfo=DD[2*i];}6556 if(DD[2*i1]=="moreinfo"){moreinfo=DD[2*i];}6557 };6558 DD=list("vmov",vmov,"version",version,"comment",comment); 6419 // for(i=1;i<=(size(DD) div 2);i++) 6420 // { 6421 // if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 6422 // // if(DD[2*i1]=="version"){version=DD[2*i];} 6423 // if(DD[2*i1]=="comment"){comment=DD[2*i];} 6424 // if(DD[2*i1]=="familyinfo"){familyinfo=DD[2*i];} 6425 // if(DD[2*i1]=="moreinfo"){moreinfo=DD[2*i];} 6426 // }; 6427 // DD=list("vmov",vmov,"comment",comment); // ,"version",version 6559 6428 int ng=size(C); 6560 6429 ideal S=F; … … 6588 6457 } 6589 6458 //if(comment>0){"System S before grobcov ="; S;} 6590 def G=grobcov(S,DD); 6459 //"T_S="; S; 6460 def G=grobcov(S); // ,DD 6461 //"T_G=";G; 6591 6462 list GG; 6592 6463 for(i=2;i<=size(G);i++) … … 6595 6466 } 6596 6467 G=GG; 6468 //"T_G=";G; 6597 6469 if(moreinfo>0){return(G);} 6598 6470 else … … 6605 6477 //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1]; 6606 6478 //string("T_EE="); EE; 6607 if( G[i][3][1][1][1]==EE)6479 if(equalideals(G[i][3][1][1],EE)) 6608 6480 { 6609 6481 t=1; … … 6617 6489 } 6618 6490 example 6619 { 6620 "EXAMPLE:"; echo = 2; 6621 6491 { "EXAMPLE:"; echo = 2; 6622 6492 if(defined(R)){kill R;} 6623 6493 ring R=(0,x,y),(r,s,y1,x1),lp; … … 6647 6517 6648 6518 proc FamElemsAtEnvCompPoints(poly F,ideal C, ideal E,list #) 6649 "USAGE: FamElemsAtEnvCompPoints(poly F,ideal C, idealE);6519 "USAGE: FamElemsAtEnvCompPoints(poly F,ideal C,poly E); 6650 6520 poly F must be the family of hypersurfaces whose 6651 6521 envelope is analyzed. It must be defined in the ring … … 6691 6561 int i; 6692 6562 int moreinfo=0; 6563 int familyinfo=0; 6564 int comment=0; 6565 int numpar=0; 6566 ideal vmov; 6693 6567 list DD=#; 6694 6568 for(i=1;i<=(size(DD) div 2);i++) 6695 6569 { 6696 6570 if(DD[2*i1]=="vmov"){vmov=DD[2*i];} 6697 // if(DD[2*i1]=="version"){version=DD[2*i];}6698 6571 if(DD[2*i1]=="comment"){comment=DD[2*i];} 6699 6572 if(DD[2*i1]=="familyinfo"){familyinfo=DD[2*i];} 6700 6573 if(DD[2*i1]=="moreinfo"){moreinfo=DD[2*i];} 6574 if(DD[2*i1]=="numpar"){numpar=DD[2*i];} 6701 6575 }; 6702 6576 ideal S=C; … … 6723 6597 //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1]; 6724 6598 //string("T_EE="); EE; 6725 if(G[i][3][1][1][1]==E E)6599 if(G[i][3][1][1][1]==E) 6726 6600 { 6727 6601 t=1; … … 6734 6608 } 6735 6609 example 6736 { 6737 "EXAMPLE:"; echo = 2; 6738 6739 if(defined(R)){kill R;} 6740 ring R=(0,x,y),(t),dp; 6741 short=0; 6742 poly F=(x5*t)^2+y^29*t^2; 6743 ideal C; 6744 6745 def Env=envelop(F,C); 6746 Env; 6610 { "EXAMPLE:"; echo = 2; 6611 if(defined(R)){kill R;} 6612 ring R=(0,x,y),(t),dp; 6613 short=0; 6614 poly F=(x5*t)^2+y^29*t^2; 6615 ideal C; 6616 6617 def Env=envelop(F,C); 6618 Env; 6747 6619 6748 6620 // E is a component of the envelope: 6749 ideal E=Env[1][1]; 6750 6751 E; 6752 6753 def A=AssocTanToEnv(F,C,E); 6754 A; 6621 def E=Env[1][1][1]; 6622 E; 6623 6624 def A=AssocTanToEnv(F,C,E); 6625 A; 6755 6626 6756 6627 // The basis of the parameter values of the associated … … 6761 6632 // element at (x,y) is 6762 6633 6763 subst(F,t,(5/12)*y);6764 6765 def FE=FamElemsAtEnvCompPoints(F,C,E);6766 FE;6767 6768 factorize(FE[2][1]);6634 subst(F,t,(5/12)*y); 6635 6636 def FE=FamElemsAtEnvCompPoints(F,C,E); 6637 FE; 6638 6639 factorize(FE[2][1]); 6769 6640 6770 6641 // Thus the unique family element passing through the envelope point (x,y) … … 6772 6643 6773 6644 // EXAMPLE: 6774 if(defined(R)){kill R;}6775 ring R=(0,x,y),(r,s,y1,x1),lp;6776 6777 poly F=(xx1)^2+(yy1)^2r;6778 ideal g=(x12*(s+r))^2+(y1s)^2s;6779 6780 def E=envelop(F,g);6781 E;6782 6783 def A=AssocTanToEnv(F,g,E[1][1][1]);6784 A;6785 6786 def M1=coef(A[2][1],x1);6787 def M2=coef(A[2][2],y1);6788 def M3=coef(A[2][3],s);6789 def M4=coef(A[2][4],r);6645 if(defined(R)){kill R;} 6646 ring R=(0,x,y),(r,s,y1,x1),lp; 6647 6648 poly F=(xx1)^2+(yy1)^2r; 6649 ideal g=(x12*(s+r))^2+(y1s)^2s; 6650 6651 def E=envelop(F,g); 6652 E; 6653 6654 def A=AssocTanToEnv(F,g,E[1][1][1]); 6655 A; 6656 6657 def M1=coef(A[2][1],x1); 6658 def M2=coef(A[2][2],y1); 6659 def M3=coef(A[2][3],s); 6660 def M4=coef(A[2][4],r); 6790 6661 6791 6662 // The parameter values corresponding to the family 6792 6663 // element tangent at point (x,y) of the envelope are: 6793 "x1=";M1[2,2]/M1[2,1];6794 6795 "y1=";M2[2,2]/M2[2,1];6796 6797 "s=";M3[2,2]/M3[2,1];6798 6799 "r=";M4[2,2]/M4[2,1];6664 "x1=";M1[2,2]/M1[2,1]; 6665 6666 "y1=";M2[2,2]/M2[2,1]; 6667 6668 "s=";M3[2,2]/M3[2,1]; 6669 6670 "r=";M4[2,2]/M4[2,1]; 6800 6671 6801 6672 // Now detect if there are other family elements passing at this point: 6802 def FE=FamElemsAtEnvCompPoints(F,g,E[1][1][1]);6803 FE;6673 def FE=FamElemsAtEnvCompPoints(F,g,E[1][1][1]); 6674 FE; 6804 6675 6805 6676 // FE[1] is the set of lpp. It has dimension 42=2. … … 6852 6723 } 6853 6724 example 6854 { 6855 "EXAMPLE:"; echo = 2; 6725 { "EXAMPLE:"; echo = 2; 6856 6726 if(defined(R)){kill R;} 6857 6727 ring R=(0,a,b,c),(x,y),dp; … … 6862 6732 } 6863 6733 6864 // AddLocus: auxil iary routine for locus0 that computes the components of the constructible:6734 // AddLocus: auxilliary routine for locus0 that computes the components of the constructible: 6865 6735 // Input: the list of locally closed sets to be added, each with its type as third argument 6866 6736 // L=[ [LC[11],..,LC[1k_1],.., [LC[r1],..,LC[rk_r] ] where … … 7128 6998 } 7129 6999 // B=reduced basis on CR 7000 //"T_PR="; PR; 7001 //"T_CR="; CR; 7002 //"T_B="; B; 7130 7003 if(rep==1){return(list(lcc,B,CR));} 7131 7004 else{return(list(lcc,B,PR));} … … 7142 7015 } 7143 7016 example 7144 { 7145 "EXAMPLE:"; echo = 2; 7146 if(defined(RE)){kill RE;} 7147 ring RE=(0,a,b,c,d,e,f),(x,y),lp; 7148 ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2; 7149 ideal A=a*eb*d; 7017 { "EXAMPLE:"; echo = 2; 7018 if(defined(RE)){kill RE;} 7019 ring RE=(0,a,b,c,d,e,f),(x,y),lp; 7020 ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2; 7021 ideal A=a*eb*d; 7150 7022 7151 7023 WLemma(F,A); … … 7228 7100 { 7229 7101 G[size(G)+1]=G0; 7102 //"T_G0="; G0; 7230 7103 N=G0[3][1][2]; 7104 //"T_N="; N; 7231 7105 for(i=1;i<=size(N);i++) 7232 7106 { … … 7236 7110 Epend[size(Epend)+1]=N[i]; 7237 7111 } 7112 //"T_i="; i; 7238 7113 } 7114 //"T_Etot="; Etot; 7115 //"T_Epend=";Epend; 7239 7116 } 7240 7117 } … … 7242 7119 } 7243 7120 example 7244 { 7245 "EXAMPLE:"; echo = 2; 7246 if(defined(RRR)){kill RRR;} 7247 ring RRR=(0,b,c,d,e,f),(x,y,t),lp; 7248 short=0; 7249 ideal S=x^2+2*c*x*y+2*d*x*t+b*y^2+2*e*y*t+f*t^2, 7121 { "EXAMPLE:"; echo = 2; 7122 if(defined(RRR)){kill RRR;} 7123 ring RRR=(0,b,c,d,e,f),(x,y,t),lp; 7124 short=0; 7125 ideal S=x^2+2*c*x*y+2*d*x*t+b*y^2+2*e*y*t+f*t^2, 7250 7126 x+c*y+d*t,c*x+b*y+e*t; 7251 7252 7127 grobcov(S); 7128 WLcgs(S); 7253 7129 } 7254 7130 … … 7312 7188 int i; int j; 7313 7189 ideal J=idminusid(E,N); 7190 //"T_J="; J; 7314 7191 for(i=1;i<=size(F);i++) 7315 7192 { … … 7422 7299 } 7423 7300 example 7424 { 7425 "EXAMPLE:"; echo = 2; 7426 if (defined(R)) {kill R;} 7427 ring R=(0,x,y),(x1,y1,x2,y2),lp; 7428 ideal F=y*x1+(x1)*y1+y, 7429 (x1)*(x1+1)+y*y1, 7430 y*x2+(x+1)*y2y, 7431 (x+1)*(x21)+y*y2, 7432 (x1x)^2+y1^2(x1x)^2y2^2; 7301 { "EAXMPLE:"; echo = 2; 7302 if (defined(R)) {kill R;} 7303 ring R=(0,x,y),(x1,y1,x2,y2),lp; 7304 ideal F=y*x1+(x1)*y1+y, 7305 (x1)*(x1+1)+y*y1, 7306 y*x2+(x+1)*y2y, 7307 (x+1)*(x21)+y*y2, 7308 (x1x)^2+y1^2(x1x)^2y2^2; 7433 7309 7434 7310 def G=grobcov(F,"rep",1); 7435 7436 7311 G; 7437 7312 7438 7313 def L=Grob1Levels(G); 7439 7440 7314 L; 7441 7315 7442 7316 def LL=Levels(L); 7443 7444 7317 LL; 7445 7318 } … … 7463 7336 setring(RP); 7464 7337 def SP=imap(RR,S); 7338 //"T_SP="; SP; 7339 //"T_typeof(SP[1])="; typeof(SP[1]); 7465 7340 ideal EP; 7466 7341 EP=SP[1]; … … 7476 7351 } 7477 7352 example 7478 { 7479 "EXAMPLE:"; echo = 2; 7353 { "EAXMPLE:"; echo = 2; 7480 7354 if(defined(R)){kill R;} 7481 7355 ring R=(0,x,y,z),(x1,y1,z1),lp; 7356 7482 7357 ideal I1=x+y*z*x1; 7483 7358 ideal I2=xy*z*y1; … … 7488 7363 intersectpar(S); 7489 7364 } 7490 7491 7365 7492 7366 proc ADGT(ideal H,ideal T,ideal H1,ideal T1,list #) … … 7614 7488 int r=m div 2; 7615 7489 if(m mod 2==0){r=r1;} 7616 for(i=1;i<=r;i++) 7490 //"L0="; L0; 7491 for(i=1;i<=r;i++) 7617 7492 { 7618 7493 G1=grobcov(F,"null",L0[2*i],"nonnull",L0[2*i+1],"rep",1); … … 7677 7552 } 7678 7553 example 7679 { 7680 "EXAMPLE:"; echo = 2; 7681 7554 { "EXAMPLE:"; echo = 2; 7682 7555 // Determine the supplementary conditions 7683 7556 // for the nondegenerate triangle A(1,0), B(1,0), C(x,y) … … 7691 7564 7692 7565 ideal H=y*x1+(x1)*y1+y, 7693 7694 7695 7566 (x1)*(x1+1)+y*y1, 7567 y*x2+(x+1)*y2y, 7568 (x+1)*(x21)+y*y2; 7696 7569 7697 7570 // Thesis T: the orthic triangle is isosceles 7698 ideal T=(x1x)^2+y1^2(x2x)^2y2^2;7571 ideal T=(x1x)^2+y1^2(x2x)^2y2^2; 7699 7572 7700 7573 // Negative Hypothesis H1: ABC is nondegenerate 7701 7574 ideal H1=y; 7702 7575 7703 7576 // Negative Thesis T1: the orthic triangle is nondegenerate 7704 7577 ideal T1=x*(y1y2)y*(x1x2)+x1*y2x2*y1; 7705 7578 … … 7710 7583 ADGT(H,T,H1,T1); 7711 7584 7712 // Now using dif ference of constructible sets for negative hypothesis and thesis7585 // Now using diference of constructible sets for negative hypothesis and thesis 7713 7586 7714 7587 ADGT(H,T,H1,T1,"neg",1); 7715 7588 7716 7589 // The results are identical using both methods for the negative propositions 7717 // 7718 // 7719 7720 // EXAMPLE 7590 //  Rabinovitch or 7591 //  DifConsLCSets 7592 7593 // EXAMPLE 2 7721 7594 7722 7595 // Automatic Theorem Proving. … … 7741 7614 7742 7615 if (defined(R1)){kill R1;} 7743 7744 7616 ring R1=(0,a,b),(x1,y1,x2,y2,x0,y0,x3,y3,r2),dp; 7745 7746 7617 short=0; 7747 7618 … … 7763 7634 (x0+a2*x3)^2+(y0+b2*y3)^2r2; 7764 7635 7765 7636 ADGT(H,T,b,1); 7766 7637 7767 7638 // Thus the nine point circle theorem is true for all real points excepting b=0. … … 7797 7668 list G2=grobcov(F,"rep",1); 7798 7669 def L2=Grob1Levels(G2); 7670 //"L2="; L2; 7799 7671 ideal FN=T1; 7800 7672 if(not(equalideals(H1,ideal(1)))) … … 7811 7683 list G3=grobcov(F,"rep",1); 7812 7684 def L3=Grob1Levels(G3); 7685 //"L3="; L3; 7813 7686 } 7814 7687 def LL=DifConsLCSets(L2,L3); 7688 //"LL="; LL; 7815 7689 LL=ConsLevels(LL); 7816 7690 LL=Levels(LL); 7691 //"LL="; LL; 7817 7692 if(rep==1){return(LL);} 7818 7693 else
Note: See TracChangeset
for help on using the changeset viewer.