source: git/Singular/LIB/atkins.lib @ 3eadab

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