Changeset 53e03a6 in git
- Timestamp:
- Dec 11, 2006, 7:51:12 PM (16 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a657104b677b4c461d018cbf3204d72d34ad66a9')
- Children:
- 688052c6836073bcb5c3e359253c9c06b18d72c7
- Parents:
- 0c0b9f11e094394b7ba5d064cb18da114c89975d
- File:
-
- 1 edited
-
Singular/LIB/AtkinsTest.lib (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/AtkinsTest.lib
r0c0b9f1 r53e03a6 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: AtkinsTest.lib,v 1. 2 2006-12-11 18:08:15Singular Exp $";2 version="$Id: AtkinsTest.lib,v 1.3 2006-12-11 18:51:12 Singular Exp $"; 3 3 category="Teaching"; 4 4 info=" … … 40 40 41 41 proc new(list L, number D) 42 //"USAGE: new(L,D);43 //RETURN: 1, if D does not already exist in L,44 //-1, if D does already exist in L45 //EXAMPLE:example new; shows an example46 //"47 42 "USAGE: new(L,D); 43 RETURN: 1, if D does not already exist in L, 44 -1, if D does already exist in L 45 EXAMPLE:example new; shows an example 46 " 47 { 48 48 number a=1; // a=1 bedeutet: D noch nicht in L vorhanden 49 49 int i; … … 58 58 59 59 return(a); 60 } 61 62 example 63 { "EXAMPLE:"; echo = 2; 60 } 61 example 62 { "EXAMPLE:"; echo = 2; 64 63 ring r = 0,x,dp; 65 64 list L=8976,-223456,556,-778,3,-55603,45,766677; 66 65 number D=-55603; 67 66 new(L,D); 68 67 } 69 68 70 69 71 70 72 71 proc bubblesort(list L) 73 //"USAGE: bubblesort(L);74 //RETURN: list L, sort in decreasing order75 //EXAMPLE:example bubblesort; shows an example76 //"77 72 "USAGE: bubblesort(L); 73 RETURN: list L, sort in decreasing order 74 EXAMPLE:example bubblesort; shows an example 75 " 76 { 78 77 number b; 79 78 int n,i,j; … … 95 94 96 95 return(L); 97 } 98 99 example 100 { "EXAMPLE:"; echo = 2; 96 } 97 example 98 { "EXAMPLE:"; echo = 2; 101 99 ring r = 0,x,dp; 102 100 list L=-567,-233,446,12,-34,8907; 103 101 bubblesort(L); 104 102 } 105 103 106 104 107 105 108 106 proc disc(number N, int k) 109 //"USAGE: disc(N,k);110 //RETURN: list L of negative discriminants D, sort in decreasing order111 //ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4N112 //NOTE: D=b^2-4*a, where 0<=b<=k and intPart((b^2)/4)+1<=a<=k for each b113 //EXAMPLE:example disc; shows an example114 //"115 107 "USAGE: disc(N,k); 108 RETURN: list L of negative discriminants D, sort in decreasing order 109 ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4N 110 NOTE: D=b^2-4*a, where 0<=b<=k and intPart((b^2)/4)+1<=a<=k for each b 111 EXAMPLE:example disc; shows an example 112 " 113 { 116 114 list L=-3,-4,-7; 117 115 number D; … … 133 131 L=bubblesort(L); 134 132 return(L); 135 } 136 137 example 138 { "EXAMPLE:"; echo = 2; 133 } 134 example 135 { "EXAMPLE:"; echo = 2; 139 136 ring R = 0,x,dp; 140 137 disc(2003,50); 141 138 } 142 139 143 140 144 141 145 142 proc Cornacchia(number d, number p) 146 //"USAGE: Cornacchia(d,p);147 //RETURN: x,y such that x^2+d*y^2=p with p prime,148 //-1, if the Diophantine equation has no solution,149 //0, if the parameters are wrong selected150 //ASSUME: 0<d<p151 //EXAMPLE:example Cornacchia; shows an example152 //"153 143 "USAGE: Cornacchia(d,p); 144 RETURN: x,y such that x^2+d*y^2=p with p prime, 145 -1, if the Diophantine equation has no solution, 146 0, if the parameters are wrong selected 147 ASSUME: 0<d<p 148 EXAMPLE:example Cornacchia; shows an example 149 " 150 { 154 151 155 152 if((d<0)||(p<d)) // (0)[Test if assumptions well-defined] … … 222 219 } 223 220 } 224 } 225 226 example 227 { "EXAMPLE:"; echo = 2; 221 } 222 example 223 { "EXAMPLE:"; echo = 2; 228 224 ring R = 0,x,dp; 229 225 Cornacchia(55,9551); 230 226 } 231 227 232 228 233 229 234 230 proc CornacchiaModified(number D, number p) 235 //"USAGE: CornacchiaModified(D,p);236 //RETURN: x,y such that x^2+|D|*y^2=p with p prime,237 //-1, if the Diophantine equation has no solution,238 //0, if the parameters are wrong selected239 //ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4p240 //EXAMPLE:example CornacchiaModified; shows an example241 //"242 231 "USAGE: CornacchiaModified(D,p); 232 RETURN: x,y such that x^2+|D|*y^2=p with p prime, 233 -1, if the Diophantine equation has no solution, 234 0, if the parameters are wrong selected 235 ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4p 236 EXAMPLE:example CornacchiaModified; shows an example 237 " 238 { 243 239 244 240 if((D>=0)||((D mod 4)==2)||((D mod 4)==3)||(absValue(D)>=4*p)) // (0)[Test if assumptions well-defined] … … 332 328 } 333 329 } 334 } 335 336 example 337 { "EXAMPLE:"; echo = 2; 330 } 331 example 332 { "EXAMPLE:"; echo = 2; 338 333 ring R = 0,x,dp; 339 334 CornacchiaModified(-107,1319); 340 335 } 341 336 342 337 343 338 344 339 proc pFactor1(number n,int B, list P) 345 //"USAGE: pFactor1(n,B,P); n to be factorized, B a bound , P a list of primes346 //RETURN: a list of factors of n or the message: no factor found347 //NOTE: Pollard's p-factorization348 //creates the product k of powers of primes (bounded by B) from349 //the list P with the idea that for a prime divisor p of n p-1|k350 //then p devides gcd(a^k-1,n) for some random a351 //EXAMPLE:example pFactor1; shows an example352 //"353 340 "USAGE: pFactor1(n,B,P); n to be factorized, B a bound , P a list of primes 341 RETURN: a list of factors of n or the message: no factor found 342 NOTE: Pollard's p-factorization 343 creates the product k of powers of primes (bounded by B) from 344 the list P with the idea that for a prime divisor p of n p-1|k 345 then p devides gcd(a^k-1,n) for some random a 346 EXAMPLE:example pFactor1; shows an example 347 " 348 { 354 349 int i; 355 350 number k=1; … … 370 365 if((d>1)&&(d<n)){return(d);} 371 366 return(n); 372 } 373 367 } 374 368 example 375 369 { "EXAMPLE:"; echo = 2; … … 385 379 386 380 proc maximum(list L) 387 //"USAGE: maximum(list L);388 //RETURN: the maximal number contained in list L389 //EXAMPLE:example maximum; shows an example390 //"391 381 "USAGE: maximum(list L); 382 RETURN: the maximal number contained in list L 383 EXAMPLE:example maximum; shows an example 384 " 385 { 392 386 number max=L[1]; 393 387 … … 402 396 403 397 return(max); 404 } 405 406 example 407 { "EXAMPLE:"; echo = 2; 398 } 399 example 400 { "EXAMPLE:"; echo = 2; 408 401 ring r = 0,x,dp; 409 402 list L=465,867,1233,4567,776544,233445,2334,556; 410 403 maximum(L); 411 404 } 412 405 413 406 414 407 415 408 proc cmod(number x, number y) 416 //"USAGE: cmod(x,y);417 //RETURN: x mod y418 //ASSUME: x,y out of Z and x,y<=2147483647419 //NOTE: this algorithm is a helping procedure to be able to calculate420 421 //EXAMPLE:example cmod; shows an example422 //"423 409 "USAGE: cmod(x,y); 410 RETURN: x mod y 411 ASSUME: x,y out of Z and x,y<=2147483647 412 NOTE: this algorithm is a helping procedure to be able to calculate 413 x mod y with x,y out of Z while working in the complex field 414 EXAMPLE:example cmod; shows an example 415 " 416 { 424 417 int rest=int(x-y*int(x/y)); 425 418 if(rest<0) … … 429 422 430 423 return(rest); 431 } 432 433 example 434 { "EXAMPLE:"; echo = 2; 424 } 425 example 426 { "EXAMPLE:"; echo = 2; 435 427 ring r = (complex,30,i),x,dp; 436 428 number x=-1004456; 437 429 number y=1233; 438 430 cmod(x,y); 439 431 } 440 432 441 433 442 434 443 435 proc sqr(number w, int k) 444 // "USAGE: sqr(w,k); 445 // RETURN: the square root of w 446 // ASSUME: w>=0 447 // NOTE: k describes the number of decimals being calculated in the real numbers, 448 // k, intPart(k/5) are inputs for the procedure "nt_solve" 449 // EXAMPLE:example sqr; shows an example 450 // " 451 { 452 poly f=var(1)^2-w; 453 def S=basering; 454 ring R=(real,k),var(1),dp; 455 poly f=imap(S,f); 456 ideal I=nt_solve(f,1.1,list(k,int(intPart(k/5)))); 457 number c=leadcoef(I[1]); 458 setring S; 459 number c=imap(R,c); 460 return(c); 461 } 462 463 example 464 { "EXAMPLE:"; echo = 2; 436 "USAGE: sqr(w,k); 437 RETURN: the square root of w 438 ASSUME: w>=0 439 NOTE: k describes the number of decimals being calculated in the real numbers, 440 k, intPart(k/5) are inputs for the procedure "nt_solve" 441 EXAMPLE:example sqr; shows an example 442 " 443 { 444 poly f=var(1)^2-w; 445 def S=basering; 446 ring R=(real,k),var(1),dp; 447 poly f=imap(S,f); 448 ideal I=nt_solve(f,1.1,list(k,int(intPart(k/5)))); 449 number c=leadcoef(I[1]); 450 setring S; 451 number c=imap(R,c); 452 return(c); 453 } 454 example 455 { "EXAMPLE:"; echo = 2; 465 456 ring R = (real,60),x,dp; 466 457 number ww=288469650108669535726081; 467 458 sqr(ww,60); 468 459 } 469 460 470 461 471 462 472 463 proc e(number z, int k) 473 // "USAGE: e(z,k); 474 // RETURN: e^z to the order k 475 // NOTE: k describes the number of summands being calculated in the exponential power series 476 // EXAMPLE:example e; shows an example 477 // " 478 { 479 number q=1; 480 number e=1; 481 482 int n; 483 for(n=1;n<=k;n++) 484 { 485 q=q*z/n; 486 e=e+q; 487 } 488 489 return(e); 490 } 491 492 example 493 { "EXAMPLE:"; echo = 2; 464 "USAGE: e(z,k); 465 RETURN: e^z to the order k 466 NOTE: k describes the number of summands being calculated in the exponential power series 467 EXAMPLE:example e; shows an example 468 " 469 { 470 number q=1; 471 number e=1; 472 473 int n; 474 for(n=1;n<=k;n++) 475 { 476 q=q*z/n; 477 e=e+q; 478 } 479 return(e); 480 } 481 482 example 483 { "EXAMPLE:"; echo = 2; 494 484 ring r = (real,30),x,dp; 495 485 number z=40.35; 496 486 e(z,1000); 497 487 } 498 488 499 489 500 490 501 491 proc jot(number t, int k) 502 //"USAGE: jot(t,k);503 //RETURN: the j-invariant of t504 //ASSUME: t is a complex number with positive imaginary part505 //NOTE: k describes the number of summands being calculated in the power series,506 //10*k is input for the procedure "e"507 //EXAMPLE:example jot; shows an example508 //"509 492 "USAGE: jot(t,k); 493 RETURN: the j-invariant of t 494 ASSUME: t is a complex number with positive imaginary part 495 NOTE: k describes the number of summands being calculated in the power series, 496 10*k is input for the procedure "e" 497 EXAMPLE:example jot; shows an example 498 " 499 { 510 500 number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j; 511 501 … … 538 528 j=(256*f+1)^3/f; 539 529 return(j); 540 541 542 example 543 530 } 531 532 example 533 { "EXAMPLE:"; echo = 2; 544 534 ring r = (complex,30,i),x,dp; 545 535 number t=(-7+i*sqr(7,250))/2; 546 536 jot(t,50); 547 537 } 548 538 549 539 550 540 551 541 proc round(number r) 552 // "USAGE: round(r); 553 // RETURN: the nearest number to r out of Z 554 // ASSUME: r should be a rational or a real number 555 // EXAMPLE:example round; shows an example 556 // " 542 "USAGE: round(r); 543 RETURN: the nearest number to r out of Z 544 ASSUME: r should be a rational or a real number 545 EXAMPLE:example round; shows an example 546 " 547 { 548 number a=absValue(r); 549 number v=r/a; 550 551 number d=10; 552 int e; 553 while(1) 554 { 555 e=e+1; 556 if(a-d^e<0) 557 557 { 558 number a=absValue(r); 559 number v=r/a; 560 561 number d=10; 562 int e; 563 while(1) 564 { 565 e=e+1; 566 if(a-d^e<0) 567 { 568 e=e-1; 569 break; 570 } 571 } 572 573 number b=a; 574 int k; 575 for(k=0;k<=e;k++) 576 { 577 while(1) 578 { 579 b=b-d^(e-k); 580 if(b<0) 581 { 582 b=b+d^(e-k); 583 break; 584 } 585 } 586 } 587 588 if(b<1/2) 589 { 590 return(v*(a-b)); 591 } 592 593 else 594 { 595 return(v*(a+1-b)); 596 } 558 e=e-1; 559 break; 597 560 } 598 599 example 600 { "EXAMPLE:"; echo = 2; 561 } 562 563 number b=a; 564 int k; 565 for(k=0;k<=e;k++) 566 { 567 while(1) 568 { 569 b=b-d^(e-k); 570 if(b<0) 571 { 572 b=b+d^(e-k); 573 break; 574 } 575 } 576 } 577 578 if(b<1/2) 579 { 580 return(v*(a-b)); 581 } 582 else 583 { 584 return(v*(a+1-b)); 585 } 586 } 587 example 588 { "EXAMPLE:"; echo = 2; 601 589 ring R = (real,50),x,dp; 602 590 number r=7357683445788723456321.6788643224; 603 591 round(r); 604 592 } 605 593 606 594 607 595 608 596 proc HilbertClassPolynomial(number D, int k) 609 //"USAGE: HilbertClassPolynomial(D,k);610 //RETURN: the monic polynomial of degree h(D) in Z[X] of which jot((D+sqr(D))/2) is a root611 //ASSUME: D is a negative discriminant612 //NOTE: k is input for the procedure "jot",613 //5*k is input for the procedure "sqr",614 //10*k describes the number of decimals being calculated in the complex numbers615 //EXAMPLE:example HilbertClassPolynomial; shows an example616 //"617 597 "USAGE: HilbertClassPolynomial(D,k); 598 RETURN: the monic polynomial of degree h(D) in Z[X] of which jot((D+sqr(D))/2) is a root 599 ASSUME: D is a negative discriminant 600 NOTE: k is input for the procedure "jot", 601 5*k is input for the procedure "sqr", 602 10*k describes the number of decimals being calculated in the complex numbers 603 EXAMPLE:example HilbertClassPolynomial; shows an example 604 " 605 { 618 606 if(D>=0) // (0)[Test if assumptions well-defined] 619 607 { … … 748 736 return(Q); 749 737 } 750 } 751 752 example 753 { "EXAMPLE:"; echo = 2; 738 } 739 example 740 { "EXAMPLE:"; echo = 2; 754 741 ring r = 0,x,dp; 755 742 number D=-23; 756 743 HilbertClassPolynomial(D,50); 757 744 } 758 745 759 746 760 747 761 748 proc RootsModp(int p, poly P) 762 //"USAGE: RootsModp(p,P);763 //RETURN: list of roots of the polynomial P modulo p with p prime764 //ASSUME: p>=3765 //NOTE: this algorithm will be called recursively, and it is understood766 //that all the operations are done in Z/pZ (excepting sqareRoot(d,p))767 //EXAMPLE:example RootsModp; shows an example768 //"769 749 "USAGE: RootsModp(p,P); 750 RETURN: list of roots of the polynomial P modulo p with p prime 751 ASSUME: p>=3 752 NOTE: this algorithm will be called recursively, and it is understood 753 that all the operations are done in Z/pZ (excepting sqareRoot(d,p)) 754 EXAMPLE:example RootsModp; shows an example 755 " 756 { 770 757 if(p<3) // (0)[Test if assumptions well-defined] 771 758 { … … 837 824 return(l); 838 825 } 839 } 840 841 example 842 { "EXAMPLE:"; echo = 2; 826 } 827 example 828 { "EXAMPLE:"; echo = 2; 843 829 ring r = 0,x,dp; 844 830 poly f=x4+2x3-5x2+x; … … 846 832 poly g=x5+112x4+655x3+551x2+1129x+831; 847 833 RootsModp(1223,g); 848 834 } 849 835 850 836 851 837 852 838 proc w(number D) 853 // "USAGE: w(D); 854 // RETURN: the number of roots of unity in the quadratic order of discriminant D 855 // ASSUME: D<0 a discriminant kongruent to 0 or 1 modulo 4 856 // EXAMPLE:example w; shows an example 857 // " 858 { 859 if((D>=0)||((D mod 4)==2)||((D mod 4)==3)) 860 { 861 ERROR("Parameter wrong selected!"); 862 } 863 864 else 865 { 866 if(D<-4) {return(2);} 867 if(D==-4){return(4);} 868 if(D==-3){return(6);} 869 } 870 } 871 872 example 873 { "EXAMPLE:"; echo = 2; 839 "USAGE: w(D); 840 RETURN: the number of roots of unity in the quadratic order of discriminant D 841 ASSUME: D<0 a discriminant kongruent to 0 or 1 modulo 4 842 EXAMPLE:example w; shows an example 843 " 844 { 845 if((D>=0)||((D mod 4)==2)||((D mod 4)==3)) 846 { 847 ERROR("Parameter wrong selected!"); 848 } 849 else 850 { 851 if(D<-4) {return(2);} 852 if(D==-4){return(4);} 853 if(D==-3){return(6);} 854 } 855 } 856 example 857 { "EXAMPLE:"; echo = 2; 874 858 ring r = 0,x,dp; 875 859 number D=-3; 876 860 w(D); 877 861 } 878 862 879 863 880 864 881 865 proc Atkin(number N, int K, int B) 882 //"USAGE: Atkin(N,K,B);883 //RETURN: 1, if N is prime,884 //-1, if N is not prime,885 //0, if the algorithm is not applicable, since there are too little discriminants886 //ASSUME: N is coprime to 6 and different from 1887 //NOTE: - K/2 is input for the procedure "disc",888 //K is input for the procedure "HilbertClassPolynomial",889 //B describes the number of recursions being calculated890 //- The basis of the the algorithm is the following theorem:891 //Let N be an integer coprime to 6 and different from 1 and E be an ellipic curve modulo N.892 //Assume that we know an integer m and a point P of E(Z/NZ) satisfying the following conditions.893 //(1) There exists a prime divisor q of m such that q>(4-th root(N)+1)^2.894 //(2) m*P=O(E)=(0:1:0).895 //(3) (m/q)*P=(x:y:t) with t element of (Z/NZ)*.896 //Then N is prime.897 //EXAMPLE:example Atkin; shows an example898 //"899 866 "USAGE: Atkin(N,K,B); 867 RETURN: 1, if N is prime, 868 -1, if N is not prime, 869 0, if the algorithm is not applicable, since there are too little discriminants 870 ASSUME: N is coprime to 6 and different from 1 871 NOTE: - K/2 is input for the procedure "disc", 872 K is input for the procedure "HilbertClassPolynomial", 873 B describes the number of recursions being calculated 874 - The basis of the the algorithm is the following theorem: 875 Let N be an integer coprime to 6 and different from 1 and E be an ellipic curve modulo N. 876 Assume that we know an integer m and a point P of E(Z/NZ) satisfying the following conditions. 877 (1) There exists a prime divisor q of m such that q>(4-th root(N)+1)^2. 878 (2) m*P=O(E)=(0:1:0). 879 (3) (m/q)*P=(x:y:t) with t element of (Z/NZ)*. 880 Then N is prime. 881 EXAMPLE:example Atkin; shows an example 882 " 883 { 900 884 if(N==1) {return(-1);} 901 885 if((N==2)||(N==3)) {return(1);} … … 1314 1298 } 1315 1299 } 1316 } 1317 1318 example 1319 { "EXAMPLE:"; echo = 2; 1300 } 1301 example 1302 { "EXAMPLE:"; echo = 2; 1320 1303 ring R = 0,x,dp; 1321 1304 printlevel=1; … … 1324 1307 Atkin(100019,100,5); 1325 1308 Atkin(10000079,100,2); 1326 1309 }
Note: See TracChangeset
for help on using the changeset viewer.
