# source:git/Singular/LIB/atkins.lib@a8a6b8

spielwiese
Last change on this file since a8a6b8 was a8a6b8, checked in by Hans Schoenemann <hannes@…>, 6 years ago
chg: version/date of libs, p1
• Property mode set to `100644`
File size: 32.4 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="version atkins.lib 4.1.1.0 Dec_2017 ";
3category="Teaching";
4info="
5LIBRARY:  atkins.lib     Procedures for teaching cryptography
6AUTHOR:                  Stefan Steidel, steidel@mathematik.uni-kl.de
7
8NOTE: The library contains auxiliary procedures to compute the elliptic
9      curve primality test of Atkin and the Atkin's Test itself.
10      The library is intended to be used for teaching purposes but not
11      for serious computations. Sufficiently high printlevel allows to
12      control each step, thus illustrating the algorithms at work.
13
14
15PROCEDURES:
16  newTest(L,D)              checks if number D already exists in list L
17  bubblesort(L)             sorts elements of the list L
18  disc(N,k)                 generates a list of negative discriminants
19  Cornacchia(d,p)           computes solution (x,y) for x^2+d*y^2=p
20  CornacchiaModified(D,p)   computes solution (x,y) for x^2+|D|*y^2=4p
21  maximum(L)                computes the maximal number contained in L
22  sqr(w,k)                  computes the square root of w w.r.t. accuracy k
23  expo(z,k)                 computes exp(z)
24  jOft(t,k)                 computes the j-invariant of t
25  round(r)                  rounds r to the nearest number out of Z
26  HilbertClassPoly(D,k)     computes the Hilbert Class Polynomial
27  rootsModp(p,P)            computes roots of the polynomial P modulo p
28  wUnit(D)                  computes the number of units in Q(sqr(D))
29  Atkin(N,K,B)              tries to prove that N is prime
30
31";
32
33LIB "crypto.lib";
34LIB "general.lib";
35LIB "ntsolve.lib";
36LIB "inout.lib";
37
38///////////////////////////////////////////////////////////////////////////////
39
40proc newTest(list L, def D)
41"USAGE: newTest(L,D);
42RETURN:  1, if D does not already exist in L,
43        -1, if D does already exist in L
44EXAMPLE:example new; shows an example
45"
46{
47  int a=1;                // a=1 means: D does not already exist in L
48  int i;
49  for(i=1;i<=size(L);i++)
50  {
51    if(D==L[i])
52    {
53      a=-1;                 // a=-1 means: D does already exist in L
54      break;
55    }
56  }
57 return(a);
58}
59example
60{ "EXAMPLE:"; echo = 2;
61    ring r = 0,x,dp;
62    list L=8976,-223456,556,-778,3,-55603,45,766677;
63    number D=-55603;
64    newTest(L,D);
65}
66
67
68proc bubblesort(list L)
69"USAGE: bubblesort(L);
70RETURN: list L, sort in decreasing order
71EXAMPLE:example bubblesort; shows an example
72"
73{
74  def b;
75  int n,i,j;
76  while(j==0)
77  {
78    i=i+1;
79    j=1;
80    for(n=1;n<=size(L)-i;n++)
81    {
82      if(L[n]<L[n+1])
83      {
84        b=L[n];
85        L[n]=L[n+1];
86        L[n+1]=b;
87        j=0;
88      }
89    }
90  }
91  return(L);
92}
93example
94{ "EXAMPLE:"; echo = 2;
95    ring r = 0,x,dp;
96    list L=-567,-233,446,12,-34,8907;
97    bubblesort(L);
98}
99
100
101proc disc(bigint N, int k)
102"USAGE: disc(N,k);
103RETURN: list L of negative discriminants D, sorted in decreasing order
104ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4N
105NOTE:   D=b^2-4*a, where 0<=b<=k and intPart((b^2)/4)+1<=a<=k for each b
106EXAMPLE:example disc; shows an example
107"
108{
109  list L=-3,-4,-7;
110  bigint D;
111  bigint B;
112  int a,b;
113  for(b=0;b<=k;b++)
114  {
115    B=b^2;
116    for(a=int(B div 4)+1;a<=k;a++)
117    {
118      D=-4*a+B;
119      if((D<0)&&((D mod 4)!=-2)&&((D mod 4)!=-1)&&(absValue(D)<4*N)&&(newTest(L,D)==1))
120      {
121        L[size(L)+1]=D;
122      }
123    }
124  }
125  L=bubblesort(L);
126  return(L);
127}
128example
129{ "EXAMPLE:"; echo = 2;
130    disc(2003,50);
131}
132
133proc Cornacchia(bigint d, bigint p)
134"USAGE: Cornacchia(d,p);
135RETURN: x,y such that x^2+d*y^2=p with p prime,
136        -1, if the Diophantine equation has no solution,
137         0, if the parameters are wrong selected
138ASSUME: 0<d<p
139EXAMPLE:example Cornacchia; shows an example
140"
141{
142  if((d<0)||(p<d))              // (0)[Test if assumptions well-defined]
143  {
144    return(0);
145    // ERROR("Parameters wrong selected! It has to be 0<d<p!");
146  }
147  else
148  {
149    bigint k,x(0),a,b,l,r,c,i;
150
151    k=Jacobi(-d,p);             // (1)[Test if residue]
152    if(k==-1)
153    {
154      return(-1);
155      // ERROR("The Diophantine equation has no solution!");
156    }
157    else
158    {
159       x(0)=squareRoot(-d,p);   // (2)[Compute square root]
160       if((p div 2>=x(0))||(p<=x(0)))
161       {
162           x(0)=-x(0) mod p+p;
163       }
164       a=p;
165       b=x(0);
166       l=intRoot(p);
167       while(b>l)               // (3)[Euclidean algorithm]
168       {
169          r=a mod b;
170          a=b;
171          b=r;
172       }
173       c=(p-b^2) div d;             // (4)[Test solution]
174       i=intRoot(c);
175       if((((p-b^2) mod d)!=0)||(c!=i^2))
176       {
177          return(-1);
178       }
179       else
180       {
181          list L=b,i;
182          return(L);
183       }
184    }
185  }
186}
187example
188{ "EXAMPLE:"; echo = 2;
189    Cornacchia(55,9551);
190}
191
192
193proc CornacchiaModified(bigint D, bigint p)
194"USAGE: CornacchiaModified(D,p);
195RETURN: x,y such that x^2+|D|*y^2=p with p prime,
196        -1, if the Diophantine equation has no solution,
197         0, if the parameters are wrong selected
198ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4p
199EXAMPLE:example CornacchiaModified; shows an example
200"
201{
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);
214      }
215      else
216      {
217        return(-1);
218        // ERROR("The Diophantine equation has no solution!");
219      }
220    }
221    else
222    {
223      bigint 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)) // D is <0
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) div absValue(D);     // (5)[Test solution]
247        bigint root_c=intRoot(c);
248        if((((4*p-b^2) mod absValue(D))!=0)||(c!=root_c^2))
249        {
250        "3";
251          return(-1);
252          // ERROR("The Diophantine equation has no solution!");
253        }
254        else
255        {
256          list L=b,root_c;
257          return(L);
258        }
259      }
260    }
261  }
262}
263example
264{ "EXAMPLE:"; echo = 2;
265    CornacchiaModified(-107,1319);
266}
267
268
269static proc pFactor1(number n,int B, list P)
270"USAGE: pFactor1(n,B,P); n to be factorized, B a bound , P a list of primes
271RETURN: a list of factors of n or the message: no factor found
272NOTE:   Pollard's p-factorization
273        creates the product k of powers of primes (bounded by B)  from
274        the list P with the idea that for a prime divisor p of n p-1|k
275        then p devides gcd(a^k-1,n) for some random a
276EXAMPLE:example pFactor1; shows an example
277"
278{
279  int i;
280  number k=1;
281  number w;
282  while(i<size(P))
283  {
284    i++;
285    w=P[i];
286    if(w>B) {break;}
287    while(w*P[i]<=B)
288    {
289      w=w*P[i];
290    }
291    k=k*w;
292  }
293  number a=random(2,2147483629);
294  number d=gcd(powerN(a,k,n)-1,n);
295  if((d>1)&&(d<n)){return(d);}
296  return(n);
297}
298example
299{ "EXAMPLE:"; echo = 2;
300   ring R = 0,z,dp;
301   list L=primList(1000);
302   pFactor1(1241143,13,L);
303   number h=10;
304   h=h^30+25;
305   pFactor1(h,20,L);
306}
307
308
309proc maximum(list L)
310"USAGE: maximum(list L);
311RETURN: the maximal number contained in list L
312EXAMPLE:example maximum; shows an example
313"
314{
315  number max=L[1];
316  int i;
317  for(i=2;i<=size(L);i++)
318  {
319    if(L[i]>max)
320    {
321      max=L[i];
322    }
323  }
324  return(max);
325}
326example
327{ "EXAMPLE:"; echo = 2;
328    ring r = 0,x,dp;
329    list L=465,867,1233,4567,776544,233445,2334,556;
330    maximum(L);
331}
332
333
334static proc cmod(number x, number y)
335"USAGE: cmod(x,y);
336RETURN: x mod y
337ASSUME: x,y out of Z and x,y<=2147483647
338NOTE:   this algorithm is a helping procedure to be able to calculate
339        x mod y with x,y out of Z while working in the complex field
340EXAMPLE:example cmod; shows an example
341"
342{
343  int rest=int(x-y*int(x/y));
344  if(rest<0)
345  {
346    rest=rest+int(y);
347  }
348  return(rest);
349}
350example
351{ "EXAMPLE:"; echo = 2;
352    ring r = (complex,30,i),x,dp;
353    number x=-1004456;
354    number y=1233;
355    cmod(x,y);
356}
357
358
359proc sqr(number w, int k)
360"USAGE: sqr(w,k);
361RETURN: the square root of w
362ASSUME: w>=0
363NOTE:   k describes the number of decimals being calculated in the real numbers,
364        k, intPart(k/5) are inputs for the procedure "nt_solve"
365EXAMPLE:example sqr; shows an example
366"
367{
368  poly f=var(1)^2-w;
369  def S=basering;
370  ring R=(real,k),var(1),dp;
371  poly f=imap(S,f);
372  ideal I=nt_solve(f,1.1,list(k,k div 5));
374  setring S;
375  number c=imap(R,c);
376  return(c);
377}
378example
379{ "EXAMPLE:"; echo = 2;
380    ring R = (real,60),x,dp;
381    number ww=288469650108669535726081;
382    sqr(ww,60);
383}
384
385
386proc expo(number z, int k)
387"USAGE: expo(z,k);
388RETURN: e^z to the order k
389NOTE:   k describes the number of summands being calculated in the exponential power series
390EXAMPLE:example expo; shows an example
391"
392{
393  number q=1;
394  number e=1;
395  int n;
396  for(n=1;n<=k;n++)
397  {
398    q=q*z/n;
399    e=e+q;
400  }
401  return(e);
402}
403example
404{ "EXAMPLE:"; echo = 2;
405    ring r = (real,30),x,dp;
406    number z=40.35;
407    expo(z,1000);
408}
409
410
411proc jOft(number t, int k)
412"USAGE: jOft(t,k);
413RETURN: the j-invariant of t
414ASSUME: t is a complex number with positive imaginary part
415NOTE:   k describes the number of summands being calculated in the power series,
416        10*k is input for the procedure @code{expo}
417EXAMPLE:example jOft; shows an example
418"
419{
420  number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j;
421
422  number pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989;
423
424  tr=repart(t);
425  ti=impart(t);
426  if(tr==-1/2){qr1=-1;}
427  else
428  {
429    if(tr==0){qr1=1;}
430    else
431    {
432      tr=tr-round(tr);
433      qr1=expo(2*i*pi*tr,10*k);
434    }
435  }
436
437  qi1=expo(-pi*ti,10*k);
438  q1=qr1*qi1^2;
439  q2=q1^2;
440
441  int n=1;
442  while(n<=k)
443  {
444    m1=m1+(-1)^n*(q1^(n*(3*n-1) div 2)+q1^(n*(3*n+1) div 2));
445    m2=m2+(-1)^n*(q2^(n*(3*n-1) div 2)+q2^(n*(3*n+1) div 2));
446    n++;
447  }
448
449  f=q1*((1+m2)/(1+m1))^24;
450
451  j=(256*f+1)^3/f;
452  return(j);
453}
454example
455{ "EXAMPLE:"; echo = 2;
456    ring r = (complex,30,i),x,dp;
457    number t=(-7+i*sqr(7,250))/2;
458    jOft(t,50);
459}
460
461
462proc round(number r)
463"USAGE: round(r);
464RETURN: the nearest number to r out of Z
465ASSUME: r should be a rational or a real number
466EXAMPLE:example round; shows an example
467"
468{
469  number a=absValue(r);
470  number v=r/a;
471
472  number d=10;
473  int e;
474  while(1)
475  {
476    e=e+1;
477    if(a-d^e<0)
478    {
479      e=e-1;
480      break;
481    }
482  }
483
484  number b=a;
485  int k;
486  for(k=0;k<=e;k++)
487  {
488    while(1)
489    {
490      b=b-d^(e-k);
491      if(b<0)
492      {
493        b=b+d^(e-k);
494        break;
495      }
496    }
497  }
498
499  if(b<1/2)
500  {
501    return(v*(a-b));
502  }
503  else
504  {
505    return(v*(a+1-b));
506  }
507}
508example
509{ "EXAMPLE:"; echo = 2;
510    ring R = (real,50),x,dp;
511    number r=7357683445788723456321.6788643224;
512    round(r);
513}
514
515
516proc HilbertClassPoly(bigint D, int k)
517"USAGE: HilbertClassPoly(D,k);
518RETURN: the monic polynomial of degree h(D) in Z[X] of which jOft((D+sqr(D))/2) is a root
519ASSUME: D is a negative discriminant
520NOTE:   k is input for the procedure "jOft",
521        5*k is input for the procedure "sqr",
522        10*k describes the number of decimals being calculated in the complex numbers
523EXAMPLE:example HilbertClassPoly; shows an example
524"
525{
526  if(D>=0)                     // (0)[Test if assumptions well-defined]
527  {
528    ERROR("Parameter wrong selected!");
529  }
530  else
531  {
532    def S=basering;
533
534    string s1,s2,s3;
535    bigint B=intRoot(absValue(D) div 3);
536
537    ring C=(complex,10*k,i),x,dp;
538    number DD=D;
539
540    poly P=1;                  // (1)[Initialize]
541    number b=cmod(DD,2);
542
543    number t,a,g,tau,j;
544    list L;
545
546    bigint a1,b1,t1,g1;
547    int step=2;
548    while(1)
549    {
550      if(step==2)              // (2)[Initialize a]
551      {
552        t=(b^2-DD)/4;
553        L=b,1;
554        a=maximum(L);
555        step=3;
556      }
557
558      if(step==3)              // (3)[Test]
559      {
560        if((cmod(t,a)!=0))
561        {
562          step=4;
563        }
564        else
565        {
566          s1=string(a);
567          s2=string(b);
568          s3=string(t);
569
570          execute("a1="+s1+";");
571          execute("b1="+s2+";");
572          execute("t1="+s3+";");
573          g1=gcd(gcd(a1,b1),t1 div a1);
574          setring C;
575          g=g1;
576
577          if(g!=1)
578          {
579            step=4;
580          }
581          else
582          {
583            tau=(-b+i*sqr(absValue(DD),5*k))/(2*a);
584            j=jOft(tau,k);
585            if((a==b)||(a^2==t)||(b==0))
586            {
587              P=P*(var(1)-repart(j));
588              step=4;
589            }
590            else
591            {
592              P=P*(var(1)^2-2*repart(j)*var(1)+repart(j)^2+impart(j)^2);
593              step=4;
594            }
595          }
596        }
597      }
598
599      if(step==4)              // (4)[Loop on a]
600      {
601        a=a+1;
602        if(a^2<=t)
603        {
604          step=3;
605          continue;
606        }
607        else
608        {
609          step=5;
610        }
611      }
612
613      if(step==5)              // (5)[Loop on b]
614      {
615        b=b+2;
616        if(b<=B)
617        {
618          step=2;
619        }
620        else
621        {
622          break;
623        }
624      }
625    }
626
627    matrix M=coeffs(P,var(1));
628
629    list liste;
630    int n;
631    for(n=1;n<=nrows(M);n++)
632    {
633      liste[n]=round(repart(number(M[n,1])));
634    }
635
636    poly Q;
637    int m;
638    for(m=1;m<=size(liste);m++)
639    {
640      Q=Q+liste[m]*var(1)^(m-1);
641    }
642
643    string s=string(Q);
644    setring S;
645    execute("poly Q="+s+";");
646    return(Q);
647  }
648}
649example
650{ "EXAMPLE:"; echo = 2;
651    ring r = 0,x,dp;
652    bigint D=-23;
653    HilbertClassPoly(D,50);
654}
655
656
657proc rootsModp(int p, poly P)
658"USAGE: rootsModp(p,P);
659RETURN: list of roots of the polynomial P modulo p with p prime
660ASSUME: p>=3
661NOTE:   this algorithm will be called recursively, and it is understood
662        that all the operations are done in Z/pZ (excepting squareRoot(d,p))
663EXAMPLE:example rootsModp; shows an example
664"
665{
666  if(p<3)                                 // (0)[Test if assumptions well-defined]
667  {
668    ERROR("Parameter wrong selected, since p<3!");
669  }
670  else
671  {
672    def S=basering;
673    ring R=p,var(1),dp;
674
675    poly P=imap(S,P);
676    number d;
677    int a;
678    list L;
679
680    poly A=gcd(var(1)^p-var(1),P);        // (1)[Isolate roots in Z/pZ]
681    if(subst(A,var(1),0)==0)
682    {
683      L[1]=0;
684      A=A/var(1);
685    }
686
687    if(deg(A)==0)                         // (2)[Small degree?]
688    {
689      return(L);
690    }
691
692    if(deg(A)==1)
693    {
694      matrix M=coeffs(A,var(1));
696      setring S;
697      list L=imap(R,L);
698      return(L);
699    }
700
701    if(deg(A)==2)
702    {
703      matrix M=coeffs(A,var(1));
705
706      ring T=0,var(1),dp;
707      number d=imap(R,d);
708      number e=squareRoot(bigint(d),bigint(p));
709      setring R;
710      number e=imap(T,e);
711
714      setring S;
715      list L=imap(R,L);
716      return(L);
717    }
718
719    poly B=1;                             // (3)[Random splitting]
720    poly C;
721    while((deg(B)==0)||(deg(B)==deg(A)))
722    {
723      a=random(0,p-1);
724      B=gcd((var(1)+a)^((p-1) div 2)-1,A);
725      C=A/B;
726    }
727
728    setring S;                            // (4)[Recurse]
729    poly B=imap(R,B);
730    poly C=imap(R,C);
731    list l=L+rootsModp(p,B)+rootsModp(p,C);
732    return(l);
733  }
734}
735example
736{ "EXAMPLE:"; echo = 2;
737    ring r = 0,x,dp;
738    poly f=x4+2x3-5x2+x;
739    rootsModp(7,f);
740    poly g=x5+112x4+655x3+551x2+1129x+831;
741    rootsModp(1223,g);
742}
743
744
745proc wUnit(bigint D)
746"USAGE: wUnit(D);
747RETURN: the number of roots of unity in the quadratic order of discriminant D
748ASSUME: D<0 a discriminant kongruent to 0 or 1 modulo 4
749EXAMPLE:example w; shows an example
750"
751{
752  if((D>=0)||((D mod 4)==2)||((D mod 4)==3))
753  {
754    ERROR("Parameter wrong selected!");
755  }
756  else
757  {
758    if(D<-4) {return(2);}
759    if(D==-4){return(4);}
760    if(D==-3){return(6);}
761  }
762}
763example
764{ "EXAMPLE:"; echo = 2;
765    bigint D=-3;
766    wUnit(D);
767}
768
769
770proc Atkin(number N, int K, int B)
771"USAGE: Atkin(N,K,B);
772RETURN:  1, if N is prime,
773        -1, if N is not prime,
774         0, if the algorithm is not applicable, since there are too few discriminants
775ASSUME: N is coprime to 6 and different from 1
776NOTE:   K/2 is input for the procedure "disc",@*
777        K is input for the procedure "HilbertClassPoly",@*
778        B describes the number of recursions being calculated.@*
779        The basis of the algorithm is the following theorem:
780          Let N be an integer coprime to 6 and different from 1 and E be an
781          ellipic curve modulo N.@* Assume that we know an integer m and a
782          point P of E(Z/NZ) satisfying the following conditions.@*
783           (1) There exists a prime divisor q of m such that q > (4-th root(N)+1)^2.@*
784           (2) m*P = O(E) = (0:1:0).@*
785           (3) (m/q)*P = (x:y:t) with t element of (Z/NZ)*.@*
786          Then N is prime.
787EXAMPLE:example Atkin; shows an example
788"
789{
790  if(N==1)           {return(-1);} // (0)[Test if assumptions well-defined]
791  if((N==2)||(N==3)) {return(1);}
792  if(gcd(N,6)!=1)
793  {
794    if(printlevel>=1) {"gcd(N,6) = "+string(gcd(N,6));pause();"";}
795    return(-1);
796  }
797  else
798  {
799    int i;                         // (1)[Initialize]
800    int n(i);
801    number N(i)=N;
802    if(printlevel>=1) {"Set i = 0, n = 0 and N(i) = N(0)= "+string(N(i))+".";pause();"";}
803
804    // declarations:
805    int j(0),j(1),j(2),j(3),j(4),k;    // running indices
806    list L;                            // all primes smaller than 1000
807    list H;                            // sequence of negative discriminants
808    number D;                          // discriminant out of H
809    list L1,L2,S,S1,S2,R;              // lists of relevant elements
810    list P,P1,P2;                      // elliptic points on E(Z/N(i)Z)
811    number m,q;                        // m=|E(Z/N(i)Z)| and q|m
812    number a,b,j,c;                    // characterize E(Z/N(i)Z)
813    number g,u;                        // g out of Z/N(i)Z, u=Jacobi(g,N(i))
814    poly T;                            // T=HilbertClassPoly(D,K)
815    matrix M;                          // M contains the coefficients of T
816
817    if(printlevel>=1) {"List H of possibly suitable discriminants will be calculated.";}
818    H=disc(bigint(N),K div 2);
819    if(printlevel>=1) {"H = "+string(H);pause();"";}
820
821    int step=2;
822    while(1)
823    {
824      if(step==2)                  // (2)[Is N(i) small??]
825      {
826        L=5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997;
827        for(j(0)=1;j(0)<=size(L);j(0)++)
828        {
829          if(N(i)==L[j(0)]){return(1);}
830          if(((N(i) mod L[j(0)])==0)&&(N(i)!=L[j(0)]))
831          {
832            if(printlevel>=1) {"N("+string(i)+") = "+string(N(i))+" is divisible by "+string(L[j(0)])+".";pause();"";}
833            step=14;
834            break;
835          }
836        }
837        if(step==2)
838        {
839          step=3;
840        }
841      }
842
843      if(step==3)                  // (3)[Choose next discriminant]
844      {
845        n(i)=n(i)+1;
846        if(n(i)==size(H)+1)
847        {
848          if(printlevel>=1) {"Algorithm is not applicable, since there are not enough suitable discriminants.";
849                             "Increase the parameter of accuracy K and start the algorithm again.";pause();"";}
850          return(0);
851        }
852        D=H[n(i)];
853        if(printlevel>=1) {"Next discriminant D will be chosen. D = "+string(D)+".";pause();"";}
854        if(Jacobi(D,N(i))!=1)
855        {
856          if(printlevel>=1) {"Jacobi(D,N("+string(i)+")) = "+string(Jacobi(D,N(i)));pause();"";}
857          continue;
858        }
859        else
860        {
861          L1=CornacchiaModified(D,N(i));
862          if(size(L1)>1)
863          {
864            if(printlevel>=1) {"The solution (x,y) of the equation x^2+|D|y^2 = 4N("+string(i)+") is";L1;pause();"";}
865            step=4;
866          }
867          else
868          {
869            if(L1[1]==-1)
870            {
871              if(printlevel>=1) {"The equation x^2+|D|y^2 = 4N("+string(i)+") has no solution.";pause();"";}
872              continue;
873            }
874            if(L1[1]==0)
875            {
876              if(printLevel>=1) {"Algorithm is not applicable for N("+string(i)+") = "+string(N(i))+",";
877                                 "since there are not enough suitable discriminants.";pause();"";}
878              step=14;
879            }
880          }
881        }
882      }
883
884      if(step==4)                  // (4)[Factor m]
885      {
886        if(printlevel>=1) {"List L2 of possible m = |E(Z/N("+string(i)+")Z)| will be calculated.";}
887        if(absValue(L1[1])^2<=4*N(i)) {L2=N(i)+1+L1[1],N(i)+1-L1[1];}
888        if(D==-4)
889        {
890          if(absValue(2*L1[2])^2<=4*N(i)) {L2[size(L2)+1]=N(i)+1+2*L1[2];
891                                           L2[size(L2)+1]=N(i)+1-2*L1[2];}
892        }
893// At this point "<=4*N(i)" has been replaced by "<=16*N(i)".
894        if(D==-3)
895        {
896          if(absValue(L1[1]+3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]+3*L1[2])/2;
897                                                  L2[size(L2)+1]=N(i)+1-(L1[1]+3*L1[2])/2;}
898          if(absValue(L1[1]-3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]-3*L1[2])/2;
899                                                  L2[size(L2)+1]=N(i)+1-(L1[1]-3*L1[2])/2;}
900        }
901///////////////////////////////////////////////////////////////
902        if(size(L2)==0)
903        {
904          if(printlevel>=1) {"Due to the theorem of Hasse there were no possible m = |E(Z/N("+string(i)+")Z)|";
905                             "found for D = "+string(D)+".";}
906          step=3;
907          continue;
908        }
909        else
910        {
911          if(printlevel>=1) {"L2 = ";L2;pause();"";}
912        }
913
914        if(printlevel>=1) {"List S of factors of all possible m will be calculated.";}
915        S=list();
916        for(j(1)=1;j(1)<=size(L2);j(1)++)
917        {
918          m=L2[j(1)];
919          if(m!=0)
920          {
921            S1=PollardRho(m,10000,1,L);
922            S2=pFactor(m,100,L);
923            S[size(S)+1]=list(m,S1+S2);
924          }
925        }
926        if(printlevel>=1) {"S=";S;pause();"";}
927        step=5;
928      }
929
930      if(step==5)                  // (5)[Does a suitable m exist??]
931      {
932        for(j(2)=1;j(2)<=size(S);j(2)++)
933        {
934          m=L2[j(2)];
935          for(j(3)=1;j(3)<=size(S[j(2)][2]);j(3)++)
936          {
937            q=S[j(2)][2][j(3)];
938// sqr(sqr(N(i),50),50) replaces intRoot(intRoot(N(i)))
939            if((q>(sqr(sqr(N(i),50),50)+1)^2) && (MillerRabin(q,5)==1))
940            {
941              step=6;
942              break;
943            }
944//////////////////////////////////////////////////////
945          }
946          if(step==6)
947          {
948            if(printlevel>=1) {"Suitable pair (m,q) has been found such that q|m,";
949                               "q > (4-th root(N("+string(i)+"))+1)^2 and q passes the Miller-Rabin-Test.";
950                               "m = "+string(m)+",";"q = "+string(q);pause();"";}
951            break;
952          }
953          else
954          {
955            step=3;
956          }
957        }
958        if(step==3)
959        {
960          if(printlevel>=1) {"No suitable pair (m,q) has been found such that q|m,";
961                             "q > (4-th root(N("+string(i)+"))+1)^2 and q passes the Miller-Rabin-Test.";
962                              pause();"";}
963          continue;
964        }
965      }
966
967      if(step==6)                  // (6)[Compute elliptic curve]
968      {
969        if(D==-4)
970        {
971          a=-1;
972          b=0;
973          if(printlevel>=1) {"Since D = -4, set a = -1 and b = 0.";pause();"";}
974        }
975        if(D==-3)
976        {
977          a=0;
978          b=-1;
979          if(printlevel>=1) {"Since D = -3, set a = 0 and b = -1.";pause();"";}
980        }
981        if(D<-4)
982        {
983          if(printlevel>=1) {"The minimal polynomial T of j((D+sqr(D))/2) in Z[X] will be calculated for D="+string(D)+".";}
984          T=HilbertClassPoly(D,K);
985          if(printlevel>=1) {"T = "+string(T);pause();"";}
986
987          M=coeffs(T,var(1));
988          T=0;
989
990          for(j(4)=1;j(4)<=nrows(M);j(4)++)
991          {
992            M[j(4),1]=leadcoef(M[j(4),1]) mod N(i);
993            T=T+M[j(4),1]*var(1)^(j(4)-1);
994          }
995          if(printlevel>=1) {"Set T = T mod N("+string(i)+").";"T = "+string(T);pause();"";}
996
997          R=rootsModp(int(N(i)),T);
998          if(deg(T)>size(R))
999          {
1000            ERROR("The polynomial T does not completely split into linear factors modulo N("+string(i)+")."
1001                  "Increase the parameter of accuracy K and start the algorithm again.");
1002          }
1003          if(printlevel>=1) {if(deg(T)>1) {"The "+string(deg(T))+" zeroes of T modulo N("+string(i)+") are";
1004                                           R;pause();"";}
1005                             if(deg(T)==1){"The zero of T modulo N("+string(i)+") is";R;pause();"";}}
1006
1007          j=R[1];
1008          c=j*exgcdN(j-1728,N(i))[1];
1009          a=-3*c mod N(i);
1010          b=2*c mod N(i);
1011          if(printlevel>=1) {"Choose the zero j = "+string(j)+" and set"; "c = j/(j-1728) mod N("+string(i)+"), a = -3c mod N("+string(i)+"), b = 2c mod N("+string(i)+").";
1012                             "a = "+string(a)+",";"b = "+string(b);pause();"";}
1013        }
1014        step=7;
1015      }
1016
1017      if(step==7)                  // (7)[Find g]
1018      {
1019        if(D==-3)
1020        {
1021          while(1)
1022          {
1023            g=random(1,2147483647) mod N(i);
1024            u=Jacobi(g,N(i));
1025            if((u==-1)&&(powerN(g,(N(i)-1)/3,N(i))!=1))
1026            {
1027              if(printlevel>=1) {"g = "+string(g);pause();"";}
1028              break;
1029            }
1030          }
1031        }
1032        else
1033        {
1034          while(1)
1035          {
1036            g=random(1,2147483647) mod N(i);
1037            u=Jacobi(g,N(i));
1038            if(u==-1)
1039            {
1040              if(printlevel>=1) {"g = "+string(g);pause();"";}
1041              break;
1042            }
1043          }
1044        }
1045        step=8;
1046      }
1047
1048      if(step==8)                  // (8)[Find P]
1049      {
1050        if(printlevel>=1) {"A random point P on the elliptic curve corresponding";
1051                           "to the equation y^2 = x^3+ax+b for";"N("+string(i)+") = "+string(N(i))+",";
1052                           "   a = "+string(a)+",";"   b = "+string(b);"will be chosen.";}
1053        P=ellipticRandomPoint(N(i),a,b);
1054        if(printlevel>=1) {"P = ("+string(P)+")";pause();"";}
1055
1056        if(size(P)==1)
1057        {
1058          step=14;
1059        }
1060        else
1061        {
1062          step=9;
1063        }
1064      }
1065
1066      if(step==9)                  // (9)[Find right curve]
1067      {
1068        if(printlevel>=1) {"The points P2 = (m/q)*P and P1 = q*P2 on the curve will be calculated.";}
1069        P2=ellipticMult(N(i),a,b,P,m/q);
1070        P1=ellipticMult(N(i),a,b,P2,q);
1071        if(printlevel>=1) {"P1 = ("+string(P1)+"),";"P2 = ("+string(P2)+")";pause();"";}
1072
1073        if((P1[1]==0)&&(P1[2]==1)&&(P1[3]==0))
1074        {
1075          step=12;
1076        }
1077        else
1078        {
1079          if(printlevel>=1) {"Since P1 != (0:1:0), it holds m != |E(Z/N("+string(i)+")Z)| for the coefficients a = "+string(a)+" and b = "+string(b)+".";
1080                             "Therefore choose new coefficients a and b.";pause();"";}
1081          step=10;
1082        }
1083      }
1084
1085      if(step==10)                 // (10)[Change coefficients]
1086      {
1087        k=k+1;
1088        if(k>=wUnit(bigint(D)))
1089        {
1090          if(printlevel>=1) {"Since k = wUnit(D) = "+string(k)+", it holds that N("+string(i)+") = "+string(N(i))+" is not prime.";pause();"";}
1091          step=14;
1092        }
1093        else
1094        {
1095          if(D<-4) {a=a*g^2 mod N(i); b=b*g^3 mod N(i);
1096                    if(printlevel>=1) {"Since D < -4, set a = a*g^2 mod N("+string(i)+") and b = b*g^3 mod N("+string(i)+").";
1097                                       "a = "+string(a)+",";"b = "+string(b)+",";"k = "+string(k);pause();"";}}
1098          if(D==-4){a=a*g mod N(i);
1099                    if(printlevel>=1) {"Since D = -4, set a = a*g mod N("+string(i)+").";"a = "+string(a)+",";
1100                                       "b = "+string(b)+",";"k = "+string(k);pause();"";}}
1101          if(D==-3){b=b*g mod N(i);
1102                    if(printlevel>=1) {"Since D = -3, set b = b*g mod N("+string(i)+").";"a = "+string(a)+",";
1103                                       "b = "+string(b)+",";"k = "+string(k);pause();"";}}
1104          step=8;
1105          continue;
1106        }
1107      }
1108
1109      if(step==11)                 // (11)[Find a new P]
1110      {
1111        if(printlevel>=1) {"A new random point P on the elliptic curve will be chosen,";
1112                           "since also P2 = (0:1:0).";}
1113        P=ellipticRandomPoint(N(i),a,b);
1114        if(printlevel>=1) {"P = ("+string(P)+")";pause();"";}
1115
1116        if(size(P)==1)
1117        {
1118          step=14;
1119        }
1120        else
1121        {
1122          if(printlevel>=1) {"The points P2 = (m/q)*P and P1 = q*P2 on the curve will be calculated.";}
1123          P2=ellipticMult(N(i),a,b,P,m/q);
1124          P1=ellipticMult(N(i),a,b,P2,q);
1125          if(printlevel>=1) {"P1 = ("+string(P1)+"),";"P2 = ("+string(P2)+")";pause();"";}
1126
1127          if((P1[1]!=0)||(P1[2]!=1)||(P1[3]!=0))
1128          {
1129            if(printlevel>=1) {"Since P1 != (0:1:0), it holds m != |E(Z/N("+string(i)+")Z)| for the coefficients a = "+string(a)+" and b = "+string(b)+".";
1130                               "Therefore choose new coefficients a and b.";pause();"";}
1131            step=10;
1132            continue;
1133          }
1134          else
1135          {
1136            step=12;
1137          }
1138        }
1139      }
1140
1141      if(step==12)                 // (12)[Check P]
1142      {
1143        if((P2[1]==0)&&(P2[2]==1)&&(P2[3]==0))
1144        {
1145          step=11;
1146          continue;
1147        }
1148        else
1149        {
1150          step=13;
1151        }
1152      }
1153
1154      if(step==13)                 // (13)[Recurse]
1155      {
1156        if(i<B)
1157        {
1158          if(printlevel>=1) {string(i+1)+". Recursion:";"";
1159                             "N("+string(i)+") = "+string(N(i))+" suffices the conditions of the underlying theorem,";
1160                             "since P1 = (0:1:0) and P2[3] in (Z/N("+string(i)+")Z)*.";"";
1161                             "Now check if also the found factor q="+string(q)+" suffices these assumptions.";
1162                             "Therefore set i = i+1, N("+string(i+1)+") = q = "+string(q)+" and restart the algorithm.";pause();"";}
1163          i=i+1;
1164          int n(i);
1165          number N(i)=q;
1166          k=0;
1167          step=2;
1168          continue;
1169        }
1170        else
1171        {
1172          if(printlevel>=1) {"N(B) = N("+string(i)+") = "+string(N(i))+" suffices the conditions of the underlying theorem,";
1173                             "since P1 = (0:1:0) and P2[3] in (Z/N("+string(i)+")Z)*.";
1174                             "In particular N = "+string(N)+" is prime.";pause();"";}
1175          return(1);
1176        }
1177      }
1178
1179      if(step==14)                 // (14)[Backtrack]
1180      {
1181        if(i>0)
1182        {
1183          if(printlevel>=1) {"Set i = i-1 and restart the algorithm for N("+string(i-1)+") = "+string(N(i-1))+" with";
1184                             "a new discriminant.";pause();"";}
1185          i=i-1;
1186          k=0;
1187          step=3;
1188        }
1189        else
1190        {
1191          if(printlevel>=1) {"N(0) = N = "+string(N)+" and therefore N is not prime.";pause();"";}
1192          return(-1);
1193        }
1194      }
1195    }
1196  }
1197}
1198example
1199{ "EXAMPLE:"; echo = 2;
1200    ring R = 0,x,dp;
1201    Atkin(7691,100,5);
1202    Atkin(3473,10,2);
1203    printlevel=1;
1204    Atkin(10000079,100,2);
1205}
1206
Note: See TracBrowser for help on using the repository browser.