# source:git/Singular/LIB/atkins.lib@7d56875

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