source: git/Singular/LIB/atkins.lib @ 05a2f6

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