Changeset 204cf0 in git
 Timestamp:
 Apr 13, 2016, 8:26:20 PM (7 years ago)
 Branches:
 (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
 Children:
 5a4188b414ab493b3e3e697845d6e6e6f3e783e8
 Parents:
 b1cef7673ad345288c0f013557cd800eaed8d1aa
 Files:

 7 added
 2 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/nfmodstd.lib
rb1cef7 r204cf0 13 13 A library for computing the Groebner basis of an ideal in the polynomial 14 14 ring over an algebraic number field Q(t) using the modular methods, where t is 15 algebraic over the field of rational numbers Q. For the case Q(t) = Q, the procedure16 is inspired by Arnold [1]. This idea is then extended17 t o the case tnot in Q using factorization as follows:15 algebraic over the field of rational numbers Q. For the case Q(t) = Q, the 16 procedure is inspired by Arnold [1]. This idea is then extended to the case 17 t not in Q using factorization as follows: 18 18 Let f be the minimal polynomial of t. 19 For I, I' ideals in Q(t)[X], Q[X,t]/<f> respectively, we map I to I' via the map sending 20 t to t + <f>. 21 We first choose a prime p such that f has at least two factors in characteristic p and 22 add each factor f_i to I' to obtain the ideal J'_i = I' + <f_i>. 23 We then compute a standard basis G'_i of J'_i for each i and combine the G'_i to G_p 24 (a standard basis of I'_p) using chinese remaindering for polynomials. The procedure is 25 repeated for many primes p, where we compute the G_p in parallel until the 26 number of primes is sufficiently large to recover the correct standard basis G' of I'. 27 Finally, by mapping G' back to Q(t)[X], a standard basis G of I is obtained. 19 For I, I' ideals in Q(t)[X], Q[X,t]/<f> respectively, we map I to I' via the 20 map sending t to t + <f>. We first choose a prime p such that f has at least 21 two factors in characteristic p and add each factor f_i to I' to obtain the 22 ideal J'_i = I' + <f_i>. We then compute a standard basis G'_i of J'_i for each 23 i and combine the G'_i to G_p (a standard basis of I'_p) using chinese remaindering 24 for polynomials. The procedure is repeated for many primes p, where we compute the 25 G_p in parallel until the number of primes is sufficiently large to recover the 26 correct standard basis G' of I'. Finally, by mapping G' back to Q(t)[X], a standard 27 basis G of I is obtained. 28 @*The procedure also works if the input is a module. For this, we consider the 29 rings A = Q(t)[X] and A' = (Q[t]/<f>)[X]. For submodules I, I' in A^m, A'^m, 30 respectively, we map I to I' via the map sending t to t + <f>. As above, we first 31 choose a prime p such that f has at least two factors in characteristic p. For each 32 factor f_{i,p} of f_p := (f mod p), we set I'_{i,p} := (I'_p mod f_{i,p}). We then 33 compute a standard basis G'_i of I'_{i,p} over F_p[t]/<f_{i,p}> for each i and 34 combine the G'_i to G_p (a standard basis of I'_p) using chinese remaindering for 35 polynomials. The procedure is repeated for many primes p as described above and we 36 finally obtain a standard basis of I. 28 37 29 38 REFERENCES: … … 33 42 PROCEDURES: 34 43 chinrempoly(l,m); chinese remaindering for polynomials 35 nfmodStd(I); standard basis of I over algebraic number field using modular methods 44 nfmodStd(I); standard basis of I over algebraic number field using modular 45 methods 36 46 "; 37 47 … … 40 50 //////////////////////////////////////////////////////////////////////////////// 41 51 42 static proc testPrime(int p, ideal I)52 static proc testPrime(int p, list L) 43 53 { 44 54 /* … … 46 56 * and leading coefficients of generating set of ideal 47 57 */ 48 int i,j; 49 poly f; 58 int i,j,k, tmp; 59 def f; 60 def I = L[1]; // list L = def I 50 61 number num; 51 62 bigint d1,d2,d3; 52 for(i = 1; i <= size(I); i++) 63 if(typeof(I)=="ideal") 64 { 65 tmp=1; 66 } 67 for(i = 1; i <= ncols(I); i++) 53 68 { 54 69 f = cleardenom(I[i]); 55 70 if(f == 0) 56 71 { 57 return(0);72 return(0); 58 73 } 59 74 num = leadcoef(I[i])/leadcoef(f); … … 68 83 return(0); 69 84 } 70 for(j = size(f); j > 0; j) 71 { 72 d3 = bigint(leadcoef(f[j])); 73 if( (d3 mod p) == 0) 74 { 75 return(0); 76 } 77 } 85 if(tmp) 86 { 87 for(j = size(f); j > 0; j) 88 { 89 d3 = bigint(leadcoef(f[j])); 90 if( (d3 mod p) == 0) 91 { 92 return(0); 93 } 94 } 95 } 96 else 97 { 98 for(j = nrows(f); j > 0; j) 99 { 100 for(k=1;k<=size(f[j]);k++) 101 { 102 d3 = bigint(leadcoef(f[j][k])); 103 if((d3 mod p) == 0) 104 { 105 return(0); 106 } 107 } 108 } 109 } 110 78 111 } 79 112 return(1); … … 103 136 nr = testfact(ur); 104 137 k=ncols(L[1]); 105 if(nr < k && k < (urnr))// set bound for k138 if(nr < k && k < (urnr))// set a bound for k 106 139 { 107 140 return(1); … … 118 151 { 119 152 /* L=list(I), I=J,f; J ideal , f minpoly */ 120 int sz,nr ,dg;121 idealJ=L[1];153 int sz,nr; 154 def J=L[1]; 122 155 sz=ncols(J); 123 poly f=J[sz]; 124 dg=deg(f); 125 if(!testPrime(p,J) or !minpolyTask(f,p)) 156 def f=J[sz]; 157 poly g; 158 if(typeof(f)=="vector") 159 { 160 g = f[1]; 161 } 162 else 163 { 164 g = f; 165 } 166 if(!testPrime(p,list(J)) or !minpolyTask(g,p)) 126 167 { 127 168 return(0); … … 206 247 static proc chinRm(list m, list ll, list lk,list l1,int uz) 207 248 { 208 poly ff,c; 209 210 for(int i=1;i<=uz;i++) 211 { 212 c = division(l1[i]*ll[i],m[i])[2][1]; 213 ff = ff + c*lk[i]; 214 } 215 return(ff); 249 if(typeof(l1[1])=="ideal" or typeof(l1[1])=="poly") 250 { 251 poly ff,c; 252 for(int i=1;i<=uz;i++) 253 { 254 c = division(l1[i]*ll[i],m[i])[2][1]; 255 ff = ff + c*lk[i]; 256 } 257 return(ff); 258 } 259 else 260 { 261 vector ff,c; 262 for(int i=1;i<=uz;i++) 263 { 264 c = vector(m[i]); 265 attrib(c,"isSB",1); 266 ff = ff + (reduce(l1[i]*ll[i],c))*lk[i]; 267 } 268 return(ff); 269 } 216 270 } 217 271 … … 220 274 proc chinrempoly(list l,list m) 221 275 "USAGE: chinrempoly(l, m); l list, m list 222 RETURN: a polynomial (resp. ideal) which is congruent to l[i] modulo m[i] for all i 276 RETURN: a polynomial (resp. ideal/module) which is congruent to l[i] modulo m[i] 277 for all i 223 278 NOTE: The procedure applies chinese remaindering to the first argument w.r.t. the 224 moduli given in the second. The elements in the first list must be of same type225 which can be polynomial or ideal. The moduli must be of type polynomial. Elements226 in the second list must be distinct and coprime.279 moduli given in the second. The elements in the first list must be of the same 280 type which can be polynomial, ideal, or module. The moduli must be of type 281 polynomial. The elements in the second list must be distinct and coprime. 227 282 SEE ALSO: chinrem 228 283 EXAMPLE: example chinrempoly; shows an example 229 284 " 230 285 { 231 int i,j,sz,uz ;286 int i,j,sz,uz, tmp; 232 287 uz = size(l); 233 sz = ncols(ideal(l[1])); 288 if(typeof(l[1])=="ideal" or typeof(l[1])=="poly") 289 { 290 sz = ncols(ideal(l[1])); 291 tmp = 1; 292 } 293 else 294 { 295 sz = ncols(module(l[1])); 296 } 234 297 poly f=1; 235 298 for(i=1;i<=uz;i++) … … 238 301 } 239 302 240 ideal I,J;241 303 list l1,ll,lk,l2; 242 304 poly c,ff; … … 247 309 } 248 310 249 for(i=1;i<=sz;i++) 250 { 251 for(j=1;j<=uz;j++) 252 { 253 I = l[j]; 254 l1[j] = I[i]; 255 } 256 J[i] = chinRm(m,ll,lk,l1,uz); 257 } 258 return(J); 311 if(tmp) 312 { 313 ideal I,J; 314 for(i=1;i<=sz;i++) 315 { 316 for(j=1;j<=uz;j++) 317 { 318 I = l[j]; 319 l1[j] = I[i]; 320 } 321 J[i] = chinRm(m,ll,lk,l1,uz); 322 } 323 return(J); 324 } 325 else 326 { 327 module I,J; 328 for(i=1;i<=sz;i++) 329 { 330 for(j=1;j<=uz;j++) 331 { 332 I = l[j]; 333 l1[j] = I[i]; 334 } 335 J[i] = chinRm(m,ll,lk,l1,uz); 336 } 337 return(J); 338 } 259 339 } 260 340 example … … 276 356 ring r = p, (x,y,a), (dp(2),dp(1)); 277 357 poly minpolynomial = a^2+1; 278 ideal kf=factorize(minpolynomial,1); //return factors without multiplicity358 ideal kf=factorize(minpolynomial,1); //return factors without multiplicity 279 359 kf; 280 360 ideal k=(a+1)*x2+y, 3xay+ a+2; … … 289 369 I=simplify(I,2); 290 370 I; 291 } 292 371 l = module(k1[2..ncols(k1)]), module(k2[2..ncols(k2)]); 372 module M = chinrempoly(l,m); 373 M; 374 } 293 375 //////////////////////////////////////////////////////////////////////////////// 294 376 … … 300 382 * size(L)>=2 301 383 */ 302 idealJ=L[1];384 def J=L[1]; 303 385 int i=size(L); 304 386 int sc=ncols(J); 305 387 int j,k; 306 polyg=leadmonom(J[1]);388 def g=leadmonom(J[1]); 307 389 for(j=1;j<=i;j++) 308 390 { … … 327 409 //////////////////////////////////////////////////////////////////////////////// 328 410 329 static proc LiftPolyCRT( idealI)411 static proc LiftPolyCRT(def I) 330 412 { 331 413 /* … … 333 415 * to modulo minpoly via CRT for poly over char p>0 334 416 */ 417 def sl; 335 418 int u,in,j; 336 list LL,Lk; 337 ideal J,K,II; 338 poly f; 339 u=ncols(I); 340 J=I[1..u1]; 341 f=I[u]; 342 K=factorize(f,1); 343 in=ncols(K); 344 for(j=1;j<=in;j++) 345 { 346 LL[j]=K[j]; 347 ideal I(j)=J,K[j]; 348 I(j)=std(I(j)); 349 if(size(I(j))==1) 350 { 351 Lk[j]=I(j); 419 list LL,Lk,T2; 420 if(typeof(I)=="ideal") 421 { 422 423 ideal J,K,II; 424 poly f; 425 u=ncols(I); 426 J=I[1..u1]; 427 f=I[u]; 428 K=factorize(f,1); 429 in=ncols(K); 430 for(j=1;j<=in;j++) 431 { 432 LL[j]=K[j]; 433 ideal I(j)=J,K[j]; 434 I(j)=std(I(j)); 435 if(size(I(j))==1) 436 { 437 Lk[j]=I(j); 438 } 439 else 440 { 441 I(j)[1]=0; 442 I(j)=simplify(I(j), 2); 443 Lk[j]=I(j); 444 } 445 } 446 if(check_leadmonom_and_size(Lk)) 447 { 448 // apply CRT for polynomials 449 II =chinrempoly(Lk,LL),f; 352 450 } 353 451 else 354 452 { 355 I(j)[1]=0; 356 I(j)=simplify(I(j), 2); 357 Lk[j]=I(j); 358 } 359 } 360 if(check_leadmonom_and_size(Lk)) 361 { 362 // apply CRT for polynomials 363 II =chinrempoly(Lk,LL),f; 364 } 365 else 366 { 367 II=0; 368 } 369 return(II); 370 } 371 372 //////////////////////////////////////////////////////////////////////////////// 373 374 static proc PtestStd(string command, list args, ideal result, int p) 453 II=0; 454 } 455 return(II); 456 } 457 else 458 { 459 module J,II; 460 vector f; 461 u=ncols(I); 462 J=I[1..u1]; 463 f=I[u]; 464 poly ff = f[1]; 465 ideal K=factorize(ff,1); 466 in=ncols(K); 467 def Ls = basering; 468 list l = ringlist(Ls); 469 if(l[3][1][1]=="c") 470 { 471 l[1] = list(l[1]) + list(list(l[2][size(l[2])])) + 472 list(list(l[3][size(l[3])]))+list(ideal(0)); 473 l[2] = delete(l[2],size(l[2])); 474 l[3] = delete(l[3],size(l[3])); 475 } 476 else 477 { 478 l[1] = list(l[1]) + list(list(l[2][size(l[2])])) + 479 list(list(l[3][size(l[3])1]))+list(ideal(0)); 480 l[2] = delete(l[2],size(l[2])); 481 l[3] = delete(l[3],size(l[3])1); 482 } 483 484 def S1 = ring(l); 485 setring S1; 486 number Num= number(imap(Ls,ff)); 487 list l = ringlist(S1); 488 l[1][4][1] = Num; 489 S1 = ring(l); 490 setring S1; 491 ideal K = imap(Ls,K); 492 def S2; 493 module II; 494 number Num; 495 /* ++++++ if minpoly is irreducible then K will be the zero ideal +++ */ 496 if(size(K)==0) 497 { 498 module M = std(imap(Ls,J)); 499 if(size(M)==1 && M[1]==[1]) 500 { 501 setring Ls; 502 return(module([1])); 503 } 504 II = normalize(M); 505 } 506 else 507 { 508 for(j=1;j<=in;j++) 509 { 510 LL[j]=K[j]; 511 Num = number(K[j]); 512 T2 = ringlist(S1); 513 T2[1][4][1] = Num; 514 S2 = ring(T2); 515 setring S2; 516 module M = std(imap(Ls,J)); 517 if(size(M)== 1 && M[1]==[1]) 518 { 519 setring Ls; 520 return(module([1])); 521 } 522 setring S1; 523 Lk[j]= imap(S2,M); 524 } 525 526 if(check_leadmonom_and_size(Lk)) 527 { 528 // apply CRT for polynomials 529 setring Ls; 530 II =chinrempoly(imap(S1,Lk),imap(S1,LL)); 531 setring S1; 532 II = normalize(imap(Ls,II)); 533 } 534 else 535 { 536 setring S1; 537 II=[0]; 538 } 539 } 540 setring Ls; 541 return(imap(S1,II)); 542 } 543 } 544 545 //////////////////////////////////////////////////////////////////////////////// 546 547 /* test if 'result' is a GB of the input ideal */ 548 static proc final_Test_minpolyzero(string command, alias list args, module result) 549 { 550 int i; 551 list arg = args; 552 attrib(result, "isSB", 1); 553 for (i = ncols(args[1]); i > 0; i) 554 { 555 if (reduce(args[1][i], result, 1) != 0) 556 { 557 return(0); 558 } 559 } 560 /* test if result is a GB */ 561 module G = std(result); 562 if (size(reduce(G, result,1))!=0) //Modstd::reduce_parallel(G, result) 563 { 564 return(0); 565 } 566 return(1); 567 } 568 569 //////////////////////////////////////////////////////////////////////////////// 570 571 /* test if 'result' is a GB of the input ideal */ 572 static proc final_Test(string command, alias list args, def result) 573 { 574 int i; 575 list arg = args; 576 // modified for module case 577 if(typeof(args[1])=="ideal") 578 { 579 /* test if args[1] is in result */ 580 attrib(result, "isSB", 1); 581 for (i = ncols(args[1]); i > 0; i) 582 { 583 if (reduce(args[1][i], result, 1) != 0) 584 { 585 return(0); 586 } 587 } 588 589 /* test if result is a GB */ 590 ideal G = std(result); 591 if (size(reduce(G, result,1))!=0) //Modstd::reduce_parallel(G, result) 592 { 593 return(0); 594 } 595 return(1); 596 } 597 else 598 { 599 /* test if args[1] is in result */ 600 def Ls = basering; 601 list l = ringlist(Ls); 602 if(l[3][1][1]=="c") 603 { 604 l[1] = list(l[1]) + list(list(l[2][size(l[2])])) + 605 list(list(l[3][size(l[3])]))+list(ideal(0)); 606 l[2] = delete(l[2],size(l[2])); 607 l[3] = delete(l[3],size(l[3])); 608 } 609 else 610 { 611 l[1] = list(l[1]) + list(list(l[2][size(l[2])])) + 612 list(list(l[3][size(l[3])1]))+list(ideal(0)); 613 l[2] = delete(l[2],size(l[2])); 614 l[3] = delete(l[3],size(l[3])1); 615 } 616 def sL = ring(l); 617 kill l; 618 setring sL; 619 list arg = imap(Ls,arg); 620 module arg_s =arg[1]; 621 list l = ringlist(sL); 622 l[1][4][1] = arg_s[ncols(arg_s)][1]; 623 arg_s = arg_s[1..ncols(arg_s)1]; 624 def cL = ring(l); 625 setring cL; 626 module ar_gs = imap(sL,arg_s); 627 def Result = imap(Ls,result); 628 attrib(Result, "isSB", 1); 629 for (i = ncols(ar_gs); i > 0; i) 630 { 631 if (reduce(ar_gs[i], Result, 1) != 0) 632 { 633 setring Ls; 634 return(0); 635 } 636 } 637 // test if result is a GB 638 module G = std(Result); 639 if (size(reduce(G,Result,1))!=0) //Modstd::reduce_parallel(G, Result) 640 { 641 setring Ls; 642 return(0); 643 } 644 setring Ls; 645 return(1); 646 } 647 } 648 649 //////////////////////////////////////////////////////////////////////////////// 650 651 static proc PtestStd_minpolyzero(string command, list args, module result, int p) 375 652 { 376 653 /* 377 654 * let G be std of I which is not yet known whether it is the correct 378 * standard basis or not. So this procedure does the first test655 * standard basis. So this procedure does the first test 379 656 */ 380 657 def br = basering; … … 390 667 def rp = ring(lbr); 391 668 setring(rp); 392 ideal Ip = fetch(br, args)[1]; 393 ideal Gp = fetch(br, result); 669 module Ip = imap(br, args)[1]; 670 int i; 671 module Gp = imap(br, result); 394 672 attrib(Gp, "isSB", 1); 395 int i;396 673 for (i = ncols(Ip); i > 0; i) 397 674 { … … 402 679 } 403 680 } 404 Ip = LiftPolyCRT(Ip);681 Ip = std(Ip); 405 682 attrib(Ip,"isSB",1); 406 683 for (i = ncols(Gp); i > 0; i) … … 418 695 //////////////////////////////////////////////////////////////////////////////// 419 696 420 static proc cleardenomIdeal(ideal I) 697 static proc PtestStd(string command, list args, def result, int p) 698 { 699 /* 700 * let G be std of I which is not yet known whether it is the correct 701 * standard basis. So this procedure does the first test 702 */ 703 def br = basering; 704 705 list lbr = ringlist(br); 706 if (typeof(lbr[1]) == "int") 707 { 708 lbr[1] = p; 709 } 710 else 711 { 712 lbr[1][1] = p; 713 } 714 def rp = ring(lbr); 715 setring(rp); 716 def Ip = imap(br, args)[1]; 717 718 int u,in,j,i; 719 list LL,Lk,T2; 720 if(typeof(Ip)!="ideal") 721 { 722 module J,II; 723 vector f; 724 u=ncols(Ip); 725 J=Ip[1..u1]; 726 f=Ip[u]; 727 poly ff = f[1]; 728 ideal K=factorize(ff,1); 729 in=ncols(K); 730 def Ls = basering; 731 list l = ringlist(Ls); 732 if(l[3][1][1]=="c") 733 { 734 l[1] = list(l[1]) + list(list(l[2][size(l[2])])) + 735 list(list(l[3][size(l[3])]))+list(ideal(0)); 736 l[2] = delete(l[2],size(l[2])); 737 l[3] = delete(l[3],size(l[3])); 738 } 739 else 740 { 741 l[1] = list(l[1]) + list(list(l[2][size(l[2])])) + 742 list(list(l[3][size(l[3])1]))+list(ideal(0)); 743 l[2] = delete(l[2],size(l[2])); 744 l[3] = delete(l[3],size(l[3])1); 745 } 746 def S1 = ring(l); 747 setring S1; 748 number Num= number(imap(Ls,ff)); 749 list l = ringlist(S1); 750 l[1][4][1] = Num; 751 S1 = ring(l); 752 setring S1; 753 ideal K = imap(Ls,K); 754 module Jp = imap(Ls,J); 755 def S2; 756 module Ip; 757 number Num; 758 /* ++++++ if the minpoly is irreducible then K = ideal(0) +++ */ 759 if(size(K)==0) 760 { 761 module M = std(Jp); 762 Ip = normalize(M); 763 } 764 else 765 { 766 for(j=1;j<=ncols(K);j++) 767 { 768 LL[j]=K[j]; 769 Num = number(K[j]); 770 T2 = ringlist(S1); 771 T2[1][4][1] = Num; 772 S2 = ring(T2); 773 setring S2; 774 module M = std(imap(Ls,J)); 775 setring S1; 776 Lk[j]= imap(S2,M); 777 } 778 if(check_leadmonom_and_size(Lk)) 779 { 780 // apply CRT for polynomials 781 setring Ls; 782 II =chinrempoly(imap(S1,Lk),imap(S1,LL)); 783 setring S1; 784 Ip = normalize(imap(Ls,II)); 785 } 786 else 787 { 788 setring S1; 789 Ip=[0]; 790 } 791 } 792 setring S1; 793 module Gp = imap(br, result); 794 attrib(Gp, "isSB", 1); 795 for (i = ncols(Jp); i > 0; i) 796 { 797 if (reduce(Jp[i], Gp, 1) != 0) 798 { 799 setring(br); 800 return(0); 801 } 802 } 803 804 attrib(Ip,"isSB",1); 805 for (i = ncols(Gp); i > 0; i) 806 { 807 if (reduce(Gp[i], Ip, 1) != 0) 808 { 809 setring(br); 810 return(0); 811 } 812 } 813 setring(br); 814 return(1); 815 } 816 else 817 { 818 ideal Gp = imap(br, result); 819 attrib(Gp, "isSB", 1); 820 for (i = ncols(Ip); i > 0; i) 821 { 822 if (reduce(Ip[i], Gp, 1) != 0) 823 { 824 setring(br); 825 return(0); 826 } 827 } 828 Ip = LiftPolyCRT(Ip); 829 attrib(Ip,"isSB",1); 830 for (i = ncols(Gp); i > 0; i) 831 { 832 if (reduce(Gp[i], Ip, 1) != 0) 833 { 834 setring(br); 835 return(0); 836 } 837 } 838 setring(br); 839 return(1); 840 } 841 } 842 843 //////////////////////////////////////////////////////////////////////////////// 844 845 static proc cleardenomIdeal(def I) 421 846 { 422 847 int t=ncols(I); … … 437 862 //////////////////////////////////////////////////////////////////////////////// 438 863 439 static proc modStdparallelized( ideal I)864 static proc modStdparallelized(def I, list #) 440 865 { 441 866 // apply modular.lib … … 443 868 intvec opt = option(get); 444 869 option(redSB); 445 I = modular("Nfmodstd::LiftPolyCRT", list(I), PrimeTestTask, Modstd::deleteUnluckyPrimes_std, 446 PtestStd, Modstd::finalTest_std,536870909); 870 if(size(#)>0) 871 { 872 I = modular("std", list(I), testPrime, Modstd::deleteUnluckyPrimes_std, 873 PtestStd_minpolyzero, final_Test_minpolyzero,536870909); 874 } 875 else 876 { 877 I = modular("Nfmodstd::LiftPolyCRT", list(I), PrimeTestTask, 878 Modstd::deleteUnluckyPrimes_std,PtestStd, final_Test,536870909); 879 } 447 880 attrib(I, "isSB", 1); 448 881 option(set,opt); … … 452 885 //////////////////////////////////////////////////////////////////////////////// 453 886 /* main procedure */ 454 proc nfmodStd( idealI, list #)455 "USAGE: nfmodStd(I, #); I ideal , # optional parameters887 proc nfmodStd(def I, list #) 888 "USAGE: nfmodStd(I, #); I ideal or module, # optional parameters 456 889 RETURN: standard basis of I over algebraic number field 457 890 NOTE: The procedure passes to @ref{modStd} if the ground field has no … … 462 895 " 463 896 { 464 list L=#; 465 def Rbs=basering; 466 poly f; 467 ideal J; 468 int n=nvars(Rbs); 469 if(size(I)==0) 470 { 471 return(ideal(0)); 472 } 473 if(npars(Rbs)==0) 474 { 475 J=modStd(I,L);//if algebraic number is in Q 476 if(size(#)>0) 477 { 478 return(cleardenomIdeal(J)); 479 } 480 return(J); 481 } 482 list rl=ringlist(Rbs); 483 f=rl[1][4][1]; 484 rl[2][n+1]=rl[1][2][1]; 485 rl[1]=rl[1][1]; 486 rl[3][size(rl[3])+1]=rl[3][size(rl[3])]; 487 rl[3][size(rl[3])1]=list("dp",1); 488 def S=ring(rl); 489 setring S; 490 poly f=imap(Rbs,f); 491 ideal I=imap(Rbs,I); 492 I = simplify(I,2);// eraze the zero generatos 493 ideal J; 494 if(f==0) 495 { 496 ERROR("minpoly must be nonzero"); 497 } 498 I=I,f; 499 J=modStdparallelized(I); 500 setring Rbs; 501 J=imap(S,J); 502 J=simplify(J,2); 503 if(size(#)>0) 504 { 505 return(cleardenomIdeal(J)); 506 } 507 return(J); 897 list L=#; 898 def Rbs=basering; 899 poly f; 900 int tmp=1; 901 if(typeof(I)!="ideal") 902 { 903 tmp =0; 904 } 905 int n=nvars(Rbs); 906 if(size(I)==0) 907 { 908 if(!tmp) 909 { 910 return(module([0])); 911 } 912 return(ideal(0)); 913 } 914 if(npars(Rbs)==0) 915 { 916 //if algebraic number is in Q 917 if(typeof(I)=="module") 918 { 919 def J = modStdparallelized(I,1); 920 return(J); 921 } 922 else 923 { 924 def J=modStd(I,L); 925 return(J); 926 } 927 } 928 def S; 929 list rl=ringlist(Rbs); 930 f=rl[1][4][1]; 931 if(tmp) 932 { 933 rl[2][n+1]=rl[1][2][1]; 934 rl[1]=rl[1][1]; 935 rl[3][size(rl[3])+1]=rl[3][size(rl[3])]; 936 rl[3][size(rl[3])1]=list("dp",1); 937 } 938 else 939 { 940 if(rl[3][1][1]!="c") 941 { 942 rl[2] = rl[2] + rl[1][2]; 943 rl[3] = insert(rl[3], rl[1][3][1],1); 944 rl[1] = rl[1][1]; 945 } 946 else 947 { 948 rl[2] = rl[2] + rl[1][2]; 949 rl[3][size(rl[3])+1] = rl[1][3][1]; 950 rl[1] = rl[1][1]; 951 } 952 } 953 S = ring(rl); 954 setring S; 955 poly f=imap(Rbs,f); 956 def I=imap(Rbs,I); 957 I = simplify(I,2); // eraze the zero generatos 958 if(f==0) 959 { 960 ERROR("minpoly must be nonzero"); 961 } 962 I=I,f; 963 def J_I=modStdparallelized(I); 964 setring Rbs; 965 def J=imap(S,J_I); 966 J=simplify(J,2); 967 return(J); 508 968 } 509 969 example 510 970 { "EXAMPLE:"; echo = 2; 511 ring r1 = (0,a),(x,y),dp;512 minpoly = a^2+1;513 ideal k =(a/2+1)*x^2+2/3y, 3*xa*y+ a/7+2;514 ideal I =nfmodStd(k);971 ring r1 = (0,a),(x,y),dp; 972 minpoly = a^2+1; 973 ideal k = (a/2+1)*x^2+2/3y, 3*xa*y+ a/7+2; 974 ideal I = nfmodStd(k); 515 975 I; 516 ring r2 =(0,a),(x,y,z),dp; 517 minpoly =a^3 +2; 518 ideal k=(a^2+a/2)*x^2+(a^2 2/3*a)*yz, (3*a^2+1)*zx(a+4/7)*y+ a+2/5; 519 ideal IJ=nfmodStd(k); 976 ring rm = (0,a),(x,y),(c,dp); 977 minpoly = a^3+2a+7; 978 module M = [(a/2+1)*x^2+2/3y, 3*xa*y+ a/7+2], [ax, y]; 979 M = nfmodStd(M); 980 M; 981 ring r2 = (0,a),(x,y,z),dp; 982 minpoly = a^3 +2; 983 ideal k = (a^2+a/2)*x^2+(a^2 2/3*a)*yz, (3*a^2+1)*zx(a+4/7)*y+ a+2/5; 984 ideal IJ = nfmodStd(k); 520 985 IJ; 521 ring r3 =0,(x,y),dp;// ring without parameter986 ring r3 = 0, (x,y), dp; // ring without parameter 522 987 ideal I = x2 + y, xy  7y + 2x; 988 ideal J1 = nfmodStd(I); 989 J1; 990 module J2 = nfmodStd(module(I)); 991 J2; 992 ring r4 = 0, (x,y), (c,dp); 993 module I = [x2, xy], [xy,0], [0,7y + 2x]; 523 994 I=nfmodStd(I); 524 995 I; 525 996 } 526 527 //////////////////////////////////////////////////////////////////////////////528 529 /*530 Benchmark Problems from531 532 Boku, Decker, Fieker, Steenpass: Groebner Bases over Algebraic Number533 Fields.534 535 // 1536 ring R = (0,a), (x,y,z), dp;537 minpoly = (a^2+1);538 poly f1 = (a+8)*x^2*y^2+5*x*y^3+(a+3)*x^3*z539 +x^2*y*z;540 poly f2 = x^5+2*y^3*z^2+13*y^2*z^3+5*y*z^4;541 poly f3 = 8*x^3+(a+12)*y^3+x*z^2+3;542 poly f4 = (a+7)*x^2*y^4+y^3*z^3+18*y^3*z^2;543 ideal I1 = f1,f2,f3,f4;544 545 // 2546 ring R = (0,a), (x,y,z), dp;547 minpoly = (a^5+a^2+2);548 poly f1 = 2*x*y^4*z^2+(a1)*x^2*y^3*z549 +(2*a)*x*y*z^2+7*y^3+(7*a+1);550 poly f2 = 2*x^2*y^4*z+(a)*x^2*y*z^2x*y^2*z^2551 +(2*a+3)*x^2*y*z12*x+(12*a)*y;552 poly f3 = (2*a)*y^5*z+x^2*y^2*zx*y^3*z553 +(a)*x*y^3+y^4+2*y^2*z;554 poly f4 = (3*a)*x*y^4*z^3+(a+1)*x^2*y^2*z555 x*y^3*z+4*y^3*z^2+(3*a)*x*y*z^3556 +4*z^2x+(a)*y;557 ideal I2 = f1,f2,f3,f4;558 559 // 3a560 ring R = (0,a), (v,w,x,y,z), dp;561 minpoly = (a^77*a+3);562 poly f1 = (a)*v+(a1)*w+x+(a+2)*y+z;563 poly f2 = v*w+(a1)*w*x+(a+2)*v*y+x*y+(a)*y*z;564 poly f3 = (a)*v*w*x+(a+5)*w*x*y+(a)*v*w*z565 +(a+2)*v*y*z+(a)*x*y*z;566 poly f4 = (a11)*v*w*x*y+(a+5)*v*w*x*z567 +(a)*v*w*y*z+(a)*v*x*y*z568 +(a)*w*x*y*z;569 poly f5 = (a+3)*v*w*x*y*z+(a+23);570 ideal I3a = f1,f2,f3,f4,f5;571 572 // 3b573 ring R = (0,a), (u,v,w,x,y,z), dp;574 minpoly = (a^77*a+3);575 poly f1 = (a)*u+(a+2)*v+w+x+y+z;576 poly f2 = u*v+v*w+w*x+x*y+(a+3)*u*z+y*z;577 poly f3 = u*v*w+v*w*x+(a+1)*w*x*y+u*v*z+u*y*z578 +x*y*z;579 poly f4 = (a1)*u*v*w*x+v*w*x*y+u*v*w*z580 +u*v*y*z+u*x*y*z+w*x*y*z;581 poly f5 = u*v*w*x*y+(a+1)*u*v*w*x*z+u*v*w*y*z582 +u*v*x*y*z+u*w*x*y*z+v*w*x*y*z;583 poly f6 = u*v*w*x*y*z+(a+2);584 ideal I3b = f1,f2,f3,f4,f5,f6;585 586 // 4587 ring R = (0,a), (w,x,y,z), dp;588 minpoly = (a^6+a^5+a^4+a^3+a^2+a+1);589 poly f1 = (a+5)*w^3*x^2*y+(a3)*w^2*x^3*y590 +(a+7)*w*x^2*y^2;591 poly f2 = (a)*w^5+(a+3)*w*x^2*y^2592 +(a^2+11)*x^2*y^2*z;593 poly f3 = (a+7)*w^3+12*x^3+4*w*x*y+(a)*z^3;594 poly f4 = 3*w^3+(a4)*x^3+x*y^2;595 ideal I4 = f1,f2,f3,f4;596 597 // 5598 ring R = (0,a), (w,x,y,z), dp;599 minpoly = (a^125*a^11+24*a^10115*a^9+551*a^8600 2640*a^7+12649*a^62640*a^5+551*a^4601 115*a^3+24*a^25*a+1);602 poly f1 = (2*a+3)*w*x^4*y^2+(a+1)*w^2*x^3*y*z603 +2*w*x*y^2*z^3+(7*a1)*x^3*z^4;604 poly f2 = 2*w^2*x^4*y+w^2*x*y^2*z^2605 +(a)*w*x^2*y^2*z^2606 +(a+11)*w^2*x*y*z^312*w*z^6607 +12*x*z^6;608 poly f3 = 2*x^5*y+w^2*x^2*y*zw*x^3*y*z609 w*x^3*z^2+(a)*x^4*z^2+2*x^2*y*z^3;610 poly f4 = 3*w*x^4*y^3+w^2*x^2*y*z^3611 w*x^3*y*z^3+(a+4)*x^3*y^2*z^3612 +3*w*x*y^3*z^3+(4*a)*y^2*z^6w*z^7613 +x*z^7;614 ideal I5 = f1,f2,f3,f4;615 616 // 6617 ring R = (0,a), (u,v,w,x,y,z), dp;618 minpoly = (a^2+5*a+1);619 poly f1 = u+v+w+x+y+z+(a);620 poly f2 = u*v+v*w+w*x+x*y+y*z+(a)*u+(a)*z;621 poly f3 = u*v*w+v*w*x+w*x*y+x*y*z+(a)*u*v622 +(a)*u*z+(a)*y*z;623 poly f4 = u*v*w*x+v*w*x*y+w*x*y*z+(a)*u*v*w624 +(a)*u*v*z+(a)*u*y*z+(a)*x*y*z;625 poly f5 = u*v*w*x*y+v*w*x*y*z+(a)*u*v*w*x626 +(a)*u*v*w*z+(a)*u*v*y*z+(a)*u*x*y*z627 +(a)*w*x*y*z;628 poly f6 = u*v*w*x*y*z+(a)*u*v*w*x*y629 +(a)*u*v*w*x*z+(a)*u*v*w*y*z630 +(a)*u*v*x*y*z+(a)*u*w*x*y*z631 +(a)*v*w*x*y*z;632 poly f7 = (a)*u*v*w*x*y*z1;633 ideal I6 = f1,f2,f3,f4,f5,f6,f7;634 635 // 7636 ring R = (0,a), (w,x,y,z), dp;637 minpoly = (a^816*a^7+19*a^6a^55*a^4+13*a^3638 9*a^2+13*a+17);639 poly f1 = (a^21)*x^2*y+2*w*x*z2*w640 +(a^2+1)*y;641 poly f2 = (a^3a3)*w^3*y+4*w*x^2*y+4*w^2*x*z642 +2*x^3*z+(a)*w^210*x^2+4*w*y10*x*z643 +(2*a^2+a);644 poly f3 = (a^2+a+11)*x*y*z+w*z^2w2*y;645 poly f4 = w*y^3+4*x*y^2*z+4*w*y*z^2+2*x*z^3646 +(2*a^3+a^2)*w*y+4*y^210*x*z10*z^2647 +(3*a^2+5);648 ideal I7 = f1,f2,f3,f4;649 650 // 8651 ring R = (0,a), (t,u,v,w,x,y,z), dp;652 minpoly = (a^7+10*a^5+5*a^3+10*a+1);653 poly f1 = v*x+w*yx*zwy;654 poly f2 = v*wu*x+x*yw*z+v+x+z;655 poly f3 = t*ww^2+x^2t;656 poly f4 = (a)*v^2u*y+y^2v*zz^2+u;657 poly f5 = t*v+v*w+(a^2a5)*x*yt*z+w*z+v+x+z658 +(a+1);659 poly f6 = t*u+u*w+(a11)*v*xt*y+w*yx*ztu660 +w+y;661 poly f7 = w^2*y^3w*x*y^3+x^2*y^3+w^2*y^2*z662 w*x*y^2*z+x^2*y^2*z+w^2*y*z^2663 w*x*y*z^2+x^2*y*z^2+w^2*z^3w*x*z^3664 +x^2*z^3;665 poly f8 = t^2*u^3+t^2*u^2*v+t^2*u*v^2+t^2*v^3666 t*u^3*xt*u^2*v*xt*u*v^2*xt*v^3*x667 +u^3*x^2+u^2*v*x^2+u*v^2*x^2668 +v^3*x^2;669 ideal I8 = f1,f2,f3,f4,f5,f6,f7,f8;670 */ 
Singular/singularlibs
rb1cef7 r204cf0 47 47 gitfan.lib gradedModules.lib \ 48 48 locnormal.lib modnormal.lib modular.lib \ 49 JMBTest.lib JMSConst.lib multigrading.lib \49 JMBTest.lib JMSConst.lib multigrading.lib nfmodsyz.lib \ 50 50 numerAlg.lib numerDecom.lib \ 51 51 orbitparam.lib parallel.lib polymake.lib\
Note: See TracChangeset
for help on using the changeset viewer.