source: git/Singular/LIB/AtkinsTest.lib @ 0c0b9f1

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