Changeset 26508d in git for Singular/LIB/atkins.lib
- Timestamp:
- Jul 20, 2007, 1:26:46 PM (16 years ago)
- Branches:
- (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
- Children:
- 7c38bcb5f2878e308ca07340e00bfe7b4fdd406e
- Parents:
- 3eadab96fecbb57aa75151e22b15ce9249b408b6
- File:
-
- 1 edited
-
Singular/LIB/atkins.lib (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/atkins.lib
r3eadab r26508d 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: atkins.lib,v 1. 3 2007-07-20 10:02:38Singular Exp $";2 version="$Id: atkins.lib,v 1.4 2007-07-20 11:26:46 Singular Exp $"; 3 3 category="Teaching"; 4 4 info=" … … 95 95 bubblesort(L); 96 96 } 97 98 99 97 100 98 proc disc(number N, int k) … … 193 191 } 194 192 195 196 197 193 proc CornacchiaModified(number D, number p) 198 194 "USAGE: CornacchiaModified(D,p); … … 204 200 " 205 201 { 206 if((D>=0)||((D mod 4)==2)||((D mod 4)==3)||(absValue(D)>=4*p)) 207 {// (0)[Test if assumptions well-defined] 208 return(0); 209 // ERROR("Parameters wrong selected!"); 210 } 211 else 212 { 213 if(p==2) // (1)[Case p=2] 214 { 215 if((D+8)==intRoot(D+8)^2) 216 { 217 return(intRoot(D+8),1); 218 } 219 else 220 { 221 return(-1); 222 // ERROR("The Diophantine equation has no solution!"); 223 } 202 if((D>=0)||((D mod 4)==2)||((D mod 4)==3)||(absValue(D)>=4*p)) 203 {// (0)[Test if assumptions well-defined] 204 return(0); 205 // ERROR("Parameters wrong selected!"); 206 } 207 else 208 { 209 if(p==2) // (1)[Case p=2] 210 { 211 if((D+8)==intRoot(D+8)^2) 212 { 213 return(intRoot(D+8),1); 224 214 } 225 215 else 226 216 { 227 number k,x(0),a,b,l,r,c; 228 k=Jacobi(D,p); // (2)[Test if residue] 229 if(k==-1) 230 { 231 return(-1); 232 // ERROR("The Diophantine equation has no solution!"); 233 } 234 else 235 { 236 x(0)=squareRoot(D,p); // (3)[Compute square root] 237 if((x(0) mod 2)!=(D mod 2)) 238 { 239 x(0)=p-x(0); 240 } 241 a=2*p; 242 b=x(0); 243 l=intRoot(4*p); 244 while(b>l) // (4)[Euclidean algorithm] 245 { 246 r=a mod b; 247 a=b; 248 b=r; 249 } 250 c=(4*p-b^2)/absValue(D);// (5)[Test solution] 251 if((((4*p-b^2) mod absValue(D))!=0)||(c!=intRoot(c)^2)) 252 { 253 return(-1); 254 // ERROR("The Diophantine equation has no solution!"); 255 } 256 else 257 { 258 list L=b,intRoot(c); 259 return(L); 260 } 261 } 262 } 263 } 217 return(-1); 218 // ERROR("The Diophantine equation has no solution!"); 219 } 220 } 221 else 222 { 223 number k,x(0),a,b,l,r,c; 224 k=Jacobi(D,p); // (2)[Test if residue] 225 if(k==-1) 226 { 227 return(-1); 228 // ERROR("The Diophantine equation has no solution!"); 229 } 230 else 231 { 232 x(0)=squareRoot(D,p); // (3)[Compute square root] 233 if((x(0) mod 2)!=(D mod 2)) 234 { 235 x(0)=p-x(0); 236 } 237 a=2*p; 238 b=x(0); 239 l=intRoot(4*p); 240 while(b>l) // (4)[Euclidean algorithm] 241 { 242 r=a mod b; 243 a=b; 244 b=r; 245 } 246 c=(4*p-b^2)/absValue(D);// (5)[Test solution] 247 if((((4*p-b^2) mod absValue(D))!=0)||(c!=intRoot(c)^2)) 248 { 249 return(-1); 250 // ERROR("The Diophantine equation has no solution!"); 251 } 252 else 253 { 254 list L=b,intRoot(c); 255 return(L); 256 } 257 } 258 } 259 } 264 260 } 265 261 example … … 268 264 CornacchiaModified(-107,1319); 269 265 } 270 271 272 266 273 267 proc pFactor1(number n,int B, list P) … … 281 275 " 282 276 { 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 277 int i; 278 number k=1; 279 number w; 280 while(i<size(P)) 281 { 282 i++; 283 w=P[i]; 284 if(w>B) {break;} 285 while(w*P[i]<=B) 286 { 287 w=w*P[i]; 288 } 289 k=k*w; 290 } 291 number a=random(2,2147483629); 292 number d=gcdN(powerN(a,k,n)-1,n); 293 if((d>1)&&(d<n)){return(d);} 294 return(n); 301 295 } 302 296 example … … 310 304 } 311 305 312 313 314 306 proc maximum(list L) 315 307 "USAGE: maximum(list L); … … 318 310 " 319 311 { 320 number max=L[1]; 321 322 int i; 323 for(i=2;i<=size(L);i++) 324 { 325 if(L[i]>max) 326 { 327 max=L[i]; 328 } 329 } 330 331 return(max); 312 number max=L[1]; 313 int i; 314 for(i=2;i<=size(L);i++) 315 { 316 if(L[i]>max) 317 { 318 max=L[i]; 319 } 320 } 321 return(max); 332 322 } 333 323 example … … 337 327 maximum(L); 338 328 } 339 340 341 329 342 330 static proc cmod(number x, number y) … … 349 337 " 350 338 { 351 int rest=int(x-y*int(x/y)); 352 if(rest<0) 353 { 354 rest=rest+int(y); 355 } 356 357 return(rest); 339 int rest=int(x-y*int(x/y)); 340 if(rest<0) 341 { 342 rest=rest+int(y); 343 } 344 return(rest); 358 345 } 359 346 example … … 365 352 } 366 353 367 368 369 static proc sqr(number w, int k) 354 proc sqr(number w, int k) 370 355 "USAGE: sqr(w,k); 371 356 RETURN: the square root of w … … 404 389 number q=1; 405 390 number e=1; 406 407 391 int n; 408 392 for(n=1;n<=k;n++) … … 432 416 " 433 417 { 434 number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j; 435 436 number pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989; 437 438 tr=repart(t); 439 ti=impart(t); 440 if(tr==-1/2){qr1=-1;} 441 if(tr==0){qr1=1;} 442 if((tr!=-1/2)&&(tr!=0)) 443 { 444 tr=tr-round(tr); 445 qr1=expo(2*i*pi*tr,10*k); 446 } 447 448 qi1=expo(-pi*ti,10*k); 449 q1=qr1*qi1^2; 450 q2=q1^2; 451 452 int n=1; 453 while(n<=k) 454 { 455 m1=m1+(-1)^n*(q1^(n*(3*n-1)/2)+q1^(n*(3*n+1)/2)); 456 m2=m2+(-1)^n*(q2^(n*(3*n-1)/2)+q2^(n*(3*n+1)/2)); 457 n=n+1; 458 } 459 460 f=q1*((1+m2)/(1+m1))^24; 461 462 j=(256*f+1)^3/f; 463 return(j); 464 } 465 418 number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j; 419 420 number pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989; 421 422 tr=repart(t); 423 ti=impart(t); 424 if(tr==-1/2){qr1=-1;} 425 if(tr==0){qr1=1;} 426 if((tr!=-1/2)&&(tr!=0)) 427 { 428 tr=tr-round(tr); 429 qr1=expo(2*i*pi*tr,10*k); 430 } 431 432 qi1=expo(-pi*ti,10*k); 433 q1=qr1*qi1^2; 434 q2=q1^2; 435 436 int n=1; 437 while(n<=k) 438 { 439 m1=m1+(-1)^n*(q1^(n*(3*n-1)/2)+q1^(n*(3*n+1)/2)); 440 m2=m2+(-1)^n*(q2^(n*(3*n-1)/2)+q2^(n*(3*n+1)/2)); 441 n++; 442 } 443 444 f=q1*((1+m2)/(1+m1))^24; 445 446 j=(256*f+1)^3/f; 447 return(j); 448 } 466 449 example 467 450 { "EXAMPLE:"; echo = 2; … … 470 453 jOft(t,50); 471 454 } 472 473 474 455 475 456 proc round(number r) … … 538 519 " 539 520 { 540 if(D>=0) // (0)[Test if assumptions well-defined] 541 { 542 ERROR("Parameter wrong selected!"); 543 } 544 545 else 546 { 547 def S=basering; 548 ring R=0,x,dp; 549 550 string s1,s2,s3; 551 number a1,b1,t1,g1; 552 number D=imap(S,D); 553 number B=intRoot(absValue(D)/3); 554 555 ring C=(complex,10*k,i),x,dp; 556 number D=imap(S,D); 557 558 poly P=1; // (1)[Initialize] 559 number b=cmod(D,2); 560 number B=imap(R,B); 561 562 number t,a,g,tau,j; 563 list L; 564 565 int step=2; 566 while(1) 567 { 568 if(step==2) // (2)[Initialize a] 569 { 570 t=(b^2-D)/4; 571 L=b,1; 572 a=maximum(L); 573 step=3; 574 } 575 576 if(step==3) // (3)[Test] 577 { 578 if((cmod(t,a)!=0)) 579 { 580 step=4; 581 } 582 583 else 584 { 585 s1=string(a); 586 s2=string(b); 587 s3=string(t); 588 589 setring R; 590 execute("a1="+s1+";"); 591 execute("b1="+s2+";"); 592 execute("t1="+s3+";"); 593 g1=gcd(gcd(a1,b1),t1/a1); 594 setring C; 595 g=imap(R,g1); 596 597 if(g!=1) 598 { 599 step=4; 600 } 601 602 else 603 { 604 tau=(-b+i*sqr(absValue(D),5*k))/(2*a); 605 j=jOft(tau,k); 606 if((a==b)||(a^2==t)||(b==0)) 607 { 608 P=P*(var(1)-repart(j)); 609 step=4; 610 } 611 612 else 613 { 614 P=P*(var(1)^2-2*repart(j)*var(1)+repart(j)^2+impart(j)^2); 615 step=4; 616 } 617 } 618 } 619 } 620 621 if(step==4) // (4)[Loop on a] 622 { 623 a=a+1; 624 if(a^2<=t) 625 { 626 step=3; 627 continue; 628 } 629 630 else 631 { 632 step=5; 633 } 634 } 635 636 if(step==5) // (5)[Loop on b] 637 { 638 b=b+2; 639 if(b<=B) 640 { 641 step=2; 642 } 643 644 else 645 { 646 break; 647 } 648 } 649 } 650 651 matrix M=coeffs(P,var(1)); 652 653 list liste; 654 int n; 655 for(n=1;n<=nrows(M);n++) 656 { 657 liste[n]=round(repart(number(M[n,1]))); 658 } 659 660 poly Q; 661 int m; 662 for(m=1;m<=size(liste);m++) 663 { 664 Q=Q+liste[m]*var(1)^(m-1); 665 } 666 667 string s=string(Q); 668 setring S; 669 execute("poly Q="+s+";"); 670 return(Q); 671 } 521 if(D>=0) // (0)[Test if assumptions well-defined] 522 { 523 ERROR("Parameter wrong selected!"); 524 } 525 else 526 { 527 def S=basering; 528 ring R=0,x,dp; 529 530 string s1,s2,s3; 531 number a1,b1,t1,g1; 532 number D=imap(S,D); 533 number B=intRoot(absValue(D)/3); 534 535 ring C=(complex,10*k,i),x,dp; 536 number D=imap(S,D); 537 538 poly P=1; // (1)[Initialize] 539 number b=cmod(D,2); 540 number B=imap(R,B); 541 542 number t,a,g,tau,j; 543 list L; 544 545 int step=2; 546 while(1) 547 { 548 if(step==2) // (2)[Initialize a] 549 { 550 t=(b^2-D)/4; 551 L=b,1; 552 a=maximum(L); 553 step=3; 554 } 555 556 if(step==3) // (3)[Test] 557 { 558 if((cmod(t,a)!=0)) 559 { 560 step=4; 561 } 562 else 563 { 564 s1=string(a); 565 s2=string(b); 566 s3=string(t); 567 568 setring R; 569 execute("a1="+s1+";"); 570 execute("b1="+s2+";"); 571 execute("t1="+s3+";"); 572 g1=gcd(gcd(a1,b1),t1/a1); 573 setring C; 574 g=imap(R,g1); 575 576 if(g!=1) 577 { 578 step=4; 579 } 580 else 581 { 582 tau=(-b+i*sqr(absValue(D),5*k))/(2*a); 583 j=jOft(tau,k); 584 if((a==b)||(a^2==t)||(b==0)) 585 { 586 P=P*(var(1)-repart(j)); 587 step=4; 588 } 589 else 590 { 591 P=P*(var(1)^2-2*repart(j)*var(1)+repart(j)^2+impart(j)^2); 592 step=4; 593 } 594 } 595 } 596 } 597 598 if(step==4) // (4)[Loop on a] 599 { 600 a=a+1; 601 if(a^2<=t) 602 { 603 step=3; 604 continue; 605 } 606 else 607 { 608 step=5; 609 } 610 } 611 612 if(step==5) // (5)[Loop on b] 613 { 614 b=b+2; 615 if(b<=B) 616 { 617 step=2; 618 } 619 else 620 { 621 break; 622 } 623 } 624 } 625 626 matrix M=coeffs(P,var(1)); 627 628 list liste; 629 int n; 630 for(n=1;n<=nrows(M);n++) 631 { 632 liste[n]=round(repart(number(M[n,1]))); 633 } 634 635 poly Q; 636 int m; 637 for(m=1;m<=size(liste);m++) 638 { 639 Q=Q+liste[m]*var(1)^(m-1); 640 } 641 642 string s=string(Q); 643 setring S; 644 execute("poly Q="+s+";"); 645 return(Q); 646 } 672 647 } 673 648 example … … 677 652 HilbertClassPoly(D,50); 678 653 } 679 680 681 654 682 655 proc rootsModp(int p, poly P) … … 689 662 " 690 663 { 691 if(p<3) // (0)[Test if assumptions well-defined] 692 { 693 ERROR("Parameter wrong selected, since p<3!"); 694 } 695 696 else 697 { 698 def S=basering; 699 ring R=p,var(1),dp; 700 701 poly P=imap(S,P); 702 number d; 703 int a; 704 list L; 705 706 poly A=gcd(var(1)^p-var(1),P); // (1)[Isolate roots in Z/pZ] 707 if(subst(A,var(1),0)==0) 708 { 709 L[1]=0; 710 A=A/var(1); 711 } 712 713 if(deg(A)==0) // (2)[Small degree?] 714 { 715 return(L); 716 } 717 718 if(deg(A)==1) 719 { 720 matrix M=coeffs(A,var(1)); 721 L[size(L)+1]=-leadcoef(M[1,1])/leadcoef(M[2,1]); 722 setring S; 723 list L=imap(R,L); 724 return(L); 725 } 726 727 if(deg(A)==2) 728 { 729 matrix M=coeffs(A,var(1)); 730 d=leadcoef(M[2,1])^2-4*leadcoef(M[1,1])*leadcoef(M[3,1]); 731 732 ring T=0,var(1),dp; 733 number d=imap(R,d); 734 number e=squareRoot(d,p); 735 setring R; 736 number e=imap(T,e); 737 738 L[size(L)+1]=(-leadcoef(M[2,1])+e)/(2*leadcoef(M[3,1])); 739 L[size(L)+1]=(-leadcoef(M[2,1])-e)/(2*leadcoef(M[3,1])); 740 setring S; 741 list L=imap(R,L); 742 return(L); 743 } 744 745 poly B=1; // (3)[Random splitting] 746 poly C; 747 while((deg(B)==0)||(deg(B)==deg(A))) 748 { 749 a=random(0,p-1); 750 B=gcd((var(1)+a)^((p-1)/2)-1,A); 751 C=A/B; 752 } 753 754 setring S; // (4)[Recurse] 755 poly B=imap(R,B); 756 poly C=imap(R,C); 757 list l=L+rootsModp(p,B)+rootsModp(p,C); 758 return(l); 759 } 664 if(p<3) // (0)[Test if assumptions well-defined] 665 { 666 ERROR("Parameter wrong selected, since p<3!"); 667 } 668 else 669 { 670 def S=basering; 671 ring R=p,var(1),dp; 672 673 poly P=imap(S,P); 674 number d; 675 int a; 676 list L; 677 678 poly A=gcd(var(1)^p-var(1),P); // (1)[Isolate roots in Z/pZ] 679 if(subst(A,var(1),0)==0) 680 { 681 L[1]=0; 682 A=A/var(1); 683 } 684 685 if(deg(A)==0) // (2)[Small degree?] 686 { 687 return(L); 688 } 689 690 if(deg(A)==1) 691 { 692 matrix M=coeffs(A,var(1)); 693 L[size(L)+1]=-leadcoef(M[1,1])/leadcoef(M[2,1]); 694 setring S; 695 list L=imap(R,L); 696 return(L); 697 } 698 699 if(deg(A)==2) 700 { 701 matrix M=coeffs(A,var(1)); 702 d=leadcoef(M[2,1])^2-4*leadcoef(M[1,1])*leadcoef(M[3,1]); 703 704 ring T=0,var(1),dp; 705 number d=imap(R,d); 706 number e=squareRoot(d,p); 707 setring R; 708 number e=imap(T,e); 709 710 L[size(L)+1]=(-leadcoef(M[2,1])+e)/(2*leadcoef(M[3,1])); 711 L[size(L)+1]=(-leadcoef(M[2,1])-e)/(2*leadcoef(M[3,1])); 712 setring S; 713 list L=imap(R,L); 714 return(L); 715 } 716 717 poly B=1; // (3)[Random splitting] 718 poly C; 719 while((deg(B)==0)||(deg(B)==deg(A))) 720 { 721 a=random(0,p-1); 722 B=gcd((var(1)+a)^((p-1)/2)-1,A); 723 C=A/B; 724 } 725 726 setring S; // (4)[Recurse] 727 poly B=imap(R,B); 728 poly C=imap(R,C); 729 list l=L+rootsModp(p,B)+rootsModp(p,C); 730 return(l); 731 } 760 732 } 761 733 example … … 767 739 rootsModp(1223,g); 768 740 } 769 770 771 741 772 742 proc wUnit(number D) … … 1225 1195 } 1226 1196 } 1227 1228 1197 example 1229 1198 { "EXAMPLE:"; echo = 2;
Note: See TracChangeset
for help on using the changeset viewer.
