[24e56cc] | 1 | //GP, last modified 23.10.06 |
---|
| 2 | /////////////////////////////////////////////////////////////////////////////// |
---|
[3eadab] | 3 | version="$Id: modstd.lib,v 1.14 2007-07-20 10:02:38 Singular Exp $"; |
---|
[24e56cc] | 4 | category="Commutative Algebra"; |
---|
| 5 | info=" |
---|
[abb4919] | 6 | LIBRARY: modstd.lib Grobner basis of ideals |
---|
[24e56cc] | 7 | AUTHORS: A. Hashemi, Amir.Hashemi@lip6.fr |
---|
| 8 | @* G. Pfister pfister@mathematik.uni-kl.de |
---|
| 9 | @* H. Schoenemann hannes@mathematik.uni-kl.de |
---|
| 10 | |
---|
| 11 | NOTE: |
---|
[6881a39] | 12 | A library for computing the Grobner basis of an ideal in the polynomial |
---|
[f86ef8] | 13 | ring over the rational numbers using modular methods. The procedures are |
---|
[6881a39] | 14 | inspired by the following paper: |
---|
[24e56cc] | 15 | Elizabeth A. Arnold: |
---|
[abb4919] | 16 | Modular Algorithms for Computing Groebner Bases , Journal of Symbolic |
---|
[24e56cc] | 17 | Computation , April 2003, Volume 35, (4), p. 403-419. |
---|
| 18 | |
---|
| 19 | |
---|
[967543] | 20 | |
---|
[24e56cc] | 21 | PROCEDURES: |
---|
[967543] | 22 | modStd(I); compute a standard basis of I using modular methods |
---|
[a2c2031] | 23 | modS(I,L); liftings to Q of standard bases of I mod p for p in L |
---|
[96b9fdf] | 24 | primeList(n); intvec of n primes <= 2134567879 in decreasing order |
---|
[24e56cc] | 25 | "; |
---|
| 26 | |
---|
[967543] | 27 | LIB "poly.lib"; |
---|
[d6374d] | 28 | LIB "crypto.lib"; |
---|
[24e56cc] | 29 | /////////////////////////////////////////////////////////////////////////////// |
---|
[967543] | 30 | proc modStd(ideal I) |
---|
[a2c2031] | 31 | "USAGE: modStd(I); |
---|
| 32 | RETURN: a standard basis of I if no warning appears; |
---|
| 33 | NOTE: the procedure computes a standard basis of I (over the |
---|
| 34 | rational numbers) by using modular methods. If a |
---|
[967543] | 35 | warning appears then the result is a standard basis with no defined |
---|
[a2c2031] | 36 | relation to I; this is a sign that not enough prime numbers have |
---|
[6881a39] | 37 | been used. For further experiments see procedure modS. |
---|
[24e56cc] | 38 | EXAMPLE: example modStd; shows an example |
---|
| 39 | " |
---|
| 40 | { |
---|
[6881a39] | 41 | def R0=basering; |
---|
| 42 | list rl=ringlist(R0); |
---|
| 43 | if((npars(R0)>0)||(rl[1]>0)) |
---|
| 44 | { |
---|
| 45 | ERROR("characteristic of basering should be zero"); |
---|
| 46 | } |
---|
| 47 | int l,j,k,q; |
---|
[967543] | 48 | int en=2134567879; |
---|
| 49 | int an=1000000000; |
---|
| 50 | intvec hi,hl,hc,hpl,hpc; |
---|
[24e56cc] | 51 | list T,TT; |
---|
[96b9fdf] | 52 | intvec L=primeList(5); |
---|
[967543] | 53 | L[6]=prime(random(an,en)); |
---|
[24e56cc] | 54 | ideal J,cT,lT,K; |
---|
[a2c2031] | 55 | ideal I0=I; |
---|
[967543] | 56 | int h=homog(I); |
---|
| 57 | if((!h)&&(ord_test(R0)==0)) |
---|
| 58 | { |
---|
| 59 | ERROR("input is not homogeneous and ordering is not local"); |
---|
| 60 | } |
---|
| 61 | if(h) |
---|
| 62 | { |
---|
| 63 | execute("ring gn="+string(L[6])+",x(1.."+string(nvars(R0))+"),dp;"); |
---|
| 64 | ideal I=fetch(R0,I); |
---|
| 65 | ideal J=std(I); |
---|
| 66 | hi=hilb(J,1); |
---|
| 67 | setring R0; |
---|
| 68 | } |
---|
[24e56cc] | 69 | for (j=1;j<=size(L);j++) |
---|
| 70 | { |
---|
| 71 | rl[1]=L[j]; |
---|
| 72 | def oro=ring(rl); |
---|
| 73 | setring oro; |
---|
| 74 | ideal I=fetch(R0,I); |
---|
| 75 | option(redSB); |
---|
[967543] | 76 | if(h) |
---|
| 77 | { |
---|
| 78 | ideal I1=std(I,hi); |
---|
| 79 | } |
---|
| 80 | else |
---|
| 81 | { |
---|
| 82 | if(ord_test(R0)==-1) |
---|
| 83 | { |
---|
| 84 | ideal I1=std(I); |
---|
| 85 | } |
---|
| 86 | else |
---|
| 87 | { |
---|
| 88 | matrix M; |
---|
| 89 | ideal I1=liftstd(I,M); |
---|
| 90 | } |
---|
| 91 | } |
---|
[24e56cc] | 92 | setring R0; |
---|
| 93 | T[j]=fetch(oro,I1); |
---|
| 94 | kill oro; |
---|
| 95 | } |
---|
| 96 | //================= delete unlucky primes ==================== |
---|
| 97 | // unlucky iff the leading ideal is wrong |
---|
[967543] | 98 | list LL=deleteUnluckyPrimes(T,L); |
---|
| 99 | T=LL[1]; |
---|
| 100 | L=LL[2]; |
---|
| 101 | lT=LL[3]; |
---|
[24e56cc] | 102 | //============ now all leading ideals are the same ============ |
---|
| 103 | for(j=1;j<=ncols(T[1]);j++) |
---|
| 104 | { |
---|
| 105 | for(k=1;k<=size(L);k++) |
---|
| 106 | { |
---|
| 107 | TT[k]=T[k][j]; |
---|
| 108 | } |
---|
[a2c2031] | 109 | J[j]=liftPoly(TT,L); |
---|
[24e56cc] | 110 | } |
---|
[6881a39] | 111 | //=========== chooses more primes up to the moment the result becomes stable |
---|
[24e56cc] | 112 | while(1) |
---|
| 113 | { |
---|
| 114 | k=0; |
---|
[967543] | 115 | q=prime(random(an,en)); |
---|
[24e56cc] | 116 | while(k<size(L)) |
---|
| 117 | { |
---|
| 118 | k++; |
---|
| 119 | if(L[k]==q) |
---|
| 120 | { |
---|
| 121 | k=0; |
---|
[967543] | 122 | q=prime(random(an,en)); |
---|
[24e56cc] | 123 | } |
---|
| 124 | } |
---|
| 125 | L[size(L)+1]=q; |
---|
| 126 | rl[1]=L[size(L)]; |
---|
| 127 | def @r=ring(rl); |
---|
| 128 | setring @r; |
---|
| 129 | ideal i=fetch(R0,I); |
---|
| 130 | option(redSB); |
---|
[967543] | 131 | if(h) |
---|
| 132 | { |
---|
| 133 | i=std(i,hi); |
---|
| 134 | } |
---|
| 135 | else |
---|
| 136 | { |
---|
| 137 | if(ord_test(R0)==-1) |
---|
| 138 | { |
---|
| 139 | i=std(i); |
---|
| 140 | } |
---|
| 141 | else |
---|
| 142 | { |
---|
| 143 | matrix M; |
---|
| 144 | i=liftstd(i,M); |
---|
| 145 | } |
---|
| 146 | } |
---|
[24e56cc] | 147 | setring R0; |
---|
| 148 | T[size(T)+1]=fetch(@r,i); |
---|
| 149 | kill @r; |
---|
| 150 | cT=lead(T[size(T)]); |
---|
| 151 | attrib(cT,"isSB",1); |
---|
| 152 | if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0)) |
---|
| 153 | { |
---|
| 154 | T=delete(T,size(T)); |
---|
[96b9fdf] | 155 | L=L[1..size(L)-1]; |
---|
[24e56cc] | 156 | k=0; |
---|
| 157 | } |
---|
| 158 | else |
---|
| 159 | { |
---|
| 160 | for(j=1;j<=ncols(T[1]);j++) |
---|
| 161 | { |
---|
| 162 | for(k=1;k<=size(L);k++) |
---|
| 163 | { |
---|
| 164 | TT[k]=T[k][j]; |
---|
| 165 | } |
---|
| 166 | K[j]=liftPoly(TT,L); |
---|
| 167 | } |
---|
| 168 | k=1; |
---|
[967543] | 169 | for(j=1;j<=size(K);j++) |
---|
| 170 | { |
---|
| 171 | if(K[j]-J[j]!=0) |
---|
| 172 | { |
---|
| 173 | k=0; |
---|
| 174 | J=K; |
---|
| 175 | break; |
---|
| 176 | } |
---|
| 177 | } |
---|
[24e56cc] | 178 | } |
---|
| 179 | if(k){break;} |
---|
| 180 | } |
---|
[967543] | 181 | //============ test for standard basis and I=J ======= |
---|
| 182 | J=std(J); |
---|
| 183 | I0=reduce(I0,J); |
---|
[0eb3af] | 184 | if(size(I0)>0) |
---|
| 185 | { |
---|
[a2c2031] | 186 | "WARNING: The input ideal is not contained |
---|
[0eb3af] | 187 | in the ideal generated by the standardbasis"; |
---|
| 188 | "list of primes used:"; |
---|
| 189 | L; |
---|
| 190 | } |
---|
[24e56cc] | 191 | attrib(J,"isSB",1); |
---|
| 192 | return(J); |
---|
| 193 | } |
---|
| 194 | example |
---|
| 195 | { "EXAMPLE:"; echo = 2; |
---|
| 196 | ring r=0,(x,y,z),dp; |
---|
| 197 | ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4; |
---|
| 198 | ideal J=modStd(I); |
---|
| 199 | J; |
---|
[a2c2031] | 200 | } |
---|
[24e56cc] | 201 | /////////////////////////////////////////////////////////////////////////////// |
---|
[4965a1] | 202 | proc modS(ideal I, intvec L, list #) |
---|
| 203 | "USAGE: modS(I,L); I ideal, L intvec of primes |
---|
[967543] | 204 | if size(#)>0 std is used instead of groebner |
---|
[6881a39] | 205 | RETURN: an ideal which is with high probability a standard basis |
---|
| 206 | NOTE: This procedure is designed for fast experiments. |
---|
| 207 | It is not tested whether the result is a standard basis. |
---|
| 208 | It is not tested whether the result generates I. |
---|
[24e56cc] | 209 | EXAMPLE: example modS; shows an example |
---|
| 210 | " |
---|
| 211 | { |
---|
| 212 | int j,k; |
---|
| 213 | list T,TT; |
---|
| 214 | def R0=basering; |
---|
| 215 | ideal J,cT,lT,K; |
---|
| 216 | ideal I0=I; |
---|
| 217 | list rl=ringlist(R0); |
---|
[6881a39] | 218 | if((npars(R0)>0)||(rl[1]>0)) |
---|
| 219 | { |
---|
| 220 | ERROR("characteristic of basering should be zero"); |
---|
| 221 | } |
---|
[24e56cc] | 222 | for (j=1;j<=size(L);j++) |
---|
| 223 | { |
---|
| 224 | rl[1]=L[j]; |
---|
| 225 | def @r=ring(rl); |
---|
| 226 | setring @r; |
---|
| 227 | ideal i=fetch(R0,I); |
---|
| 228 | option(redSB); |
---|
[967543] | 229 | if(size(#)>0) |
---|
| 230 | { |
---|
| 231 | i=std(i); |
---|
| 232 | } |
---|
| 233 | else |
---|
| 234 | { |
---|
[d56af3f] | 235 | i=groebner(i); |
---|
[967543] | 236 | } |
---|
[24e56cc] | 237 | setring R0; |
---|
| 238 | T[j]=fetch(@r,i); |
---|
| 239 | kill @r; |
---|
| 240 | } |
---|
| 241 | //================= delete unlucky primes ==================== |
---|
| 242 | // unlucky iff the leading ideal is wrong |
---|
[967543] | 243 | list LL=deleteUnluckyPrimes(T,L); |
---|
| 244 | T=LL[1]; |
---|
| 245 | L=LL[2]; |
---|
| 246 | //============ now all leading ideals are the same ============ |
---|
| 247 | for(j=1;j<=ncols(T[1]);j++) |
---|
| 248 | { |
---|
| 249 | for(k=1;k<=size(L);k++) |
---|
| 250 | { |
---|
| 251 | TT[k]=T[k][j]; |
---|
| 252 | } |
---|
[a2c2031] | 253 | J[j]=liftPoly(TT,L); |
---|
[967543] | 254 | } |
---|
| 255 | attrib(J,"isSB",1); |
---|
| 256 | return(J); |
---|
| 257 | } |
---|
| 258 | example |
---|
| 259 | { "EXAMPLE:"; echo = 2; |
---|
[4965a1] | 260 | intvec L=3,5,11,13,181; |
---|
[967543] | 261 | ring r=0,(x,y,z),dp; |
---|
| 262 | ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4; |
---|
| 263 | ideal J=modS(I,L); |
---|
| 264 | J; |
---|
[a2c2031] | 265 | } |
---|
| 266 | /////////////////////////////////////////////////////////////////////////////// |
---|
[96b9fdf] | 267 | proc deleteUnluckyPrimes(list T,intvec L) |
---|
| 268 | "USAGE: deleteUnluckyPrimes(T,L);T list of polys, L intvec of primes |
---|
| 269 | RETURN: list L,T with T list of polys, L intvec of primes |
---|
[967543] | 270 | EXAMPLE: example deleteUnluckyPrimes; shows an example |
---|
| 271 | NOTE: works only for homogeneous ideals with global orderings or |
---|
| 272 | for ideals with local orderings |
---|
| 273 | " |
---|
| 274 | { |
---|
| 275 | int j,k; |
---|
| 276 | intvec hl,hc; |
---|
| 277 | ideal cT,lT; |
---|
| 278 | |
---|
| 279 | lT=lead(T[size(T)]); |
---|
| 280 | attrib(lT,"isSB",1); |
---|
| 281 | hl=hilb(lT,1); |
---|
| 282 | for (j=1;j<size(T);j++) |
---|
[24e56cc] | 283 | { |
---|
| 284 | cT=lead(T[j]); |
---|
[967543] | 285 | attrib(cT,"isSB",1); |
---|
| 286 | hc=hilb(cT,1); |
---|
| 287 | if(hl==hc) |
---|
[24e56cc] | 288 | { |
---|
[967543] | 289 | for(k=1;k<=size(lT);k++) |
---|
| 290 | { |
---|
| 291 | if(lT[k]<cT[k]){lT=cT;break;} |
---|
| 292 | if(lT[k]>cT[k]){break;} |
---|
| 293 | } |
---|
| 294 | } |
---|
| 295 | else |
---|
| 296 | { |
---|
| 297 | if(hc<hl){lT=cT;hl=hilb(lT,1);} |
---|
[24e56cc] | 298 | } |
---|
| 299 | } |
---|
| 300 | j=1; |
---|
| 301 | attrib(lT,"isSB",1); |
---|
| 302 | while(j<=size(T)) |
---|
| 303 | { |
---|
| 304 | cT=lead(T[j]); |
---|
| 305 | attrib(cT,"isSB",1); |
---|
| 306 | if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0)) |
---|
| 307 | { |
---|
| 308 | T=delete(T,j); |
---|
[96b9fdf] | 309 | if(j==1) { L=L[2..size(L)]; } |
---|
| 310 | else |
---|
| 311 | { |
---|
| 312 | if (j==size(L)) { L=L[1..size(L)-1]; } |
---|
| 313 | else { L=L[1..j-1],L[j+1..size(L)]; } |
---|
| 314 | } |
---|
[24e56cc] | 315 | j--; |
---|
| 316 | } |
---|
| 317 | j++; |
---|
| 318 | } |
---|
[967543] | 319 | return(list(T,L,lT)); |
---|
[24e56cc] | 320 | } |
---|
| 321 | example |
---|
| 322 | { "EXAMPLE:"; echo = 2; |
---|
[967543] | 323 | list L=2,3,5,7,11; |
---|
| 324 | ring r=0,(y,x),Dp; |
---|
| 325 | ideal I1=y2x,y6; |
---|
| 326 | ideal I2=yx2,y3x,x5,y6; |
---|
| 327 | ideal I3=y2x,x3y,x5,y6; |
---|
| 328 | ideal I4=y2x,x3y,x5; |
---|
| 329 | ideal I5=y2x,yx3,x5,y6; |
---|
| 330 | list T=I1,I2,I3,I4,I5; |
---|
| 331 | list TT=deleteUnluckyPrimes(T,L); |
---|
| 332 | TT; |
---|
[a2c2031] | 333 | } |
---|
[abb4919] | 334 | /////////////////////////////////////////////////////////////////////////////// |
---|
[96b9fdf] | 335 | proc liftPoly(list T, intvec L) |
---|
| 336 | "USAGE: liftPoly(T,L); T list of polys, L intvec of primes |
---|
[24e56cc] | 337 | RETURN: poly p in Q[x] such that p mod L[i]=T[i] |
---|
| 338 | EXAMPLE: example liftPoly; shows an example |
---|
| 339 | " |
---|
[1f535d] | 340 | { |
---|
[e3c414] | 341 | int i; |
---|
| 342 | list TT; |
---|
| 343 | for(i=size(T);i>0;i--) |
---|
[573731] | 344 | { TT[i]=ideal(T[i]); } |
---|
| 345 | T=TT; |
---|
[1f535d] | 346 | ideal hh=chinrem(T,L); |
---|
| 347 | poly h=hh[1]; |
---|
| 348 | poly p=lead(h); |
---|
| 349 | poly result; |
---|
| 350 | number n; |
---|
| 351 | number N=L[1]; |
---|
| 352 | for(i=size(L);i>1;i--) |
---|
| 353 | { |
---|
| 354 | N=N*L[i]; |
---|
| 355 | } |
---|
| 356 | while(h!=0) |
---|
| 357 | { |
---|
| 358 | n=Farey(N,leadcoef(h)); |
---|
| 359 | result=result+n*p; |
---|
| 360 | h=h-lead(h); |
---|
| 361 | p=leadmonom(h); |
---|
| 362 | } |
---|
| 363 | return(result); |
---|
| 364 | } |
---|
| 365 | example |
---|
| 366 | { "EXAMPLE:"; echo = 2; |
---|
| 367 | ring R = 0,(x,y),dp; |
---|
| 368 | intvec L=32003,181,241,499; |
---|
| 369 | list T=ideal(x2+7000x+13000),ideal(x2+100x+147y+40),ideal(x2+120x+191y+10),ideal(x2+x+67y+100); |
---|
| 370 | liftPoly(T,L); |
---|
| 371 | } |
---|
| 372 | /////////////////////////////////////////////////////////////////////////// |
---|
| 373 | proc liftPoly1(list T, intvec L) |
---|
| 374 | "USAGE: liftPoly1(T,L); T list of polys, L intvec of primes |
---|
| 375 | RETURN: poly p in Q[x] such that p mod L[i]=T[i] |
---|
| 376 | EXAMPLE: example liftPoly; shows an example |
---|
| 377 | " |
---|
[24e56cc] | 378 | { |
---|
| 379 | poly result; |
---|
| 380 | int i; |
---|
| 381 | poly p; |
---|
| 382 | list TT; |
---|
| 383 | number n; |
---|
[a2c2031] | 384 | |
---|
[24e56cc] | 385 | number N=L[1]; |
---|
| 386 | for(i=2;i<=size(L);i++) |
---|
| 387 | { |
---|
| 388 | N=N*L[i]; |
---|
| 389 | } |
---|
| 390 | while(1) |
---|
| 391 | { |
---|
| 392 | p=leadmonom(T[1]); |
---|
| 393 | for(i=2;i<=size(T);i++) |
---|
| 394 | { |
---|
| 395 | if(leadmonom(T[i])>p) |
---|
| 396 | { |
---|
| 397 | p=leadmonom(T[i]); |
---|
| 398 | } |
---|
| 399 | } |
---|
| 400 | if (p==0) {return(result);} |
---|
| 401 | for(i=1;i<=size(T);i++) |
---|
| 402 | { |
---|
| 403 | if(p==leadmonom(T[i])) |
---|
| 404 | { |
---|
| 405 | TT[i]=leadcoef(T[i]); |
---|
| 406 | T[i]=T[i]-lead(T[i]); |
---|
| 407 | } |
---|
| 408 | else |
---|
| 409 | { |
---|
| 410 | TT[i]=0; |
---|
| 411 | } |
---|
| 412 | } |
---|
[967543] | 413 | n=chineseR(TT,L,N); |
---|
[a2c2031] | 414 | n=Farey(N,n); |
---|
[24e56cc] | 415 | result=result+n*p; |
---|
| 416 | } |
---|
| 417 | } |
---|
| 418 | example |
---|
| 419 | { "EXAMPLE:"; echo = 2; |
---|
| 420 | ring R = 0,(x,y),dp; |
---|
[4965a1] | 421 | intvec L=32003,181,241,499; |
---|
[24e56cc] | 422 | list T=x2+7000x+13000,x2+100x+147y+40,x2+120x+191y+10,x2+x+67y+100; |
---|
[1f535d] | 423 | liftPoly1(T,L); |
---|
[967543] | 424 | } |
---|
[a2c2031] | 425 | /////////////////////////////////////////////////////////////////////////////// |
---|
[96b9fdf] | 426 | proc fareyIdeal(ideal I,intvec L) |
---|
[967543] | 427 | { |
---|
| 428 | poly result,p; |
---|
| 429 | int i,j; |
---|
| 430 | number n; |
---|
| 431 | number N=L[1]; |
---|
| 432 | for(i=2;i<=size(L);i++) |
---|
| 433 | { |
---|
| 434 | N=N*L[i]; |
---|
| 435 | } |
---|
| 436 | |
---|
| 437 | for(i=1;i<=size(I);i++) |
---|
| 438 | { |
---|
| 439 | p=I[i]; |
---|
| 440 | result=lead(p); |
---|
| 441 | while(1) |
---|
| 442 | { |
---|
| 443 | if (p==0) {break;} |
---|
| 444 | p=p-lead(p); |
---|
[a2c2031] | 445 | n=Farey(N,leadcoef(p)); |
---|
[967543] | 446 | result=result+n*leadmonom(p); |
---|
| 447 | } |
---|
| 448 | I[i]=result; |
---|
| 449 | } |
---|
| 450 | return(I); |
---|
[abb4919] | 451 | } |
---|
[24e56cc] | 452 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 453 | proc Farey (number P, number N) |
---|
| 454 | "USAGE: Farey (P,N); P, N number; |
---|
[a2c2031] | 455 | RETURN: a rational number a/b such that a/b=N mod P |
---|
[24e56cc] | 456 | and |a|,|b|<(P/2)^{1/2} |
---|
| 457 | " |
---|
| 458 | { |
---|
| 459 | if (P<0){P=-P;} |
---|
| 460 | if (N<0){N=N+P;} |
---|
| 461 | number A,B,C,D,E; |
---|
| 462 | E=P; |
---|
| 463 | B=1; |
---|
| 464 | while (N!=0) |
---|
| 465 | { |
---|
[967543] | 466 | if (2*N^2<P) |
---|
[a2c2031] | 467 | { |
---|
[967543] | 468 | return(N/B); |
---|
| 469 | } |
---|
[62a98a] | 470 | D=E mod N; |
---|
| 471 | C=A-(E-E mod N)/N*B; |
---|
| 472 | E=N; |
---|
| 473 | N=D; |
---|
| 474 | A=B; |
---|
| 475 | B=C; |
---|
[24e56cc] | 476 | } |
---|
| 477 | return(0); |
---|
| 478 | } |
---|
| 479 | example |
---|
| 480 | { "EXAMPLE:"; echo = 2; |
---|
| 481 | ring R = 0,x,dp; |
---|
| 482 | Farey(32003,12345); |
---|
[abb4919] | 483 | } |
---|
[967543] | 484 | /////////////////////////////////////////////////////////////////////////////// |
---|
[96b9fdf] | 485 | proc chineseR(list T,intvec L,number N) |
---|
[62a98a] | 486 | "USAGE: chineseR(T,L,N); |
---|
| 487 | RETURN: x such that x = T[i] mod L[i], N=product(L[i]) |
---|
[abb4919] | 488 | NOTE: chinese remainder theorem |
---|
[967543] | 489 | EXAMPLE:example chineseR; shows an example |
---|
[24e56cc] | 490 | " |
---|
| 491 | { |
---|
| 492 | number x; |
---|
| 493 | if(size(L)==1) |
---|
| 494 | { |
---|
| 495 | x=T[1] mod L[1]; |
---|
| 496 | return(x); |
---|
| 497 | } |
---|
| 498 | int i; |
---|
| 499 | int n=size(L); |
---|
| 500 | list M; |
---|
| 501 | for(i=1;i<=n;i++) |
---|
| 502 | { |
---|
| 503 | M[i]=N/L[i]; |
---|
| 504 | } |
---|
| 505 | list S=eexgcdN(M); |
---|
| 506 | for(i=1;i<=n;i++) |
---|
| 507 | { |
---|
| 508 | x=x+S[i]*M[i]*T[i]; |
---|
| 509 | } |
---|
| 510 | x=x mod N; |
---|
| 511 | return(x); |
---|
| 512 | } |
---|
| 513 | example |
---|
| 514 | { "EXAMPLE:"; echo = 2; |
---|
| 515 | ring R = 0,x,dp; |
---|
[96b9fdf] | 516 | chineseR(list(24,15,7),intvec(2,3,5),30); |
---|
[abb4919] | 517 | } |
---|
[a2c2031] | 518 | |
---|
[24e56cc] | 519 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 520 | proc primeList(int n) |
---|
| 521 | "USAGE: primeList(n); |
---|
[96b9fdf] | 522 | RETURN: the intvec of n greatest primes <= 2134567879 |
---|
[24e56cc] | 523 | EXAMPLE:example primList; shows an example |
---|
| 524 | " |
---|
| 525 | { |
---|
[96b9fdf] | 526 | intvec L=0:n; |
---|
[24e56cc] | 527 | int i; |
---|
| 528 | int p=2134567879; |
---|
| 529 | for(i=1;i<=n;i++) |
---|
| 530 | { |
---|
| 531 | L[i]=p; |
---|
| 532 | p=prime(p-1); |
---|
| 533 | } |
---|
| 534 | return(L); |
---|
| 535 | } |
---|
| 536 | example |
---|
| 537 | { "EXAMPLE:"; echo = 2; |
---|
[96b9fdf] | 538 | intvec L=primeList(10); |
---|
[24e56cc] | 539 | size(L); |
---|
| 540 | L[size(L)]; |
---|
[a2c2031] | 541 | } |
---|
[24e56cc] | 542 | /////////////////////////////////////////////////////////////////////////////// |
---|
[f86ef8] | 543 | proc pStd(int p,ideal i) |
---|
| 544 | "USAGE: pStd(p,i);p integer, i ideal; |
---|
| 545 | RETURN: an ideal G which is the groebner base for i |
---|
| 546 | EXAMPLE: example pStd; shows an example |
---|
| 547 | " |
---|
| 548 | { |
---|
| 549 | def r=basering; |
---|
| 550 | list rl=ringlist(r); |
---|
| 551 | rl[1]=p; |
---|
| 552 | def r1=ring(rl); |
---|
| 553 | setring r1; |
---|
| 554 | option(redSB); |
---|
| 555 | ideal j=fetch(r,i); |
---|
| 556 | ideal GP=groebner(j); |
---|
| 557 | setring r; |
---|
| 558 | ideal G=fetch(r1,GP); |
---|
| 559 | attrib(G,"isSB",1); |
---|
| 560 | matrix Z=transmat(p,i,G); |
---|
| 561 | matrix G1=gstrich1(p,Z,i,G); |
---|
| 562 | ideal g1=G1; |
---|
| 563 | ideal g22=reduce(g1,G); |
---|
| 564 | matrix G22=transpose(matrix(g22)); |
---|
| 565 | matrix M=redmat(G,G1,G22); |
---|
| 566 | matrix Z2=-M*Z; |
---|
| 567 | kill r1; |
---|
| 568 | number c=p; |
---|
| 569 | matrix G0=transpose(matrix(G)); |
---|
| 570 | G0= MmodN(G0+ (c)* G22,c^2); |
---|
| 571 | matrix GF=fareyMatrix(G0,c^2); |
---|
| 572 | Z=MmodN(Z+(c)*Z2,c^2); |
---|
| 573 | matrix C=transpose(G); |
---|
| 574 | int n=3; |
---|
| 575 | while(GF<>C) |
---|
| 576 | { |
---|
| 577 | C=GF; |
---|
| 578 | G1= gstrich2(c,Z,i,G0,n); |
---|
| 579 | g1=G1; |
---|
| 580 | g22=reduce(g1,G); |
---|
| 581 | G22=transpose(matrix(g22)); |
---|
| 582 | M=redmat(G,G1,G22); |
---|
| 583 | Z2=-M*Z; |
---|
| 584 | Z=MmodN(Z+(c^(n-1))*Z2,c^n); |
---|
| 585 | G0= MmodN(G0+ (c^(n-1))* G22,c^n); |
---|
| 586 | GF=fareyMatrix(G0,c^n); |
---|
| 587 | n++; |
---|
| 588 | } |
---|
| 589 | return(ideal(GF)); |
---|
| 590 | } |
---|
| 591 | example |
---|
| 592 | { "EXAMPLE:"; echo = 2; |
---|
| 593 | ring r=0,(x,y,z),dp; |
---|
| 594 | ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4; |
---|
| 595 | ideal J=pStd(32003,I); |
---|
| 596 | J; |
---|
| 597 | } |
---|
| 598 | /////////////////////////////////////////////////////////////////////////// |
---|
| 599 | proc transmat(int p,ideal i,ideal G) |
---|
| 600 | "USAGE: transmat(p,I,G); p integer, I,G ideal; |
---|
| 601 | RETURN: the transformationmatrix Z for the ideal i mod p and the groebner base for i mod p |
---|
| 602 | EXAMPLE: example transmit; shows an example |
---|
| 603 | " |
---|
| 604 | { |
---|
| 605 | def r=basering; |
---|
| 606 | int n=nvars(r); |
---|
| 607 | list rl=ringlist(r); |
---|
| 608 | rl[1]=p; |
---|
| 609 | def r1=ring(rl); |
---|
| 610 | setring r1; |
---|
| 611 | ideal i=fetch(r,i); |
---|
| 612 | ideal G=fetch(r,G); |
---|
| 613 | attrib(G,"isSB",1); |
---|
| 614 | ring rhelp=p,x(1..n),dp; |
---|
| 615 | list lhelp=ringlist(rhelp); |
---|
| 616 | list l=lhelp[3]; |
---|
| 617 | setring r; |
---|
| 618 | rl[3]=l; |
---|
| 619 | def r2=ring(rl); |
---|
| 620 | setring r2; |
---|
| 621 | ideal i=fetch(r,i); |
---|
| 622 | option(redSB); |
---|
| 623 | ideal j=std(i); |
---|
| 624 | matrix T=lift(i,j); |
---|
| 625 | setring r1; |
---|
| 626 | matrix T=fetch(r2,T); |
---|
| 627 | ideal j=fetch(r2,j); |
---|
| 628 | matrix M=lift(j,G); |
---|
| 629 | matrix Z=transpose(T*M); |
---|
| 630 | setring r; |
---|
| 631 | matrix Z=fetch(r1,Z); |
---|
| 632 | return(Z); |
---|
| 633 | } |
---|
| 634 | example |
---|
| 635 | { "EXAMPLE:"; echo = 2; |
---|
| 636 | ring r=0,(x,y,z),dp; |
---|
| 637 | ideal i=3x3+x2+1,11y5+y3+2,5z4+z2+4; |
---|
| 638 | ideal g=x3-60x2-60, z4-36z2+37, y5+33y3+66; |
---|
| 639 | int p=181; |
---|
| 640 | matrix Z=transmat(p,i,g); |
---|
| 641 | Z; |
---|
| 642 | } |
---|
| 643 | |
---|
| 644 | /////////////////////////////////////////////////////////////////////////// |
---|
| 645 | proc gstrich1(int p, matrix Z, ideal i, ideal gp) |
---|
| 646 | "USAGE: gstrich1 (p,Z,i,gp); p integer, Z matrix, i,gp ideals; |
---|
| 647 | RETURN: a matrix G such that (Z*F-GP)/p, where F and GP are the matrices of the ideals i and gp |
---|
| 648 | " |
---|
| 649 | { |
---|
| 650 | matrix F=transpose(matrix(i)); |
---|
| 651 | matrix GP=transpose(matrix(gp)); |
---|
| 652 | matrix G=(Z*F-GP)/p; |
---|
| 653 | return(G); |
---|
| 654 | } |
---|
| 655 | /////////////////////////////////////////////////////////////////////////// |
---|
| 656 | proc gstrich2(number p, matrix Z, ideal i, ideal gp, int n) |
---|
| 657 | "USAGE: gstrich2 (p,Z,i,gp,n); p,n integer, Z matrix, i,gp ideals; |
---|
| 658 | RETURN: a matrix G such that (Z*F-GP)/(p^(n-1)), where F and GP are the matrices of the ideals i and gp |
---|
| 659 | " |
---|
| 660 | { |
---|
| 661 | matrix F=transpose(matrix(i)); |
---|
| 662 | matrix GP=transpose(matrix(gp)); |
---|
| 663 | matrix G=(Z*F-GP)/(p^(n-1)); |
---|
| 664 | return(G); |
---|
| 665 | } |
---|
| 666 | /////////////////////////////////////////////////////////////////////////// |
---|
| 667 | proc redmat(ideal i, matrix h, matrix g) |
---|
| 668 | "USAGE: redmat(i,h,g); i ideal , h,g matrices; |
---|
| 669 | RETURN: a matrix M such that i=M*h+g |
---|
| 670 | " |
---|
| 671 | { |
---|
| 672 | matrix c=h-g; |
---|
| 673 | ideal f=transpose(c); |
---|
| 674 | matrix N=lift(i,f); |
---|
| 675 | matrix M=transpose(N); |
---|
| 676 | return(M); |
---|
| 677 | } |
---|
| 678 | /////////////////////////////////////////////////////////////////////////// |
---|
| 679 | proc fareyMatrix(matrix m,number N) |
---|
| 680 | "USAGE: fareyMatrix(m,y); m matrix, y integer; |
---|
| 681 | RETURN: a matrix k of the matrix m with Farey rational numbers a/b as coefficients |
---|
| 682 | EXAMPLE: example fareyMatrix; shows an example |
---|
| 683 | " |
---|
| 684 | { |
---|
| 685 | ideal I=m; |
---|
| 686 | poly result,p; |
---|
| 687 | int i,j; |
---|
| 688 | number n; |
---|
| 689 | for(i=1;i<=size(I);i++) |
---|
| 690 | { |
---|
| 691 | p=I[i]; |
---|
| 692 | result=lead(p); |
---|
| 693 | while(1) |
---|
| 694 | { |
---|
| 695 | if (p==0) {break;} |
---|
| 696 | p=p-lead(p); |
---|
| 697 | n=Farey(N,leadcoef(p)); |
---|
| 698 | result=result+n*leadmonom(p); |
---|
| 699 | } |
---|
| 700 | I[i]=result; |
---|
| 701 | } |
---|
| 702 | matrix k=transpose(I); |
---|
| 703 | return(k); |
---|
| 704 | } |
---|
| 705 | example |
---|
| 706 | {"EXAMPLE:"; echo = 2; |
---|
| 707 | ring r=0,(x,y,z),dp; |
---|
| 708 | matrix m[3][1]=x3+682794673x2+682794673,z4+204838402z2+819353608, y5+186216729y3+372433458; |
---|
| 709 | int p=32003; |
---|
| 710 | matrix b=fareyMatrix(m,p^2); |
---|
| 711 | b; |
---|
| 712 | } |
---|
| 713 | /////////////////////////////////////////////////////////////////////////// |
---|
| 714 | proc MmodN(matrix Z,number N) |
---|
| 715 | "USAGE: MmodN(Z,N);Z matrix, N number; |
---|
| 716 | RETURN: the matrix Z mod N |
---|
| 717 | EXAMPLE: example MmodN; |
---|
| 718 | " |
---|
| 719 | { |
---|
| 720 | int i,j,k; |
---|
| 721 | poly m,p; |
---|
| 722 | number c; |
---|
| 723 | for(i=1;i<=nrows(Z);i++) |
---|
| 724 | { |
---|
| 725 | for(j=1;j<=ncols(Z);j++) |
---|
| 726 | { |
---|
| 727 | for(k=1;k<=size(Z[i,j]);k++) |
---|
| 728 | { |
---|
| 729 | m=leadmonom(Z[i,j][k]); |
---|
| 730 | c=leadcoef(Z[i,j][k]) mod N; |
---|
| 731 | p=p+c*m; |
---|
| 732 | } |
---|
| 733 | Z[i,j]=p; |
---|
| 734 | p=0; |
---|
| 735 | } |
---|
| 736 | } |
---|
| 737 | return(Z); |
---|
| 738 | } |
---|
| 739 | example |
---|
| 740 | { "EXAMPLE:"; echo = 2; |
---|
| 741 | ring r = 0,(x,y,z),dp; |
---|
| 742 | matrix m[3][1]= x3+10668x2+10668, z4-12801z2+12802, y5-8728y3+14547; |
---|
| 743 | number p=32003; |
---|
| 744 | matrix b=MmodN(m,p^2); |
---|
| 745 | b; |
---|
| 746 | } |
---|
| 747 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 748 | /* |
---|
| 749 | ring r=0,(x,y,z),lp; |
---|
| 750 | poly s1 = 5x3y2z+3y3x2z+7xy2z2; |
---|
| 751 | poly s2 = 3xy2z2+x5+11y2z2; |
---|
| 752 | poly s3 = 4xyz+7x3+12y3+1; |
---|
| 753 | poly s4 = 3x3-4y3+yz2; |
---|
| 754 | ideal i = s1, s2, s3, s4; |
---|
| 755 | |
---|
| 756 | ring r=0,(x,y,z),lp; |
---|
| 757 | poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7; |
---|
| 758 | poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y; |
---|
| 759 | poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z; |
---|
| 760 | poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y; |
---|
| 761 | ideal i = s1, s2, s3, s4; |
---|
| 762 | |
---|
| 763 | ring r=0,(x,y,z),lp; |
---|
| 764 | poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz; |
---|
| 765 | poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4; |
---|
| 766 | poly s3 = 8x3 + 12y3 + xz2 + 3; |
---|
| 767 | poly s4 = 7x2y4 + 18xy3z2 + y3z3; |
---|
| 768 | ideal i = s1, s2, s3, s4; |
---|
| 769 | |
---|
| 770 | int n = 6; |
---|
| 771 | ring r = 0,(x(1..n)),lp; |
---|
| 772 | ideal i = cyclic(n); |
---|
| 773 | ring s=0,(x(1..n),t),lp; |
---|
| 774 | ideal i=imap(r,i); |
---|
| 775 | i=homog(i,t); |
---|
| 776 | |
---|
| 777 | ring r=0,(x(1..4),s),(dp(4),dp); |
---|
| 778 | poly s1 =1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4); |
---|
| 779 | poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4); |
---|
| 780 | poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4); |
---|
| 781 | poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4); |
---|
| 782 | poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4); |
---|
| 783 | ideal i = s1, s2, s3, s4, s5; |
---|
| 784 | |
---|
| 785 | ring r=0,(x,y,z),ds; |
---|
| 786 | int a =16; |
---|
| 787 | int b =15; |
---|
| 788 | int c =4; |
---|
| 789 | int t =1; |
---|
| 790 | poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2; |
---|
| 791 | ideal i= jacob(f); |
---|
| 792 | |
---|
| 793 | ring r=0,(x,y,z),ds; |
---|
| 794 | int a =25; |
---|
| 795 | int b =25; |
---|
| 796 | int c =5; |
---|
| 797 | int t =1; |
---|
| 798 | poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2; |
---|
| 799 | ideal i= jacob(f),f; |
---|
| 800 | |
---|
| 801 | ring r=0,(x,y,z),ds; |
---|
| 802 | int a=10; |
---|
| 803 | poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a; |
---|
| 804 | ideal i= jacob(f); |
---|
| 805 | |
---|
| 806 | ring r=0,(x,y,z),ds; |
---|
| 807 | int a =6; |
---|
| 808 | int b =8; |
---|
| 809 | int c =10; |
---|
| 810 | int alpha =5; |
---|
| 811 | int beta= 5; |
---|
| 812 | int t= 1; |
---|
| 813 | poly f =x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)+x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2; |
---|
| 814 | ideal i= jacob(f); |
---|
| 815 | |
---|
| 816 | */ |
---|
| 817 | |
---|
[24e56cc] | 818 | /* |
---|
| 819 | ring r=0,(x,y,z),lp; |
---|
| 820 | poly s1 = 5x3y2z+3y3x2z+7xy2z2; |
---|
| 821 | poly s2 = 3xy2z2+x5+11y2z2; |
---|
| 822 | poly s3 = 4xyz+7x3+12y3+1; |
---|
| 823 | poly s4 = 3x3-4y3+yz2; |
---|
| 824 | ideal i = s1, s2, s3, s4; |
---|
| 825 | |
---|
| 826 | ring r=0,(x,y,z),lp; |
---|
| 827 | poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7; |
---|
| 828 | poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y; |
---|
| 829 | poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z; |
---|
| 830 | poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y; |
---|
| 831 | ideal i = s1, s2, s3, s4; |
---|
| 832 | |
---|
| 833 | ring r=0,(x,y,z),lp; |
---|
| 834 | poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz; |
---|
| 835 | poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4; |
---|
| 836 | poly s3 = 8x3 + 12y3 + xz2 + 3; |
---|
| 837 | poly s4 = 7x2y4 + 18xy3z2 + y3z3; |
---|
| 838 | ideal i = s1, s2, s3, s4; |
---|
| 839 | |
---|
[967543] | 840 | int n = 6; |
---|
| 841 | ring r = 0,(x(1..n)),lp; |
---|
| 842 | ideal i = cyclic(n); |
---|
| 843 | ring s=0,(x(1..n),t),lp; |
---|
| 844 | ideal i=imap(r,i); |
---|
| 845 | i=homog(i,t); |
---|
[24e56cc] | 846 | |
---|
| 847 | ring r=0,(x(1..4),s),(dp(4),dp); |
---|
| 848 | poly s1 =1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4); |
---|
| 849 | poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4); |
---|
| 850 | poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4); |
---|
| 851 | poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4); |
---|
| 852 | poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4); |
---|
| 853 | ideal i = s1, s2, s3, s4, s5; |
---|
[967543] | 854 | |
---|
| 855 | ring r=0,(x,y,z),ds; |
---|
| 856 | int a =16; |
---|
| 857 | int b =15; |
---|
| 858 | int c =4; |
---|
| 859 | int t =1; |
---|
| 860 | poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2; |
---|
| 861 | ideal i= jacob(f); |
---|
| 862 | |
---|
| 863 | ring r=0,(x,y,z),ds; |
---|
| 864 | int a =25; |
---|
| 865 | int b =25; |
---|
| 866 | int c =5; |
---|
| 867 | int t =1; |
---|
| 868 | poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2; |
---|
| 869 | ideal i= jacob(f),f; |
---|
| 870 | |
---|
| 871 | ring r=0,(x,y,z),ds; |
---|
| 872 | int a=10; |
---|
| 873 | poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a; |
---|
| 874 | ideal i= jacob(f); |
---|
| 875 | |
---|
| 876 | ring r=0,(x,y,z),ds; |
---|
| 877 | int a =6; |
---|
| 878 | int b =8; |
---|
| 879 | int c =10; |
---|
| 880 | int alpha =5; |
---|
| 881 | int beta= 5; |
---|
| 882 | int t= 1; |
---|
| 883 | poly f =x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)+x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2; |
---|
| 884 | ideal i= jacob(f); |
---|
| 885 | |
---|
[24e56cc] | 886 | */ |
---|
[3eadab] | 887 | |
---|