Changeset 2e7410 in git
 Timestamp:
 Sep 24, 2010, 4:52:16 PM (13 years ago)
 Branches:
 (u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')
 Children:
 f0e080284bc1a4ba286ba66df791c785e808feb8
 Parents:
 1492bbce4e131a22de468dd41b9ced71d0fe3f93
 Location:
 Singular/LIB
 Files:

 1 added
 1 edited
 1 moved
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/assprimeszerodim.lib

Property
mode
changed from
100644
to100755
r1492bb r2e7410 1 1 //////////////////////////////////////////////////////////////////////////////// 2 2 version="$Id$"; 3 category = "Commutative algebra";3 category = "Commutative Algebra"; 4 4 info=" 5 LIBRARY: ass Primes.libassociated primes of a zerodimensional ideal5 LIBRARY: assprimeszerodim.lib associated primes of a zerodimensional ideal 6 6 7 7 AUTHORS: N. Idrees nazeranjawwad@gmail.com … … 11 11 OVERVIEW: 12 12 13 A library for computing the associated primes 13 A library for computing the associated primes and the radical of a 14 14 zerodimensional ideal in the polynomial ring over the rational numbers, 15 Q[x_1,...,x_n] using modular computations.15 Q[x_1,...,x_n], using modular computations. 16 16 17 17 SEE ALSO: primdec_lib 18 18 19 19 PROCEDURES: 20 zeroR(I); computes the radical of I 20 21 assPrimes(I); computes the associated primes of I 21 zeroR(I); compute the radical of I22 22 "; 23 23 … … 30 30 "USAGE: zeroR(I,[n]); I ideal, optional: n number of processors (for parallel computing) 31 31 ASSUME: I is zerodimensional in Q[variables] 32 NOTE: Parallelization is just applicable using 32bit Singular version since 33 MPlinks are not compatible with 64bit Singular version. 32 34 RETURN: the radical of I 33 35 EXAMPLE: example zeroR; shows an example … … 64 66 65 67 // Initialize the list of primes  66 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10) 67 intvec L = primeList(10); 68 intvec L = primeList(I,10); 68 69 L[5] = prime(random(100000000,536870912)); 69 70 … … 141 142 } 142 143 143 // hier deleteUnlucky einbauen, streichen, wenn der Grad nicht stimmt144 // insert deleteUnluckyPrimes 144 145 G = chinrem(CO1,CO2); 145 146 N = CO2[1]; … … 160 161 161 162 j = size(L) + 1; 162 163 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10,L) 164 L = primeList(10,L); 163 L = primeList(I,10,L); 165 164 166 165 if(n > 1) … … 198 197 199 198 proc assPrimes(list #) 200 "USAGE: assPrimes(I,[ n],[a]); I ideal or module,199 "USAGE: assPrimes(I,[a],[n]); I ideal or module, 201 200 optional: n number of processors (for parallel computing), a 202  a =1: method of Eisenbud/Hunecke/Vasconcelos203  a =2: method of Gianni/Trager/Zacharias204  a =3: mathod of Monico205 assPrimes(I) chooses n=a=1201  a = 1: method of Eisenbud/Hunecke/Vasconcelos 202  a = 2: method of Gianni/Trager/Zacharias 203  a = 3: mathod of Monico 204 assPrimes(I) chooses n = a = 1 206 205 ASSUME: I is zerodimensional over Q[variables] 206 NOTE: Parallelization is just applicable using 32bit Singular version since 207 MPlinks are not compatible with 64bit Singular version. 207 208 RETURN: a list pr of associated primes of I: 208 209 EXAMPLE: example assPrimes; shows an example … … 214 215 215 216 ideal I; 216 if(typeof(#[1]) =="ideal")217 if(typeof(#[1]) == "ideal") 217 218 { 218 219 I = #[1]; … … 223 224 I = Ann(M); 224 225 } 226 227 // Initialize optional parameter  228 if(size(#) > 1) 229 { 230 if(size(#) == 2) 231 { 232 int alg = #[2]; 233 int n = 1; 234 } 235 if(size(#) == 3) 236 { 237 int alg = #[2]; 238 int n = #[3]; 239 if((n > 1) && (1  system("with","MP"))) 240 { 241 "========================================================================"; 242 "There is no MP available on your system. Since this is necessary to "; 243 "parallelize the algorithm, the computation will be done without forking."; 244 "========================================================================"; 245 n = 1; 246 } 247 } 248 } 249 else 250 { 251 int alg = 1; 252 int n = 1; 253 } 225 254 226 255 def SPR = basering; 227 I = std(I);256 I = modStd(I,n); 228 257 int d = vdim(I); 229 258 if(d == 1) { ERROR("Ideal is not zerodimensional."); } … … 263 292 "CPUtime for radicalcheck is "+string(timer  T)+" seconds."; 264 293 } 265 // Initialize optional parameter 266 if(size(#) > 1)267 {268 if(size(#) == 2)269 {270 int alg = #[2];271 int n = 1;272 }273 if(size(#) == 3)274 {275 int alg = #[2];276 int n = #[3];277 if((n > 1) && (1  system("with","MP")))278 {279 "========================================================================";280 "There is no MP available on your system. Since this is necessary to ";281 "parallelize the algorithm, the computation will be done without forking.";282 "========================================================================";283 n = 1;284 }285 }286 }287 else288 {289 int alg = 1;290 int n = 1;291 }292 294 293 295 export(SPR); … … 297 299 export(I_for_fork); // I available for each link 298 300 301 // Initialize the list of primes  302 intvec L = primeList(I,10); 303 L[5] = prime(random(1000000000,2134567879)); 304 299 305 list Re; 300 306 … … 308 314 int j = 1; 309 315 int index = 1; 310 311 312 // Initialize the list of primes 313 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10)314 intvec L = primeList(10);315 L[5] = prime(random(1000000000,2134567879));316 316 317 317 // If there is more than one processor available, we parallelize the  … … 359 359 for(i = 1; i <= n; i++) 360 360 { 361 if(status(l(i), "read", "ready")) 361 if(status(l(i), "read", "ready")) // ask if link l(i) is ready otherwise sleep for t seconds 362 362 { 363 P = read(l(i)); 363 P = read(l(i)); // read the result from l(i) 364 364 CO1[index] = P[1]; 365 365 CO2[index] = bigint(P[2]); … … 378 378 } 379 379 } 380 if(k == n) 380 if(k == n) // k describes the number of closed links 381 381 { 382 382 j++; … … 401 401 "CPUtime for computing list in assPrimes is "+string(timer  tt)+" seconds."; 402 402 } 403 //hier deleteUnlucky einbauen, streichen, wenn der Grad nicht stimmt 403 404 // insert deleteUnluckyPrimes 404 405 G = chinrem(CO1,CO2); 405 406 N = CO2[1]; … … 461 462 { 462 463 H[1][i] = quickSubst(H[1][i],f,I); 463 464 if(typeof(#[1])=="ideal") 465 { 466 Re[i1] = #[1] + ideal(H[1][i]); 467 } 468 else 469 { 470 Re[i1] = M + H[1][i]*freemodule(nrows(M)); 471 } 464 Re[i1] = I + ideal(H[1][i]); 472 465 } 473 466 … … 514 507 j = size(L) + 1; 515 508 516 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10,L) 517 L = primeList(10,L); 509 setring(SPR); 510 L = primeList(I,10,L); 511 setring rHelp; 518 512 519 513 if(n > 1) … … 746 740 //////////////////////////////////////////////////////////////////////////////// 747 741 748 static proc primeList(int n, list #)749 "USAGE: primeList(n); (resp. primeList(n,L); )750 RETURN: the intvec of n greatest primes <= 536870912 (resp. n greatest751 primes < L[size(L)] union with L)752 EXAMPLE:example primList; shows an example753 "754 {755 intvec L;756 int i,p;757 if(size(#) == 0)758 {759 p = prime(536870912);760 L[1] = p;761 }762 else763 {764 L = #[1];765 p = prime(L[size(L)]1);766 L[size(L)+1] = p;767 }768 if(p == 2){ERROR("no more primes");}769 for(i = 2; i <= n; i++)770 {771 p = prime(p1);772 L[size(L)+1] = p;773 }774 return(L);775 }776 example777 { "EXAMPLE:"; echo = 2;778 intvec L=primeList(10);779 size(L);780 L[size(L)];781 L=primeList(5,L);782 size(L);783 L[size(L)];784 }785 786 ////////////////////////////////////////////////////////////////////////////////787 788 742 static proc zeroRadP(ideal I, int p) 789 743 { … … 874 828 875 829 //Test for zeroR 876 ring R = 0, (x,y),dp;830 ring R = 0,(x,y),dp; 877 831 ideal I = xy42xy2+x, x2x, y42y2+1; 878 832 879 833 //Cyclic 6 880 ring R =0,x(1..6),dp;881 ideal I =cyclic(6);834 ring R = 0,x(1..6),dp; 835 ideal I = cyclic(6); 882 836 883 837 //Amrhei 884 ring R=0,(a,b,c,d,e,f),dp; 885 ideal I= 886 2fb+2ec+d2+a2+a, 887 2fc+2ed+2ba+b, 888 2fd+e2+2ca+c+b2, 889 2fe+2da+d+2cb, 890 f2+2ea+e+2db+c2, 891 2fa+f+2eb+2dc; 838 ring R = 0,(a,b,c,d,e,f),dp; 839 ideal I = 2fb+2ec+d2+a2+a, 840 2fc+2ed+2ba+b, 841 2fd+e2+2ca+c+b2, 842 2fe+2da+d+2cb, 843 f2+2ea+e+2db+c2, 844 2fa+f+2eb+2dc; 892 845 893 846 894 847 //BeckerNiermann 895 848 ring R = 0,(x,y,z),dp; 896 ideal I= 897 x2+xy2z2xy+y4+y2+z2, 898 x3y2+xy2z+xyz32xy+y4, 899 2x2y+xy4+yz43; 849 ideal I = x2+xy2z2xy+y4+y2+z2, 850 x3y2+xy2z+xyz32xy+y4, 851 2x2y+xy4+yz43; 900 852 901 853 //Roczen 902 ring R=0,(a,b,c,d,e,f,g,h,k,o),dp; 903 ideal I= 904 o+1, 905 k4+k, 906 hk, 907 h4+h, 908 gk, 909 gh, 910 g3+h3+k3+1, 911 fk, 912 f4+f, 913 eh, 914 ef, 915 f3h3+e3k3+e3+f3+h3+k3+1, 916 e3g+f3g+g, 917 e4+e, 918 dh3+dk3+d, 919 dg, 920 df, 921 de, 922 d3+e3+f3+1, 923 e2g2+d2h2+c, 924 f2g2+d2k2+b, 925 f2h2+e2k2+a; 854 ring R = 0,(a,b,c,d,e,f,g,h,k,o),dp; 855 ideal I = o+1, k4+k, hk, h4+h, gk, gh, g3+h3+k3+1, 856 fk, f4+f, eh, ef, f3h3+e3k3+e3+f3+h3+k3+1, 857 e3g+f3g+g, e4+e, dh3+dk3+d, dg, df, de, 858 d3+e3+f3+1, e2g2+d2h2+c, f2g2+d2k2+b, 859 f2h2+e2k2+a; 926 860 927 861 //4 bodies with equal masses, before symmetrisation. 928 862 //We are looking for the real positive solutions 929 ring R=0,(B,D,F,b,d,f),dp; 930 931 ring R=0,(B,b,D,d,F,f),dp; 932 ideal I=(bd)*(BD)2*F+2, 933 (bd)*(B+D2*F)+2*(BD), 934 (bd)^22*(b+d)+f+1, 935 B^2*b^31,D^2*d^31,F^2*f^31; 936 937 ring R=0,(B,b,D,d,F,f),dp; 938 ideal I=(bd)*(BD)F+1, 939 (bd)*(B+DF)+(BD), 940 (bd)^2(b+d)+f+1, 941 B^2*b^31,D^2*d^31,F^2*f^31; 863 ring R = 0,(B,b,D,d,F,f),dp; 864 ideal I = (bd)*(BD)2*F+2, 865 (bd)*(B+D2*F)+2*(BD), 866 (bd)^22*(b+d)+f+1, 867 B^2*b^31,D^2*d^31,F^2*f^31; 868 869 ring R = 0,(B,b,D,d,F,f),dp; 870 ideal I = (bd)*(BD)F+1, 871 (bd)*(B+DF)+(BD), 872 (bd)^2(b+d)+f+1, 873 B^2*b^31,D^2*d^31,F^2*f^31; 942 874 943 875 944 876 //prime 945 ring R=0,(x,y,z,t,u),dp; 946 947 ideal I=2x22y2+2z22t2+2u21, 877 ring R = 0,(x,y,z,t,u),dp; 878 ideal I = 2x22y2+2z22t2+2u21, 948 879 2x32y3+2z32t3+2u31, 949 2x42y4+2z42t4+2u41,950 2x52y5+2z52t5+2u51,951 2x62y6+2z62t6+2u61;880 2x42y4+2z42t4+2u41, 881 2x52y5+2z52t5+2u51, 882 2x62y6+2z62t6+2u61; 952 883 953 884 */ 
Property
mode
changed from

Singular/LIB/surfacesignature.lib
r1492bb r2e7410 2 2 category="Singularities"; 3 3 info=" 4 LIBRARY: surfacesignature.lib signature of surface singularity 5 AUTHORS: Gerhard Pfister, pfister@mathematik.unikl.de 6 Muhammad Ahsan Banyamin, ahsanbanyamin@gmail.com 7 8 OVERVIEW: 9 A library for computing the signature of irreducible surface singularity. 10 It contains also a procedure for computing the newton pairs of a surface 11 singularity. 12 13 THEORY: The signature of a surface singularity is defined in 14 [Durfee,A.:The Signature of Smoothings of Complex Surface 15 Singularities, Math. Ann.,232,8598 (1978)]. 16 The algorithm we use has been proposed in 17 [Nemethi,A.:The signature of f(x,y)+z^n, Proceedings of Singularity 18 Conference (C.T.C Wall's 60th birthday meeting),Liverpool 1996, 19 London Math.Soc. LN 263(1999),131149]. 4 LIBRARY: signaturesurfsing.lib signature of surface singularity 5 6 AUTHORS: Gerhard Pfister pfister@mathematik.unikl.de 7 @* Muhammad Ahsan Banyamin ahsanbanyamin@gmail.com 8 @* Stefan Steidel steidel@mathematik.unikl.de 9 10 OVERVIEW: 11 12 A library for computing the signature of irreducible surface singularity. 13 The signature of a surface singularity is defined in [Durfee, A.: The 14 Signature of Smoothings of Complex Surface Singularities, Math. Ann., 232, 15 8598 (1978)]. The algorithm we use has been proposed in [Nemethi, A.: The 16 signature of f(x,y)+z^N, Proceedings of Singularity Conference (C.T.C. Wall's 17 60th birthday meeting), Liverpool 1996, London Math.Soc. LN 263(1999), 18 131149]. 20 19 21 20 PROCEDURES: 22 brieskornSign(a,b,c); signature of Brieskorn singularity x^a+y^b+z^c=0 23 newtonpairs(f); newton pairs of surface singularity 24 signature(N,f); signature of singularity z^+f(x,y)=0,f irreducible 21 brieskornSign(a1,a2,a3); signature of Brieskorn singularity x^a1+y^a2+z^a3 22 signature(N,f); signature of singularity z^N+f(x,y)=0, f irreducible 25 23 "; 26 24 27 25 LIB "hnoether.lib"; 28 /////////////////////////////////////////////////////////////////////////////// 29 proc signature(int N,poly f) 26 LIB "alexpoly.lib"; 27 LIB "gmssing.lib"; 28 29 /////////////////////////////////////////////////////////////////////////////// 30 // sigma(z^N + f) in terms of Puiseux pairs of f for f irreducible  31 32 static proc exponentSequence(poly f) 33 //=== computes the sequence a_1,...,a_s of exponents as described in [Nemethi] 34 //=== using the Puiseux pairs (m_1, n_1),...,(m_s, n_s) of f: 35 //===  a_1 = m_1, 36 //===  a_i = m_i  n_i * (m_[i1]  n_[i1] * a_[i1]). 37 //=== 38 //=== Return: list of two intvecs: 39 //=== 1st entry: A = (a_1,...,a_s) 40 //=== 2nd entry: N = (n_1,...,n_s) 41 { 42 def R = basering; 43 ring S = 0,(x,y),dp; 44 poly f = fetch(R,f); 45 list puiseuxPairs = invariants(f); 46 setring R; 47 48 intvec M = puiseuxPairs[1][3]; 49 intvec N = puiseuxPairs[1][4]; 50 51 int i; 52 int a = M[1]; 53 intvec A = a; 54 for(i = 2; i <= size(M); i++) 55 { 56 a = M[i]  N[i] * (M[i1]  N[i1] * a); 57 A[size(A)+1] = a; 58 } 59 60 return(list(A,N)); 61 } 62 example 63 { "EXAMPLE:"; echo = 2; 64 ring r = 0,(x,y),dp; 65 exponentSequence(y4+2x3y2+x6+x5y); 66 } 67 68 /////////////////////////////////////////////////////////////////////////////// 69 70 proc brieskornSign(a1,a2,a3) 71 "USAGE: brieskornSign(a1,a2,a3); a1,a2,a3 = integers 72 RETURN: signature of Brieskorn singularity x^a1+y^a2+z^a3 73 EXAMPLE: example brieskornSign; shows an example 30 74 " 31 USAGE: signature(N,f); N=integer, f=irreducible poly in 2 variables 32 RETURN: signature of surface singularity defined by z^N+f=0 33 EXAMPLE: example signature; shows an example 75 { 76 int a_temp, t, k1, k2, k3, s_t, sigma; 77 number s; 78 79 if(a1 > a2) { a_temp = a1; a1 = a2; a2 = a_temp; } 80 if(a2 > a3) { a_temp = a2; a2 = a3; a3 = a_temp; } 81 if(a1 > a2) { a_temp = a1; a1 = a2; a2 = a_temp; } 82 83 for(t = 0; t <= 2; t++) 84 { 85 s_t = 0; 86 for(k1 = 1; k1 <= a11; k1++) 87 { 88 for(k2 = 1; k2 <= a21; k2++) 89 { 90 for(k3 = 1; k3 <= a31; k3++) 91 { 92 s = number(k1)/a1 + number(k2)/a2 + number(k3)/a3; 93 if(t < s) 94 { 95 if(s < t+1) 96 { 97 s_t = s_t + 1; 98 } 99 else 100 { 101 break; 102 } 103 } 104 } 105 if(k3 == 1) { break; } 106 } 107 if(k2 == 1) { break; } 108 } 109 sigma = sigma + (1)^t * s_t; 110 } 111 return(sigma); 112 } 113 example 114 { "EXAMPLE:"; echo = 2; 115 ring R = 0,x,dp; 116 brieskornSign(11,3,5); 117 } 118 119 /////////////////////////////////////////////////////////////////////////////// 120 121 static proc signatureP(int N,poly f) 122 "USAGE: signatureP(N,f); N = integer, f = irreducible poly in 2 variables 123 RETURN: signature of surface singularity defined by z^N + f(x,y) = 0 124 EXAMPLE: example signatureP; shows an example 34 125 " 35 126 { 36 def R=basering; 37 ring S=0,(x,y),dp; 38 poly f = fetch(R,f); 39 list L=newtonpairs(f); 40 setring R; 41 if(size(L)>2) 42 { 43 return(Generalcase(N,L)); 44 } 45 if(size(L)==2) 46 { 47 return(signtwopairs(N,L)); 48 } 49 return(brieskornSign(L[1][2],L[1][1],N)); 50 } 51 example 127 int i, d, prod, sigma; 128 list L = exponentSequence(f); 129 int s = size(L[2]); 130 131 if(s == 1) 132 { 133 return(brieskornSign(L[1][1], L[2][1], N)); 134 } 135 136 prod = 1; 137 sigma = brieskornSign(L[1][s], L[2][s], N); 138 for(i = s  1; i >= 1; i) 139 { 140 prod = prod * L[2][i+1]; 141 d = gcd(N, prod); 142 sigma = sigma + d * brieskornSign(L[1][i], L[2][i], N/d); 143 } 144 145 return(sigma); 146 } 147 example 52 148 { "EXAMPLE:"; echo = 2; 53 149 ring r = 0,(x,y),dp; 54 150 int N = 3; 55 poly f= x1521x14+8x13y6x1316x12y+20x11y2x12+8x11y36x10y2 56 +24x9y3+4x9y216x8y3+26x7y46x6y4+8x5y5+4x3y6y8; 57 signature(N,f); 58 } 59 60 /////////////////////////////////////////////////////////////////////////// 61 proc newtonpairs(poly f) 62 "USAGE: newtonpairs(f); f= irreducible bivariate poly 63 RETURN: newton pairs of curve singularity f(x,y)=0 64 EXAMPLE: example newtonpairs; shows an example 65 " 66 { 67 def R=basering; 68 ring S=0,(x,y),dp; 69 poly f = fetch(R,f); 70 list M=invariants(f); 71 setring R; 72 list L=M[1][3]; 73 list K=M[1][4]; 74 int i,j; 75 intvec v; 76 list N1,N2,P,T; 77 N1[1]=M[1][4][1]; 78 N2[1]=M[1][3][1]; 79 for(i=2;i<=size(M[1][3]);i++) 80 { 81 N1[i]=M[1][4][i]; 82 N2[i]=M[1][3][i]N1[i]*M[1][3][i1]; 83 } 84 P[1]=N1; 85 P[2]=N2; 86 for(j=1;j<=size(P[1]);j++) 87 { 88 v=P[1][j],P[2][j]; 89 T[j]=v; 90 } 91 return(T); 92 } 93 example 94 { "EXAMPLE:"; echo = 2; 95 ring r=0,(x,y),dp; 96 newtonpairs(y4+2x3y2+x6+x5y); 97 } 98 99 /////////////////////////////////////////////////////////////////////////// 100 proc brieskornSign(a,b,c) 101 "USAGE: brieskornSign(a,b,c);a,b,c=integers 102 RETURN: signature of Brieskorn singularity x^a+y^b+z^c 103 EXAMPLE: example brieskornSign; shows an example 104 " 105 { 106 int i,j,k,d; 107 for(i=1;i<=a1;i++) 108 { 109 for(j=1;j<=b1;j++) 110 { 111 for(k=1;k<=c1;k++) 112 { 113 d=d+signa(number(i)/a+number(j)/b+number(k)/c); 114 } 115 } 116 } 117 return(d); 118 } 119 example 120 { "EXAMPLE:"; echo = 2; 121 ring R=0,x,dp; 122 brieskornSign(11,3,5); 123 } 124 /////////////////////////////////////////////////////////////////////////// 125 static proc signsin(number n) 126 "USAGE: signsin(n); n=number 127 RETURN: the sign of sin(n) (sin(n) = the value of the sine of n) 128 EXAMPLE: example signSin; shows an example 129 " 130 { 131 def r=basering; 132 int a; 133 int b=10; 134 ring s=(real,10),x,dp; 135 number m=imap(r,n); 136 number pi = number_pi(11); 137 number t=m*pi; 138 poly f=sin(b); 139 poly h=subst(f,x,t); 140 if(h>0){a=1;} 141 if(h<0){a=1;} 142 setring r; 143 return(a); 144 } 145 example 146 { "EXAMPLE:"; echo = 2; 147 ring r=0,x,dp; 148 signsin(11/3); 149 } 150 /////////////////////////////////////////////////////////////////////////// 151 static proc split1(number n) 152 "USAGE: split1(n); n=number 153 RETURN: integral and fractional parts of number n 154 EXAMPLE: example split1; shows an example 155 " 156 { 157 number a,b; 158 int r; 159 a=numerator(n); 160 b=denominator(n); 161 int z=int(number(a)); 162 int y=int(number(b)); 163 r=z mod y; 164 int q=(zr) div y; 165 number n1=q; 166 number n2=nn1; 167 list l=n1,n2; 168 return(l); 169 } 170 example 171 { "EXAMPLE:"; echo = 2; 172 ring r=0,x,dp; 173 split1(11/3); 174 } 175 /////////////////////////////////////////////////////////////////////////// 176 static proc sin( int n) 177 "USAGE: sin(n); n=integer 178 RETURN: approximate value of the sine of n by using maclaurin series 179 EXAMPLE: example sin; shows an example 180 " 181 { 182 def R=basering; 183 ring S=0,x,dp; 184 int i; 185 poly f; 186 int z; 187 for(i=1;i<=n;i=i+2) 188 { 189 f=f+(1)^z*x^i/factorial(i) ; 190 z++; 191 } 192 setring R; 193 map phi=S,var(1); 194 poly f=phi(f); 195 return(f); 196 } 197 example 198 { "EXAMPLE:"; echo = 2; 199 ring r=0,x,dp; 200 sin(10); 201 } 202 /////////////////////////////////////////////////////////////////////////// 203 static proc cos(int n) 204 "USAGE: cos(n); n=integer 205 RETURN: approximate value of the cosine of n by using maclaurin series 206 EXAMPLE: example cos; shows an example 207 " 208 { 209 def R=basering; 210 ring S=0,x,dp; 211 poly f; 212 int i; 213 int z; 214 for(i=0;i<=n;i=i+2) 215 { 216 // f=f+(1)^z*x^i/factorial(i,0); 217 f=f+(1)^z*x^i/factorial(i) ; 218 z++; 219 } 220 setring R; 221 map phi=S,var(1); 222 poly f=phi(f); 223 return(f); 224 } 225 example 226 { "EXAMPLE:"; echo = 2; 227 ring r=0,x,dp; 228 cos(10); 229 } 230 /////////////////////////////////////////////////////////////////////////// 231 static proc signcos(number n) 232 "USAGE: signcos(n); n=number 233 RETURN: the sign of cosin(n) (cosin(n) = the value of the cosine of n) 234 EXAMPLE: example signcos; shows an example 235 " 236 { 237 def r=basering; 238 int a; 239 int b=10; 240 ring s=(real,10),x,dp; 241 number m=imap(r,n); 242 number pi = number_pi(11); 243 number t=m*pi; 244 poly f=cos(b); 245 poly h=subst(f,x,t); 246 if(h>0){a=1;} 247 if(h<0){a=1;} 248 setring r; 249 return(a); 250 } 251 example 252 { "EXAMPLE:"; echo = 2; 253 ring r=0,x,dp; 254 signcos(11/3); 255 } 256 /////////////////////////////////////////////////////////////////////////// 257 static proc signa( number u) 258 "USAGE: signa(u); u=number 259 RETURN: the signa of a number 260 EXAMPLE: example signa; shows an example 261 " 262 { 263 list l=split1(u); 264 int z; 265 if( l[1] mod 2==0 ) 266 { z=signsin(l[2]); } 267 else 268 { z=signcos(l[1]);} 269 return(z); 270 } 271 example 272 { "EXAMPLE:"; echo = 2; 273 ring r=0,x,dp; 274 signa(11/3); 275 } 276 /////////////////////////////////////////////////////////////////////////// 277 static proc prods(list L) 278 "USAGE: product(L); L=list of intvec 279 RETURN: product of first components of Newton pairs in L 280 EXAMPLE: example product; shows an example 281 " 282 { 283 int a=L[2][1]; 284 int i; 285 for(i=1;i<=size(L)2;i++) 286 { 287 a=a*L[i+2][1]; 288 } 289 return(a); 290 } 291 example 292 { "EXAMPLE:"; echo = 2; 293 list L=intvec(2,3),intvec(2,1); 294 prods(L); 295 } 296 /////////////////////////////////////////////////////////////////////////// 297 static proc Generalcase(int N, list L) 298 "USAGE: Generalcase(N,f);N=integer,list L of intvec 299 RETURN: signature of surface singularity with Newton pairs in L 300 ASSUME: number of newton pairs greater than 2 301 EXAMPLE: example Generalcase; shows an example 302 " 303 { 304 int i,j,k,n,m,t,p; 305 int a=L[1][2]; 306 int b=prods(L); 307 int d=gcd(N,b); 308 int q=d*brieskornSign(a,L[1][1],N/d); 309 list A; 310 list B; 311 list S; 312 for(i=1;i<=size(L)1;i++) 313 { 314 a=L[i+1][2]+L[i+1][1]*L[i][1]*a; 315 A[i]=a; 316 S[i]=L[i+1][1]; 317 } 318 for(m=1;m<=size(L)2;m++) 319 { 320 B[m]=L[m+2][1]; 321 } 322 list C; 323 int c=B[size(B)]; 324 C[size(B)]=c; 325 for(j=1;j<=size(B)1;j++) 326 { 327 c=c*B[size(B)j]; 328 C[size(B)j]=c; 329 } 330 list D; 331 D[size(L)1]=1; 332 for(k=1;k<=size(L)2;k++) 333 { 334 D[k]=gcd(N,C[k]); 335 } 336 for(n=1;n<=size(L)1;n++) 337 { 338 t=t+D[n]*brieskornSign(A[n],S[n],N/D[n]); 339 } 340 return(q+t); 341 } 342 example 343 { "EXAMPLE:"; echo = 2; 344 int N=2; 345 list L=intvec(2,3),intvec(2,1),intvec(2,1); 346 Generalcase(N,L); 347 } 348 /////////////////////////////////////////////////////////////////////////// 349 static proc signtwopairs(int N,list L) 350 "USAGE: signtwopairs(N,f);N=integer,L=list of intvec 351 RETURN: signature of surface singularity with Newton pairs in L 352 ASSUME: number of newton pairs equal to 2 353 EXAMPLE: example signtwopairs; shows an example 354 " 355 { 356 int a1,a2,b,d1,q,t; 357 a1=L[1][2]; 358 b=prods(L); 359 d1=gcd(N,b); 360 q=d1*brieskornSign(a1,L[1][1],N/d1); 361 a2=L[2][2]+L[2][1]*L[1][1]*a1; 362 t=brieskornSign(a2,L[2][1],N); 363 return(q+t); 364 } 365 example 366 { "EXAMPLE:"; echo = 2; 367 int N=2; 368 list L=intvec(2,3),intvec(2,1); 369 signtwopairs(N,L); 370 } 371 /////////////////////////////////////////////////////////////////////////// 372 static proc DedekindSum(number b, number c, int a) 151 poly f = x1521x14+8x13y6x1316x12y+20x11y2x12+8x11y36x10y2 152 +24x9y3+4x9y216x8y3+26x7y46x6y4+8x5y5+4x3y6y8; 153 signatureP(N,f); 154 } 155 156 /////////////////////////////////////////////////////////////////////////////// 157 // sigma(z^N + f) in terms of the imbedded resolution graph of f  158 159 static proc dedekindSum(number b, number c, int a) 373 160 { 374 161 number s,d,e; … … 385 172 return(s); 386 173 } 387 /////////////////////////////////////////////////////////////////////////// 174 175 /////////////////////////////////////////////////////////////////////////////// 176 177 static proc isRupture(intvec v) 178 //=== decides whether the exceptional divisor given by the row v in the 179 //=== incidence matrix of the resolution graph intersects at least 3 other divisors 180 { 181 int i,j; 182 for(i=1;i<=size(v);i++) 183 { 184 if(v[i]<0){return(0);} 185 if(v[i]!=0){j++;} 186 } 187 return(j>=4); 188 } 189 190 /////////////////////////////////////////////////////////////////////////////// 191 192 static proc sumExcepDiv(intmat N, list M, int K, int n) 193 //=== computes part of the formulae for eta(g,K), g defining an 194 //=== isolated curve singularity 195 //=== N the incidence matrix of the resolution graph of g 196 //=== M list of total multiplicities 197 //=== n = nrows(N) 198 { 199 int i,j,m,d; 200 for(i=1;i<=n;i++) 201 { 202 if(N[i,i]>0) 203 { 204 m=gcd(K,M[i]); 205 for(j=1;j<=n;j++) 206 { 207 if((i!=j)&&(N[i,j]!=0)) 208 { 209 if(m==1){break;} 210 m=gcd(m,M[j]); 211 } 212 } 213 d=d+m1; 214 } 215 } 216 return(d); 217 } 218 219 /////////////////////////////////////////////////////////////////////////////// 220 221 static proc sumEdges(intmat N, list M, int K, int n) 222 //=== computes part of the formulae for eta(g,K), g defining an 223 //=== isolated curve singularity 224 //=== N the incidence matrix of the resolution graph of g 225 //=== M list of total multiplicities 226 //=== n = nrows(N) 227 { 228 int i,j,d; 229 for(i=1;i<=n1;i++) 230 { 231 for(j=i+1;j<=n;j++) 232 { 233 if(N[i,j]==1) 234 { 235 d=d+gcd(K,gcd(M[i],M[j]))1; 236 } 237 } 238 } 239 return(d); 240 } 241 242 /////////////////////////////////////////////////////////////////////////////// 243 244 static proc etaRes(list L, int K) 245 //=== L total multiplicities 246 //=== etainvariant in terms of the imbedded resolution graph of f 247 { 248 int i,j,d; 249 intvec v; 250 number e; 251 intmat N = L[1]; // incidence matrix of the resolution graph 252 int n = ncols(L[1]); // number of vertices in the resolution graph 253 int a = ncols(L[2]); // number of branches 254 list M; // total multiplicities 255 for(i=1;i<=n;i++) 256 { 257 d=L[2][i,1]; 258 for(j=2;j<=a;j++) 259 { 260 d=d+L[2][i,j]; 261 } 262 if(d==0){d=1;} 263 M[i]=d; 264 } 265 for(i=1;i<=n;i++) 266 { 267 v=N[i,1..n]; 268 if(isRupture(v)) // the divisor intersects more then two others 269 { 270 for(j=1;j<=n;j++) 271 { 272 if((i!=j)&&(v[j]!=0)) 273 { 274 e=e+dedekindSum(M[j],K,M[i]); 275 } 276 } 277 } 278 } 279 if(a==1) 280 { 281 //the irreducible case 282 return(4*e); 283 } 284 return(a1+4*e+sumEdges(N,M,K,n)sumExcepDiv(N,M,K,n)); 285 } 286 287 /////////////////////////////////////////////////////////////////////////////// 288 // sigma(z^N + f) in terms of the spectral pairs of f  289 290 static proc fracPart(number n) 291 //=== computes the fractional part n2 of n 292 //=== i.e. n2 is not in Z but nn2 is in Z 293 { 294 number a,b; 295 int r; 296 a = numerator(n); 297 b = denominator(n); 298 int z = int(number(a)); 299 int y = int(number(b)); 300 r = z mod y; 301 int q = (zr) div y; 302 number n1 = q; 303 number n2 = nn1; 304 return(n2); 305 } 306 307 /////////////////////////////////////////////////////////////////////////////// 308 309 static proc etaSpec(list L, int N) 310 //=== L spectral numbers 311 //=== etainvariant in terms of the spectral pairs of f 312 { 313 int i; 314 number e, h; 315 316 int n = ncols(L[1]); 317 318 if((n mod 2) == 0) 319 // 0 is not a spectral number, thus f is irreducible 320 { 321 for(i = n/2+1; i <= n; i++) 322 { 323 e = e + (1  2 * fracPart(N * number(L[1][i]))) * L[3][i]; 324 } 325 return(2*e); 326 } 327 else 328 // 0 is a spectral number, thus f is reducible 329 { 330 // sum of Hodge numbers in eta function 331 for(i = 1; i <= n; i++) 332 { 333 if((L[2][i] == 2)&&((denominator(leadcoef(N*L[1][i]))==1)(denominator(leadcoef(N*L[1][i]))==1))) 334 { 335 h = h + L[3][i]; 336 } 337 } 338 339 // summand coming from spectral number 0 in eta function 340 h = h + L[3][(n+1)/2]; 341 342 // sum coming from nonzero spectral numbers in eta function 343 for(i = (n+3)/2; i <= n; i++) 344 { 345 if(!((denominator(leadcoef(N*L[1][i]))==1)(denominator(leadcoef(N*L[1][i]))==1))) 346 { 347 e = e + (1  2 * fracPart(N * number(L[1][i]))) * L[3][i]; 348 } 349 } 350 return(h + 2*e); 351 } 352 } 353 354 /////////////////////////////////////////////////////////////////////////////// 355 // Consolidation of the three former variants  356 357 proc signature(int N, poly f, list #) 358 "USAGE: signature(N,f); N = integer, f = reduced poly in 2 variables, # empty or 1,2,3 359 @*  if # is empty or #[1] = 2 then resolution of singularities is used 360 @*  if #[1] = 1 then f has to be analytically irreducible and Puiseux expansions are used 361 @*  if #[1] = 3 then spectral pairs are used 362 RETURN: signature of surface singularity defined by z^N + f(x,y) = 0 363 EXAMPLE: example signature; shows an example 364 " 365 { 366 if(size(#) == 0) 367 { 368 list L = totalmultiplicities(f); 369 return(etaRes(L,N)  N*etaRes(L,1)); 370 } 371 372 if(#[1] == 1) 373 { 374 return(signatureP(N,f)); 375 } 376 377 if(#[1] == 2) 378 { 379 list L = totalmultiplicities(f); 380 return(etaRes(L,N)  N*etaRes(L,1)); 381 } 382 383 if(#[1] == 3) 384 { 385 def R = basering; 386 def Rds = changeord("ds"); 387 setring Rds; 388 poly f = imap(R,f); 389 list L = sppairs(f); 390 setring R; 391 list L = imap(Rds,L); 392 return(etaSpec(L,N)  N*etaSpec(L,1)); 393 } 394 } 395 example 396 { "EXAMPLE:"; echo = 2; 397 ring r = 0,(x,y),dp; 398 int N = 3; 399 poly f = x1521x14+8x13y6x1316x12y+20x11y2x12+8x11y36x10y2 400 +24x9y3+4x9y216x8y3+26x7y46x6y4+8x5y5+4x3y6y8; 401 signature(N,f,1); 402 signature(N,f,2); 403 } 404 405 /////////////////////////////////////////////////////////////////////////////// 388 406 389 407 /* 390 408 Further examples 391 409 392 ring r = 0,(x,y),dp; 393 int N; 394 poly f; 395 396 N = 5; 397 poly f= x1521x14+8x13y6x1316x12y+20x11y2x12+8x11y36x10y2 //3 characteristic pairs 398 +24x9y3+4x9y216x8y3+26x7y46x6y4+8x5y5+4x3y6y8; 399 400 N=6; 401 f= y4+2x3y2+x6+x5y; //2 characteristic pairs 402 403 N=7; 404 f=x5+y11; //1 characteristc pair 410 ring r = 0,(x,y),dp; 411 int N; 412 poly f,g,g1,g2,g3; 413 414 415 // irreducible polynomials 416 417 N = 5; 418 f = x1521x14+8x13y6x1316x12y+20x11y2x12+8x11y36x10y2 419 +24x9y3+4x9y216x8y3+26x7y46x6y4+8x5y5+4x3y6y8; 420 g = f^3 + x17y17; 421 422 N = 6; 423 f = y4+2x3y2+x6+x5y; 424 g1 = f^2 + x5y5; 425 g2 = f^3 + x11y11; 426 g3 = f^3 + x17y17; 427 428 N = 7; 429 f = x5+y11; 430 g1 = f^3 + x11y11; 431 g2 = f^3 + x17y17; 432 433 N = 6; 434 // k0 = 30, k1 = 35, k2 = 71 435 f = x71+6x65+15x59630x52y6+20x53+6230x46y6+910x39y12+15x47 436 7530x40y6+14955x33y12285x26y18+6x41+1230x34y6+4680x27y12 437 +1830x20y18+30x13y24+x355x28y6+10x21y1210x14y18+5x7y24y30; 438 439 // k0 = 16, k1 = 24, k2 = 28, k3 = 30, k4 = 31 440 f = x31781x30+16x29y3010x292464x28y+104x27y22805x287024x27y 441 5352x26y2+368x25y3+366x277136x26y984x25y28000x24y3 442 +836x23y4+34x26320x25y6464x24y2+6560x23y38812x22y4+1392x21y5 443 12x25+256x24y1296x23y21536x22y3+4416x21y48864x20y5+1752x19y6 444 x24+16x23y88x22y216x21y3404x20y4+3056x19y56872x18y6+1648x17y7 445 +8x21y296x20y3+524x19y41472x18y5+3464x17y63808x16y7+1290x15y8 446 28x18y4+240x17y5976x16y6+2208x15y72494x14y8+816x13y9+56x15y6 447 320x14y7+844x13y81216x12y9+440x11y1070x12y8+240x11y9344x10y10 448 +240x9y11+56x9y1096x8y11+52x7y1228x6y12+16x5y13+8x3y14y16; 449 450 451 // reducible polynomials 452 453 N = 12; 454 f = ((y2x3)^2  4x5y  x7)*(x2y3); 455 456 f = 2x3y32y5+x4xy2; 457 458 f = x3y3+x6y+xy6x4y4; 459 */
Note: See TracChangeset
for help on using the changeset viewer.