# source:git/Singular/LIB/AtkinsTest.lib@53e03a6

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