source: git/Singular/LIB/atkins.lib @ f8df35

fieker-DuValspielwiese
Last change on this file since f8df35 was f8df35, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: more bigint instead of number
  • Property mode set to 100644
File size: 32.4 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version atkins.lib 4.0.0.1 Jun_2014 ";
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)==-1)||(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          return(-1);
251          // ERROR("The Diophantine equation has no solution!");
252        }
253        else
254        {
255          list L=b,root_c;
256          return(L);
257        }
258      }
259    }
260  }
261}
262example
263{ "EXAMPLE:"; echo = 2;
264    CornacchiaModified(-107,1319);
265}
266
267
268static proc pFactor1(number n,int B, list P)
269"USAGE: pFactor1(n,B,P); n to be factorized, B a bound , P a list of primes
270RETURN: a list of factors of n or the message: no factor found
271NOTE:   Pollard's p-factorization
272        creates the product k of powers of primes (bounded by B)  from
273        the list P with the idea that for a prime divisor p of n p-1|k
274        then p devides gcd(a^k-1,n) for some random a
275EXAMPLE:example pFactor1; shows an example
276"
277{
278  int i;
279  number k=1;
280  number w;
281  while(i<size(P))
282  {
283    i++;
284    w=P[i];
285    if(w>B) {break;}
286    while(w*P[i]<=B)
287    {
288      w=w*P[i];
289    }
290    k=k*w;
291  }
292  number a=random(2,2147483629);
293  number d=gcd(powerN(a,k,n)-1,n);
294  if((d>1)&&(d<n)){return(d);}
295  return(n);
296}
297example
298{ "EXAMPLE:"; echo = 2;
299   ring R = 0,z,dp;
300   list L=primList(1000);
301   pFactor1(1241143,13,L);
302   number h=10;
303   h=h^30+25;
304   pFactor1(h,20,L);
305}
306
307
308proc maximum(list L)
309"USAGE: maximum(list L);
310RETURN: the maximal number contained in list L
311EXAMPLE:example maximum; shows an example
312"
313{
314  number max=L[1];
315  int i;
316  for(i=2;i<=size(L);i++)
317  {
318    if(L[i]>max)
319    {
320      max=L[i];
321    }
322  }
323  return(max);
324}
325example
326{ "EXAMPLE:"; echo = 2;
327    ring r = 0,x,dp;
328    list L=465,867,1233,4567,776544,233445,2334,556;
329    maximum(L);
330}
331
332
333static proc cmod(number x, number y)
334"USAGE: cmod(x,y);
335RETURN: x mod y
336ASSUME: x,y out of Z and x,y<=2147483647
337NOTE:   this algorithm is a helping procedure to be able to calculate
338        x mod y with x,y out of Z while working in the complex field
339EXAMPLE:example cmod; shows an example
340"
341{
342  int rest=int(x-y*int(x/y));
343  if(rest<0)
344  {
345    rest=rest+int(y);
346  }
347  return(rest);
348}
349example
350{ "EXAMPLE:"; echo = 2;
351    ring r = (complex,30,i),x,dp;
352    number x=-1004456;
353    number y=1233;
354    cmod(x,y);
355}
356
357
358proc sqr(number w, int k)
359"USAGE: sqr(w,k);
360RETURN: the square root of w
361ASSUME: w>=0
362NOTE:   k describes the number of decimals being calculated in the real numbers,
363        k, intPart(k/5) are inputs for the procedure "nt_solve"
364EXAMPLE:example sqr; shows an example
365"
366{
367  poly f=var(1)^2-w;
368  def S=basering;
369  ring R=(real,k),var(1),dp;
370  poly f=imap(S,f);
371  ideal I=nt_solve(f,1.1,list(k,k div 5));
372  number c=leadcoef(I[1]);
373  setring S;
374  number c=imap(R,c);
375  return(c);
376}
377example
378{ "EXAMPLE:"; echo = 2;
379    ring R = (real,60),x,dp;
380    number ww=288469650108669535726081;
381    sqr(ww,60);
382}
383
384
385proc expo(number z, int k)
386"USAGE: expo(z,k);
387RETURN: e^z to the order k
388NOTE:   k describes the number of summands being calculated in the exponential power series
389EXAMPLE:example expo; shows an example
390"
391{
392  number q=1;
393  number e=1;
394  int n;
395  for(n=1;n<=k;n++)
396  {
397    q=q*z/n;
398    e=e+q;
399  }
400  return(e);
401}
402example
403{ "EXAMPLE:"; echo = 2;
404    ring r = (real,30),x,dp;
405    number z=40.35;
406    expo(z,1000);
407}
408
409
410proc jOft(number t, int k)
411"USAGE: jOft(t,k);
412RETURN: the j-invariant of t
413ASSUME: t is a complex number with positive imaginary part
414NOTE:   k describes the number of summands being calculated in the power series,
415        10*k is input for the procedure @code{expo}
416EXAMPLE:example jOft; shows an example
417"
418{
419  number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j;
420
421  number pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989;
422
423  tr=repart(t);
424  ti=impart(t);
425  if(tr==-1/2){qr1=-1;}
426  else
427  {
428    if(tr==0){qr1=1;}
429    else
430    {
431      tr=tr-round(tr);
432      qr1=expo(2*i*pi*tr,10*k);
433    }
434  }
435
436  qi1=expo(-pi*ti,10*k);
437  q1=qr1*qi1^2;
438  q2=q1^2;
439
440  int n=1;
441  while(n<=k)
442  {
443    m1=m1+(-1)^n*(q1^(n*(3*n-1) div 2)+q1^(n*(3*n+1) div 2));
444    m2=m2+(-1)^n*(q2^(n*(3*n-1) div 2)+q2^(n*(3*n+1) div 2));
445    n++;
446  }
447
448  f=q1*((1+m2)/(1+m1))^24;
449
450  j=(256*f+1)^3/f;
451  return(j);
452}
453example
454{ "EXAMPLE:"; echo = 2;
455    ring r = (complex,30,i),x,dp;
456    number t=(-7+i*sqr(7,250))/2;
457    jOft(t,50);
458}
459
460
461proc round(number r)
462"USAGE: round(r);
463RETURN: the nearest number to r out of Z
464ASSUME: r should be a rational or a real number
465EXAMPLE:example round; shows an example
466"
467{
468  number a=absValue(r);
469  number v=r/a;
470
471  number d=10;
472  int e;
473  while(1)
474  {
475    e=e+1;
476    if(a-d^e<0)
477    {
478      e=e-1;
479      break;
480    }
481  }
482
483  number b=a;
484  int k;
485  for(k=0;k<=e;k++)
486  {
487    while(1)
488    {
489      b=b-d^(e-k);
490      if(b<0)
491      {
492        b=b+d^(e-k);
493        break;
494      }
495    }
496  }
497
498  if(b<1/2)
499  {
500    return(v*(a-b));
501  }
502  else
503  {
504    return(v*(a+1-b));
505  }
506}
507example
508{ "EXAMPLE:"; echo = 2;
509    ring R = (real,50),x,dp;
510    number r=7357683445788723456321.6788643224;
511    round(r);
512}
513
514
515proc HilbertClassPoly(bigint D, int k)
516"USAGE: HilbertClassPoly(D,k);
517RETURN: the monic polynomial of degree h(D) in Z[X] of which jOft((D+sqr(D))/2) is a root
518ASSUME: D is a negative discriminant
519NOTE:   k is input for the procedure "jOft",
520        5*k is input for the procedure "sqr",
521        10*k describes the number of decimals being calculated in the complex numbers
522EXAMPLE:example HilbertClassPoly; shows an example
523"
524{
525  if(D>=0)                     // (0)[Test if assumptions well-defined]
526  {
527    ERROR("Parameter wrong selected!");
528  }
529  else
530  {
531    def S=basering;
532
533    string s1,s2,s3;
534    bigint B=intRoot(absValue(D) div 3);
535
536    ring C=(complex,10*k,i),x,dp;
537    number DD=D;
538
539    poly P=1;                  // (1)[Initialize]
540    number b=cmod(DD,2);
541
542    number t,a,g,tau,j;
543    list L;
544
545    bigint a1,b1,t1,g1;
546    int step=2;
547    while(1)
548    {
549      if(step==2)              // (2)[Initialize a]
550      {
551        t=(b^2-DD)/4;
552        L=b,1;
553        a=maximum(L);
554        step=3;
555      }
556
557      if(step==3)              // (3)[Test]
558      {
559        if((cmod(t,a)!=0))
560        {
561          step=4;
562        }
563        else
564        {
565          s1=string(a);
566          s2=string(b);
567          s3=string(t);
568
569          execute("a1="+s1+";");
570          execute("b1="+s2+";");
571          execute("t1="+s3+";");
572          g1=gcd(gcd(a1,b1),t1 div a1);
573          setring C;
574          g=g1;
575
576          if(g!=1)
577          {
578            step=4;
579          }
580          else
581          {
582            tau=(-b+i*sqr(absValue(DD),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  }
647}
648example
649{ "EXAMPLE:"; echo = 2;
650    ring r = 0,x,dp;
651    bigint D=-23;
652    HilbertClassPoly(D,50);
653}
654
655
656proc rootsModp(int p, poly P)
657"USAGE: rootsModp(p,P);
658RETURN: list of roots of the polynomial P modulo p with p prime
659ASSUME: p>=3
660NOTE:   this algorithm will be called recursively, and it is understood
661        that all the operations are done in Z/pZ (excepting squareRoot(d,p))
662EXAMPLE:example rootsModp; shows an example
663"
664{
665  if(p<3)                                 // (0)[Test if assumptions well-defined]
666  {
667    ERROR("Parameter wrong selected, since p<3!");
668  }
669  else
670  {
671    def S=basering;
672    ring R=p,var(1),dp;
673
674    poly P=imap(S,P);
675    number d;
676    int a;
677    list L;
678
679    poly A=gcd(var(1)^p-var(1),P);        // (1)[Isolate roots in Z/pZ]
680    if(subst(A,var(1),0)==0)
681    {
682      L[1]=0;
683      A=A/var(1);
684    }
685
686    if(deg(A)==0)                         // (2)[Small degree?]
687    {
688      return(L);
689    }
690
691    if(deg(A)==1)
692    {
693      matrix M=coeffs(A,var(1));
694      L[size(L)+1]=-leadcoef(M[1,1])/leadcoef(M[2,1]);
695      setring S;
696      list L=imap(R,L);
697      return(L);
698    }
699
700    if(deg(A)==2)
701    {
702      matrix M=coeffs(A,var(1));
703      d=leadcoef(M[2,1])^2-4*leadcoef(M[1,1])*leadcoef(M[3,1]);
704
705      ring T=0,var(1),dp;
706      number d=imap(R,d);
707      number e=squareRoot(bigint(d),bigint(p));
708      setring R;
709      number e=imap(T,e);
710
711      L[size(L)+1]=(-leadcoef(M[2,1])+e)/(2*leadcoef(M[3,1]));
712      L[size(L)+1]=(-leadcoef(M[2,1])-e)/(2*leadcoef(M[3,1]));
713      setring S;
714      list L=imap(R,L);
715      return(L);
716    }
717
718    poly B=1;                             // (3)[Random splitting]
719    poly C;
720    while((deg(B)==0)||(deg(B)==deg(A)))
721    {
722      a=random(0,p-1);
723      B=gcd((var(1)+a)^((p-1) div 2)-1,A);
724      C=A/B;
725    }
726
727    setring S;                            // (4)[Recurse]
728    poly B=imap(R,B);
729    poly C=imap(R,C);
730    list l=L+rootsModp(p,B)+rootsModp(p,C);
731    return(l);
732  }
733}
734example
735{ "EXAMPLE:"; echo = 2;
736    ring r = 0,x,dp;
737    poly f=x4+2x3-5x2+x;
738    rootsModp(7,f);
739    poly g=x5+112x4+655x3+551x2+1129x+831;
740    rootsModp(1223,g);
741}
742
743
744proc wUnit(bigint D)
745"USAGE: wUnit(D);
746RETURN: the number of roots of unity in the quadratic order of discriminant D
747ASSUME: D<0 a discriminant kongruent to 0 or 1 modulo 4
748EXAMPLE:example w; shows an example
749"
750{
751  if((D>=0)||((D mod 4)==2)||((D mod 4)==3))
752  {
753    ERROR("Parameter wrong selected!");
754  }
755  else
756  {
757    if(D<-4) {return(2);}
758    if(D==-4){return(4);}
759    if(D==-3){return(6);}
760  }
761}
762example
763{ "EXAMPLE:"; echo = 2;
764    bigint D=-3;
765    wUnit(D);
766}
767
768
769proc Atkin(number N, int K, int B)
770"USAGE: Atkin(N,K,B);
771RETURN:  1, if N is prime,
772        -1, if N is not prime,
773         0, if the algorithm is not applicable, since there are too few discriminants
774ASSUME: N is coprime to 6 and different from 1
775NOTE:   K/2 is input for the procedure "disc",@*
776        K is input for the procedure "HilbertClassPoly",@*
777        B describes the number of recursions being calculated.@*
778        The basis of the algorithm is the following theorem:
779          Let N be an integer coprime to 6 and different from 1 and E be an
780          ellipic curve modulo N.@* Assume that we know an integer m and a
781          point P of E(Z/NZ) satisfying the following conditions.@*
782           (1) There exists a prime divisor q of m such that q > (4-th root(N)+1)^2.@*
783           (2) m*P = O(E) = (0:1:0).@*
784           (3) (m/q)*P = (x:y:t) with t element of (Z/NZ)*.@*
785          Then N is prime.
786EXAMPLE:example Atkin; shows an example
787"
788{
789  if(N==1)           {return(-1);} // (0)[Test if assumptions well-defined]
790  if((N==2)||(N==3)) {return(1);}
791  if(gcd(N,6)!=1)
792  {
793    if(printlevel>=1) {"gcd(N,6) = "+string(gcd(N,6));pause();"";}
794    return(-1);
795  }
796  else
797  {
798    int i;                         // (1)[Initialize]
799    int n(i);
800    number N(i)=N;
801    if(printlevel>=1) {"Set i = 0, n = 0 and N(i) = N(0)= "+string(N(i))+".";pause();"";}
802
803    // declarations:
804    int j(0),j(1),j(2),j(3),j(4),k;    // running indices
805    list L;                            // all primes smaller than 1000
806    list H;                            // sequence of negative discriminants
807    number D;                          // discriminant out of H
808    list L1,L2,S,S1,S2,R;              // lists of relevant elements
809    list P,P1,P2;                      // elliptic points on E(Z/N(i)Z)
810    number m,q;                        // m=|E(Z/N(i)Z)| and q|m
811    number a,b,j,c;                    // characterize E(Z/N(i)Z)
812    number g,u;                        // g out of Z/N(i)Z, u=Jacobi(g,N(i))
813    poly T;                            // T=HilbertClassPoly(D,K)
814    matrix M;                          // M contains the coefficients of T
815
816    if(printlevel>=1) {"List H of possibly suitable discriminants will be calculated.";}
817    H=disc(bigint(N),K div 2);
818    if(printlevel>=1) {"H = "+string(H);pause();"";}
819
820    int step=2;
821    while(1)
822    {
823      if(step==2)                  // (2)[Is N(i) small??]
824      {
825        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;
826        for(j(0)=1;j(0)<=size(L);j(0)++)
827        {
828          if(N(i)==L[j(0)]){return(1);}
829          if(((N(i) mod L[j(0)])==0)&&(N(i)!=L[j(0)]))
830          {
831            if(printlevel>=1) {"N("+string(i)+") = "+string(N(i))+" is divisible by "+string(L[j(0)])+".";pause();"";}
832            step=14;
833            break;
834          }
835        }
836        if(step==2)
837        {
838          step=3;
839        }
840      }
841
842      if(step==3)                  // (3)[Choose next discriminant]
843      {
844        n(i)=n(i)+1;
845        if(n(i)==size(H)+1)
846        {
847          if(printlevel>=1) {"Algorithm is not applicable, since there are not enough suitable discriminants.";
848                             "Increase the parameter of accuracy K and start the algorithm again.";pause();"";}
849          return(0);
850        }
851        D=H[n(i)];
852        if(printlevel>=1) {"Next discriminant D will be chosen. D = "+string(D)+".";pause();"";}
853        if(Jacobi(D,N(i))!=1)
854        {
855          if(printlevel>=1) {"Jacobi(D,N("+string(i)+")) = "+string(Jacobi(D,N(i)));pause();"";}
856          continue;
857        }
858        else
859        {
860          L1=CornacchiaModified(D,N(i));
861          if(size(L1)>1)
862          {
863            if(printlevel>=1) {"The solution (x,y) of the equation x^2+|D|y^2 = 4N("+string(i)+") is";L1;pause();"";}
864            step=4;
865          }
866          else
867          {
868            if(L1[1]==-1)
869            {
870              if(printlevel>=1) {"The equation x^2+|D|y^2 = 4N("+string(i)+") has no solution.";pause();"";}
871              continue;
872            }
873            if(L1[1]==0)
874            {
875              if(printLevel>=1) {"Algorithm is not applicable for N("+string(i)+") = "+string(N(i))+",";
876                                 "since there are not enough suitable discriminants.";pause();"";}
877              step=14;
878            }
879          }
880        }
881      }
882
883      if(step==4)                  // (4)[Factor m]
884      {
885        if(printlevel>=1) {"List L2 of possible m = |E(Z/N("+string(i)+")Z)| will be calculated.";}
886        if(absValue(L1[1])^2<=4*N(i)) {L2=N(i)+1+L1[1],N(i)+1-L1[1];}
887        if(D==-4)
888        {
889          if(absValue(2*L1[2])^2<=4*N(i)) {L2[size(L2)+1]=N(i)+1+2*L1[2];
890                                           L2[size(L2)+1]=N(i)+1-2*L1[2];}
891        }
892// At this point "<=4*N(i)" has been replaced by "<=16*N(i)".
893        if(D==-3)
894        {
895          if(absValue(L1[1]+3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]+3*L1[2])/2;
896                                                  L2[size(L2)+1]=N(i)+1-(L1[1]+3*L1[2])/2;}
897          if(absValue(L1[1]-3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]-3*L1[2])/2;
898                                                  L2[size(L2)+1]=N(i)+1-(L1[1]-3*L1[2])/2;}
899        }
900///////////////////////////////////////////////////////////////
901        if(size(L2)==0)
902        {
903          if(printlevel>=1) {"Due to the theorem of Hasse there were no possible m = |E(Z/N("+string(i)+")Z)|";
904                             "found for D = "+string(D)+".";}
905          step=3;
906          continue;
907        }
908        else
909        {
910          if(printlevel>=1) {"L2 = ";L2;pause();"";}
911        }
912
913        if(printlevel>=1) {"List S of factors of all possible m will be calculated.";}
914        S=list();
915        for(j(1)=1;j(1)<=size(L2);j(1)++)
916        {
917          m=L2[j(1)];
918          if(m!=0)
919          {
920            S1=PollardRho(m,10000,1,L);
921            S2=pFactor(m,100,L);
922            S[size(S)+1]=list(m,S1+S2);
923          }
924        }
925        if(printlevel>=1) {"S=";S;pause();"";}
926        step=5;
927      }
928
929      if(step==5)                  // (5)[Does a suitable m exist??]
930      {
931        for(j(2)=1;j(2)<=size(S);j(2)++)
932        {
933          m=L2[j(2)];
934          for(j(3)=1;j(3)<=size(S[j(2)][2]);j(3)++)
935          {
936            q=S[j(2)][2][j(3)];
937// sqr(sqr(N(i),50),50) replaces intRoot(intRoot(N(i)))
938            if((q>(sqr(sqr(N(i),50),50)+1)^2) && (MillerRabin(q,5)==1))
939            {
940              step=6;
941              break;
942            }
943//////////////////////////////////////////////////////
944          }
945          if(step==6)
946          {
947            if(printlevel>=1) {"Suitable pair (m,q) has been found such that q|m,";
948                               "q > (4-th root(N("+string(i)+"))+1)^2 and q passes the Miller-Rabin-Test.";
949                               "m = "+string(m)+",";"q = "+string(q);pause();"";}
950            break;
951          }
952          else
953          {
954            step=3;
955          }
956        }
957        if(step==3)
958        {
959          if(printlevel>=1) {"No suitable pair (m,q) has been found such that q|m,";
960                             "q > (4-th root(N("+string(i)+"))+1)^2 and q passes the Miller-Rabin-Test.";
961                              pause();"";}
962          continue;
963        }
964      }
965
966      if(step==6)                  // (6)[Compute elliptic curve]
967      {
968        if(D==-4)
969        {
970          a=-1;
971          b=0;
972          if(printlevel>=1) {"Since D = -4, set a = -1 and b = 0.";pause();"";}
973        }
974        if(D==-3)
975        {
976          a=0;
977          b=-1;
978          if(printlevel>=1) {"Since D = -3, set a = 0 and b = -1.";pause();"";}
979        }
980        if(D<-4)
981        {
982          if(printlevel>=1) {"The minimal polynomial T of j((D+sqr(D))/2) in Z[X] will be calculated for D="+string(D)+".";}
983          T=HilbertClassPoly(D,K);
984          if(printlevel>=1) {"T = "+string(T);pause();"";}
985
986          M=coeffs(T,var(1));
987          T=0;
988
989          for(j(4)=1;j(4)<=nrows(M);j(4)++)
990          {
991            M[j(4),1]=leadcoef(M[j(4),1]) mod N(i);
992            T=T+M[j(4),1]*var(1)^(j(4)-1);
993          }
994          if(printlevel>=1) {"Set T = T mod N("+string(i)+").";"T = "+string(T);pause();"";}
995
996          R=rootsModp(int(N(i)),T);
997          if(deg(T)>size(R))
998          {
999            ERROR("The polynomial T does not completely split into linear factors modulo N("+string(i)+")."
1000                  "Increase the parameter of accuracy K and start the algorithm again.");
1001          }
1002          if(printlevel>=1) {if(deg(T)>1) {"The "+string(deg(T))+" zeroes of T modulo N("+string(i)+") are";
1003                                           R;pause();"";}
1004                             if(deg(T)==1){"The zero of T modulo N("+string(i)+") is";R;pause();"";}}
1005
1006          j=R[1];
1007          c=j*exgcdN(j-1728,N(i))[1];
1008          a=-3*c mod N(i);
1009          b=2*c mod N(i);
1010          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)+").";
1011                             "a = "+string(a)+",";"b = "+string(b);pause();"";}
1012        }
1013        step=7;
1014      }
1015
1016      if(step==7)                  // (7)[Find g]
1017      {
1018        if(D==-3)
1019        {
1020          while(1)
1021          {
1022            g=random(1,2147483647) mod N(i);
1023            u=Jacobi(g,N(i));
1024            if((u==-1)&&(powerN(g,(N(i)-1)/3,N(i))!=1))
1025            {
1026              if(printlevel>=1) {"g = "+string(g);pause();"";}
1027              break;
1028            }
1029          }
1030        }
1031        else
1032        {
1033          while(1)
1034          {
1035            g=random(1,2147483647) mod N(i);
1036            u=Jacobi(g,N(i));
1037            if(u==-1)
1038            {
1039              if(printlevel>=1) {"g = "+string(g);pause();"";}
1040              break;
1041            }
1042          }
1043        }
1044        step=8;
1045      }
1046
1047      if(step==8)                  // (8)[Find P]
1048      {
1049        if(printlevel>=1) {"A random point P on the elliptic curve corresponding";
1050                           "to the equation y^2 = x^3+ax+b for";"N("+string(i)+") = "+string(N(i))+",";
1051                           "   a = "+string(a)+",";"   b = "+string(b);"will be chosen.";}
1052        P=ellipticRandomPoint(N(i),a,b);
1053        if(printlevel>=1) {"P = ("+string(P)+")";pause();"";}
1054
1055        if(size(P)==1)
1056        {
1057          step=14;
1058        }
1059        else
1060        {
1061          step=9;
1062        }
1063      }
1064
1065      if(step==9)                  // (9)[Find right curve]
1066      {
1067        if(printlevel>=1) {"The points P2 = (m/q)*P and P1 = q*P2 on the curve will be calculated.";}
1068        P2=ellipticMult(N(i),a,b,P,m/q);
1069        P1=ellipticMult(N(i),a,b,P2,q);
1070        if(printlevel>=1) {"P1 = ("+string(P1)+"),";"P2 = ("+string(P2)+")";pause();"";}
1071
1072        if((P1[1]==0)&&(P1[2]==1)&&(P1[3]==0))
1073        {
1074          step=12;
1075        }
1076        else
1077        {
1078          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)+".";
1079                             "Therefore choose new coefficients a and b.";pause();"";}
1080          step=10;
1081        }
1082      }
1083
1084      if(step==10)                 // (10)[Change coefficients]
1085      {
1086        k=k+1;
1087        if(k>=wUnit(bigint(D)))
1088        {
1089          if(printlevel>=1) {"Since k = wUnit(D) = "+string(k)+", it holds that N("+string(i)+") = "+string(N(i))+" is not prime.";pause();"";}
1090          step=14;
1091        }
1092        else
1093        {
1094          if(D<-4) {a=a*g^2 mod N(i); b=b*g^3 mod N(i);
1095                    if(printlevel>=1) {"Since D < -4, set a = a*g^2 mod N("+string(i)+") and b = b*g^3 mod N("+string(i)+").";
1096                                       "a = "+string(a)+",";"b = "+string(b)+",";"k = "+string(k);pause();"";}}
1097          if(D==-4){a=a*g mod N(i);
1098                    if(printlevel>=1) {"Since D = -4, set a = a*g mod N("+string(i)+").";"a = "+string(a)+",";
1099                                       "b = "+string(b)+",";"k = "+string(k);pause();"";}}
1100          if(D==-3){b=b*g mod N(i);
1101                    if(printlevel>=1) {"Since D = -3, set b = b*g mod N("+string(i)+").";"a = "+string(a)+",";
1102                                       "b = "+string(b)+",";"k = "+string(k);pause();"";}}
1103          step=8;
1104          continue;
1105        }
1106      }
1107
1108      if(step==11)                 // (11)[Find a new P]
1109      {
1110        if(printlevel>=1) {"A new random point P on the elliptic curve will be chosen,";
1111                           "since also P2 = (0:1:0).";}
1112        P=ellipticRandomPoint(N(i),a,b);
1113        if(printlevel>=1) {"P = ("+string(P)+")";pause();"";}
1114
1115        if(size(P)==1)
1116        {
1117          step=14;
1118        }
1119        else
1120        {
1121          if(printlevel>=1) {"The points P2 = (m/q)*P and P1 = q*P2 on the curve will be calculated.";}
1122          P2=ellipticMult(N(i),a,b,P,m/q);
1123          P1=ellipticMult(N(i),a,b,P2,q);
1124          if(printlevel>=1) {"P1 = ("+string(P1)+"),";"P2 = ("+string(P2)+")";pause();"";}
1125
1126          if((P1[1]!=0)||(P1[2]!=1)||(P1[3]!=0))
1127          {
1128            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)+".";
1129                               "Therefore choose new coefficients a and b.";pause();"";}
1130            step=10;
1131            continue;
1132          }
1133          else
1134          {
1135            step=12;
1136          }
1137        }
1138      }
1139
1140      if(step==12)                 // (12)[Check P]
1141      {
1142        if((P2[1]==0)&&(P2[2]==1)&&(P2[3]==0))
1143        {
1144          step=11;
1145          continue;
1146        }
1147        else
1148        {
1149          step=13;
1150        }
1151      }
1152
1153      if(step==13)                 // (13)[Recurse]
1154      {
1155        if(i<B)
1156        {
1157          if(printlevel>=1) {string(i+1)+". Recursion:";"";
1158                             "N("+string(i)+") = "+string(N(i))+" suffices the conditions of the underlying theorem,";
1159                             "since P1 = (0:1:0) and P2[3] in (Z/N("+string(i)+")Z)*.";"";
1160                             "Now check if also the found factor q="+string(q)+" suffices these assumptions.";
1161                             "Therefore set i = i+1, N("+string(i+1)+") = q = "+string(q)+" and restart the algorithm.";pause();"";}
1162          i=i+1;
1163          int n(i);
1164          number N(i)=q;
1165          k=0;
1166          step=2;
1167          continue;
1168        }
1169        else
1170        {
1171          if(printlevel>=1) {"N(B) = N("+string(i)+") = "+string(N(i))+" suffices the conditions of the underlying theorem,";
1172                             "since P1 = (0:1:0) and P2[3] in (Z/N("+string(i)+")Z)*.";
1173                             "In particular N = "+string(N)+" is prime.";pause();"";}
1174          return(1);
1175        }
1176      }
1177
1178      if(step==14)                 // (14)[Backtrack]
1179      {
1180        if(i>0)
1181        {
1182          if(printlevel>=1) {"Set i = i-1 and restart the algorithm for N("+string(i-1)+") = "+string(N(i-1))+" with";
1183                             "a new discriminant.";pause();"";}
1184          i=i-1;
1185          k=0;
1186          step=3;
1187        }
1188        else
1189        {
1190          if(printlevel>=1) {"N(0) = N = "+string(N)+" and therefore N is not prime.";pause();"";}
1191          return(-1);
1192        }
1193      }
1194    }
1195  }
1196}
1197example
1198{ "EXAMPLE:"; echo = 2;
1199    ring R = 0,x,dp;
1200    Atkin(7691,100,5);
1201    Atkin(3473,10,2);
1202    printlevel=1;
1203    Atkin(10000079,100,2);
1204}
1205
Note: See TracBrowser for help on using the repository browser.