# source:git/Singular/LIB/atkins.lib@0dd77c2

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