# source:git/Singular/LIB/crypto.lib@f076709

spielwiese
Last change on this file since f076709 was f076709, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix gcdN, lcmN git-svn-id: file:///usr/local/Singular/svn/trunk@14306 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 52.5 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id\$";
3category="Teaching";
4info="
5LIBRARY:  crypto.lib     Procedures for teaching cryptography
6AUTHOR:                  Gerhard Pfister, pfister@mathematik.uni-kl.de
7
8OVERVIEW:
9     The library contains procedures to compute the discrete logarithm,
10      primality-tests, factorization included elliptic curves.
11      The library is intended to be used for teaching purposes but not
12      for serious computations. Sufficiently high printlevel allows to
13      control each step, thus illustrating the algorithms at work.
14
15
16PROCEDURES:
17 decimal(s);                number corresponding to the hexadecimal number s
18 exgcdN(a,n)                compute s,t,d such that d=gcd(a,n)=s*a+t*n
19 eexgcdN(L)                 T with sum_i L[i]*T[i]=T[n+1]=gcd(L[1],...,L[n])
20 gcdN(a,b)                  compute gcd(a,b)
21 lcmN(a,b)                  compute lcm(a,b)
22 powerN(m,d,n)              compute m^d mod n
23 chineseRem(T,L)            compute x such that x = T[i] mod L[i]
24 Jacobi(a,n)                the generalized Legendre symbol of a and n
25 primList(n)                the list of all primes <=n
26 primL(q)                   first primes p_1,...,p_r such that q<p_1*...*p_r
27 intPart(x)                 the integral part of a rational number
28 intRoot(m)                 the integral part of the square root of m
29 squareRoot(a,p)            the square root of a in Z/p, p prime
30 solutionsMod2(M)           basis solutions of Mx=0 over Z/2
31 powerX(q,i,I)              q-th power of the i-th variable modulo I
32 babyGiant(b,y,p)           discrete logarithm x: b^x=y mod p
33 rho(b,y,p)                 discrete logarithm x: b^x=y mod p
34 MillerRabin(n,k)           probabilistic primaly-test of Miller-Rabin
35 SolowayStrassen(n,k)       probabilistic primaly-test of Soloway-Strassen
36 PocklingtonLehmer(N,[])    primaly-test of Pocklington-Lehmer
37 PollardRho(n,k,a,[])       Pollard's rho factorization
38 pFactor(n,B,P)             Pollard's p-factorization
40 isOnCurve(N,a,b,P)         P is on the curve y^2z=x^3+a*xz^2+b*z^3  over Z/N
42 ellipticMult(N,a,b,P,k)    k*P on elliptic curves
43 ellipticRandomCurve(N)     generates y^2z=x^3+a*xz^2+b*z^3  over Z/N randomly
44 ellipticRandomPoint(N,a,b) random point on y^2z=x^3+a*xz^2+b*z^3  over Z/N
45 countPoints(N,a,b)         number of points of y^2=x^3+a*x+b  over Z/N
46 ellipticAllPoints(N,a,b)   points of y^2=x^3+a*x+b  over Z/N
47 ShanksMestre(q,a,b,[])     number of points of y^2=x^3+a*x+b  over Z/N
48 Schoof(N,a,b)              number of points of y^2=x^3+a*x+b  over Z/N
49 generateG(a,b,m)           m-th division polynomial of y^2=x^3+a*x+b  over Z/N
50 factorLenstraECM(N,S,B,[]) Lenstra's factorization
51 ECPP(N)                    primaly-test of Goldwasser-Kilian
52
53              [parameters in square brackets are optional]
54";
55
56LIB "poly.lib";
57
58///////////////////////////////////////////////////////////////////////////////
59
60
61//=============================================================================
62//=========================== basic prozedures ================================
63//=============================================================================
64
65proc decimal(string s)
66"USAGE:  decimal(s); s = string
67RETURN: the (decimal) number corresponding to the hexadecimal number s
68EXAMPLE:example decimal; shows an example
69"
70{
71   int n=size(s);
72   int i;
73   bigint k;
74   bigint t=16;
75   bigint m=0;
76   for(i=1;i<=n;i++)
77   {
78      if(s[i]=="1"){k=1;}
79      if(s[i]=="2"){k=2;}
80      if(s[i]=="3"){k=3;}
81      if(s[i]=="4"){k=4;}
82      if(s[i]=="5"){k=5;}
83      if(s[i]=="6"){k=6;}
84      if(s[i]=="7"){k=7;}
85      if(s[i]=="8"){k=8;}
86      if(s[i]=="9"){k=9;}
87      if(s[i]=="a"){k=10;}
88      if(s[i]=="b"){k=11;}
89      if(s[i]=="c"){k=12;}
90      if(s[i]=="d"){k=13;}
91      if(s[i]=="e"){k=14;}
92      if(s[i]=="f"){k=15;}
93      m=m*t+k;
94   }
95   return(m);
96}
97example
98{ "EXAMPLE:"; echo = 2;
99   string s  ="8edfe37dae96cfd2466d77d3884d4196";
100   decimal(s);
101}
102
103proc exgcdN(number a, number n)
104"USAGE:  exgcdN(a,n);
105RETURN: a list s,t,d of numbers satisfying d=gcd(a,n)=s*a+t*n
106EXAMPLE:example exgcdN; shows an example
107"
108{
109  number x=a mod n;
110  if(x==0){return(list(0,1,n))}
111  list l=exgcdN(n,x);
112  return(list(l[2],l[1]-(a-x)*l[2]/n,l[3]))
113}
114example
115{ "EXAMPLE:"; echo = 2;
116   ring R = 0,x,dp;
117   exgcdN(24,15);
118}
119
120proc eexgcdN(list L)
121"USAGE:  eexgcdN(L);
122RETURN: list T such that sum_i L[i]*T[i]=T[n+1]=gcd(L[1],...,L[n])
123EXAMPLE:example eexgcdN; shows an example
124"
125{
126   if(size(L)==2){return(exgcdN(L[1],L[2]));}
127   number p=L[size(L)];
128   L=delete(L,size(L));
129   list T=eexgcdN(L);
130   list S=exgcdN(T[size(T)],p);
131   int i;
132   for(i=1;i<=size(T)-1;i++)
133   {
134      T[i]=T[i]*S[1];
135   }
136   p=T[size(T)];
137   T[size(T)]=S[2];
138   T[size(T)+1]=S[3];
139   return(T);
140}
141example
142{ "EXAMPLE:"; echo = 2;
143   ring R = 0,x,dp;
144   eexgcdN(list(24,15,21));
145}
146
147proc gcdN(number a, number b)
148"USAGE:  gcdN(a,b);
149RETURN: gcd(a,b)
150EXAMPLE:example gcdN; shows an example
151"
152{
153   //if((a mod b)==0){return(b)}
154   //return(gcdN(b,a mod b));
155   return(gcd(a,b));
156}
157example
158{ "EXAMPLE:"; echo = 2;
159   ring R = 0,x,dp;
160   gcdN(24,15);
161}
162
163proc lcmN(number a, number b)
164"USAGE:  lcmN(a,b);
165RETURN: lcm(a,b);
166EXAMPLE:example lcmN; shows an example
167"
168{
169   //number d=gcdN(a,b);
170   //return(a*b/d);
171   return (a*b/gcd(a,b));
172}
173example
174{ "EXAMPLE:"; echo = 2;
175   ring R = 0,x,dp;
176   lcmN(24,15);
177}
178
179proc powerN(number m, number d, number n)
180"USAGE:  powerN(m,d,n);
181RETURN: m^d mod n
182EXAMPLE:example powerN; shows an example
183"
184{
185   if(d==0){return(1)}
186   int i;
187   if(n==0)
188   {
189      for(i=12;i>=2;i--)
190      {
191         if((d mod i)==0){return(powerN(m,d/i,n)^i);}
192      }
193      return(m*powerN(m,d-1,n));
194   }
195   for(i=12;i>=2;i--)
196   {
197      if((d mod i)==0){return(powerN(m,d/i,n)^i mod n);}
198   }
199   return(m*powerN(m,d-1,n) mod n);
200}
201example
202{ "EXAMPLE:"; echo = 2;
203   ring R = 0,x,dp;
204   powerN(24,15,7);
205}
206
207proc chineseRem(list T,list L)
208"USAGE:  chineseRem(T,L);
209RETURN: x such that x = T[i] mod L[i]
210NOTE:   chinese remainder theorem
211EXAMPLE:example chineseRem; shows an example
212"
213{
214   int i;
215   int n=size(L);
216   number N=1;
217   for(i=1;i<=n;i++)
218   {
219      N=N*L[i];
220   }
221   list M;
222   for(i=1;i<=n;i++)
223   {
224      M[i]=N/L[i];
225   }
226   list S=eexgcdN(M);
227   number x;
228   for(i=1;i<=n;i++)
229   {
230      x=x+S[i]*M[i]*T[i];
231   }
232   x=x mod N;
233   return(x);
234}
235example
236{ "EXAMPLE:"; echo = 2;
237   ring R = 0,x,dp;
238   chineseRem(list(24,15,7),list(2,3,5));
239}
240
241proc Jacobi(number a, number n)
242"USAGE:  Jacobi(a,n);
243RETURN: the generalized Legendre symbol
244NOTE: if n is an odd prime then Jacobi(a,n)=0,1,-1 if n|a, a=x^2 mod n,else
245EXAMPLE:example Jacobi; shows an example
246"
247{
248   int i;
249   int z=1;
250   number t=1;
251   number k;
252
253   if((((n-1)/2) mod 2)!=0){z=-1;}
254   if(a<0){return(z*Jacobi(-a,n));}
255   a=a mod n;
256   if(n==1){return(1);}
257   if(a==0){return(0);}
258
259   while(a!=0)
260   {
261      while((a mod 2)==0)
262      {
263         a=a/2;
264         if(((n mod 8)==3)||((n mod 8)==5)){t=-t;}
265      }
266      k=a;a=n;n=k;
267      if(((a mod 4)==3)&&((n mod 4)==3)){t=-t;}
268      a=a mod n;
269   }
270   if (n==1){return(t);}
271   return(0);
272}
273example
274{ "EXAMPLE:"; echo = 2;
275   ring R = 0,x,dp;
276   Jacobi(13580555397810650806,5792543);
277}
278
279proc primList(int n)
280"USAGE:  primList(n);
281RETURN: the list of all primes <=n
282EXAMPLE:example primList; shows an example
283"
284{
285   int i,j;
286   list re;
287   re[1]=2;
288   re[2]=3;
289   for(i=5;i<=n;i=i+2)
290   {
291     j=1;
292     while(j<=size(re))
293     {
294        if((i mod re[j])==0){break;}
295        j++;
296     }
297     if(j==size(re)+1){re[size(re)+1]=i;}
298   }
299   return(re);
300}
301example
302{ "EXAMPLE:"; echo = 2;
303   list L=primList(100);
304   size(L);
305   L[size(L)];
306}
307
308proc primL(number q)
309"USAGE:  primL(q);
310RETURN: list of the first primes p_1,...,p_r such that q>p_1*...*p_(r-1)
311        and q<p_1*...*p_r
312EXAMPLE:example primL; shows an example
313"
314{
315   int i,j;
316   list re;
317   re[1]=2;
318   re[2]=3;
319   number s=6;
320   i=3;
321   while(s<=q)
322   {
323     i=i+2;
324     j=1;
325     while(j<=size(re))
326     {
327        if((i mod re[j])==0){break;}
328        j++;
329     }
330     if(j==size(re)+1)
331     {
332        re[size(re)+1]=i;
333        s=s*i;
334     }
335   }
336   return(re);
337}
338example
339{ "EXAMPLE:"; echo = 2;
340   ring R = 0,x,dp;
341   primL(20);
342}
343
344proc intPart(number x)
345"USAGE:  intPart(x);
346RETURN: the integral part of a rational number
347EXAMPLE:example intPart; shows an example
348"
349{
350   return((numerator(x)-(numerator(x) mod denominator(x)))/denominator(x));
351}
352example
353{ "EXAMPLE:"; echo = 2;
354   ring R = 0,x,dp;
355   intPart(7/3);
356}
357
358proc intRoot(number m)
359"USAGE:  intRoot(m);
360RETURN: the integral part of the square root of m
361EXAMPLE:example intRoot; shows an example
362"
363{
364   number x=1;
365   number t=x^2;
366   number s=(x+1)^2;
367   while(((m>t)&&(m>s))||((m<t)&&(m<s)))
368   {
369      x=intPart(x/2+m/(2*x));  //Newton step
370      t=x^2;
371      if(t>m)
372      {
373         s=(x-1)^2;
374      }
375      else
376      {
377         s=(x+1)^2;
378      }
379   }
380   if(t>m){return(x-1);}
381   if(s==m){return(x+1);}
382   return(x);
383}
384example
385{ "EXAMPLE:"; echo = 2;
386   ring R = 0,x,dp;
387   intRoot(20);
388}
389
390proc squareRoot(number a, number p)
391"USAGE:  squareRoot(a,p);
392RETURN: the square root of a in Z/p, p prime
393NOTE:   assumes the Jacobi symbol is 1 or p=2.
394EXAMPLE:example squareRoot; shows an example
395"
396{
397   if(p==2){return(a);}
398   if((a mod p)==0){return(0);}
399   if(powerN(a,p-1,p)!=1)
400   {
401      "p is not prime";
402      return(number(-5));
403   }
404   number n=random(1,2147483647) mod p;
405   if(n==0){n=n+1;}
406   number j=Jacobi(n,p);
407   if(j==0)
408   {
409      "p is not prime";
410      return(number(-5));
411   }
412   if(j==1)
413   {
414      return(squareRoot(a,p));
415   }
416   number q=p-1;
417   number e=0;
418   number two=2;
419   number z,m,t;
420   while((q mod 2)==0)
421   {
422      e=e+1;
423      q=q/2;
424   }
425   number y=powerN(n,q,p);
426   number r=e;
427   number x=powerN(a,(q-1)/2,p);
428   number b=a*x^2 mod p;
429   x=a*x mod p;
430
431   while(((b-1) mod p)!=0)
432   {
433      m=0;z=b;
434      while(((z-1) mod p)!=0)
435      {
436         z=z^2 mod p;
437         m=m+1;
438      }
439      t=powerN(y,powerN(two,r-m-1,p),p);
440      y=t^2 mod p;
441      r=m;
442      x=x*t mod p;
443      b=b*y mod p;
444   }
445   return(x);
446}
447example
448{ "EXAMPLE:"; echo = 2;
449   ring R = 0,x,dp;
450   squareRoot(8315890421938608,32003);
451}
452
453
454proc solutionsMod2(matrix M)
455"USAGE:  solutionsMod2(M);
456RETURN: an intmat containing a basis of the vector space of solutions of the
457        linear system of equations defined by M over the prime field of
458        characteristic 2
459EXAMPLE:example solutionsMod2; shows an example
460"
461{
462   def R=basering;
463   ring Rhelp=2,z,(c,dp);
464   matrix M=imap(R,M);
465   matrix S=syz(M);
466   setring(R);
467   matrix S=imap(Rhelp,S);
468   int i,j;
469   intmat v[nrows(S)][ncols(S)];
470   for(i=1;i<=nrows(S);i++)
471   {
472      for(j=1;j<=ncols(S);j++)
473      {
474         if(S[i,j]==1){v[i,j]=1;}
475      }
476   }
477   return(v);
478}
479example
480{ "EXAMPLE:"; echo = 2;
481   ring R = 0,x,dp;
482   matrix M[3][3]=1,2,3,4,5,6,7,6,5;
483   solutionsMod2(M);
484}
485
486proc powerX(int q, int i, ideal I)
487"USAGE:  powerX(q,i,I);
488RETURN: the q-th power of the i-th variable modulo I
489ASSUME: I is a standard basis
490EXAMPLE:example powerX; shows an example
491"
492{
493   if(q<=181){return(reduce(var(i)^int(q),I));}
494   if((q mod 5)==0){return(reduce(powerX(q div 5,i,I)^5,I));}
495   if((q mod 4)==0){return(reduce(powerX(q div 4,i,I)^4,I));}
496   if((q mod 3)==0){return(reduce(powerX(q div 3,i,I)^3,I));}
497   if((q mod 2)==0){return(reduce(powerX(q div 2,i,I)^2,I));}
498   return(reduce(var(i)*powerX(q-1,i,I),I));
499}
500example
501{ "EXAMPLE:"; echo = 2;
502   ring R = 0,(x,y),dp;
503   powerX(100,2,std(ideal(x3-1,y2-x)));
504}
505
506//======================================================================
507//=========================== Discrete Logarithm =======================
508//======================================================================
509
510//============== Shank's baby step - giant step ========================
511
512proc babyGiant(number b, number y, number p)
513"USAGE:  babyGiant(b,y,p);
514RETURN: the discrete logarithm x: b^x=y mod p
515NOTE:   This procedure works based on Shank's baby step - giant step method.
516EXAMPLE:example babyGiant; shows an example
517"
518{
519   int i,j,m;
520   list l;
521   number h=1;
522   number x;
523
524//choose m minimal such that m^2>p
525   for(i=1;i<=p;i++){if(i^2>p) break;}
526   m=i;
527
528//baby-step:  compute the table y*b^i for 1<=i<=m
529   for(i=1;i<=m;i++){l[i]=y*b^i mod p;}
530
531//giant-step: compute b^(m+j), 1<=j<=m and search in the baby-step table
532//for an i with y*b^i=b^(m*j). If found then x=m*j-i
533   number g=b^m mod p;
534   while(j<m)
535   {
536      j++;
537      h=h*g mod p;
538      for(i=1;i<=m;i++)
539      {
540         if(h==l[i])
541         {
542            x=m*j-i;
543            j=m;
544            break;
545         }
546     }
547   }
548   return(x);
549}
550example
551{ "EXAMPLE:"; echo = 2;
552   ring R = 0,z,dp;
553   number b=2;
554   number y=10;
555   number p=101;
556   babyGiant(b,y,p);
557}
558
559//==============  Pollards rho  =================================
560
561proc rho(number b, number y, number p)
562"USAGE:  rho(b,y,p);
563RETURN: the discrete logarithm x=log_b(y): b^x=y mod p
564NOTE: Pollard's rho:
565       choose random f_0 in 0,...,p-2 ,e_0=0, define x_0=b^f_0, define
566       x_i=y^e_i*b^f_i as below. For i large enough there is i with
567       x_(i/2)=x_i. Let s:=e_(i/2)-e_i mod p-1 and t:=f_i-f_(i/2) mod p-1,
568       d=gcd(s,p-1)=u*s+v*(p-1) then x=tu/d +j*(p-1)/d for some j (to be
569       found by trying)
570EXAMPLE:example rho; shows an example
571"
572{
573   int i=1;
574   int j;
575   number s,t;
576   list e,f,x;
577
578   e[1]=0;
579   f[1]=random(0,2147483629) mod (p-1);
580   x[1]=powerN(b,f[1],p);
581   while(i)
582   {
583      if((x[i] mod 3)==1)
584      {
585         x[i+1]=y*x[i] mod p;
586         e[i+1]=e[i]+1 mod (p-1);
587         f[i+1]=f[i];
588      }
589      if((x[i] mod 3)==2)
590      {
591         x[i+1]=x[i]^2 mod p;
592         e[i+1]=e[i]*2 mod (p-1);
593         f[i+1]=f[i]*2 mod (p-1);
594      }
595      if((x[i] mod 3)==0)
596      {
597         x[i+1]=x[i]*b mod p;
598         e[i+1]=e[i];
599         f[i+1]=f[i]+1 mod (p-1);
600      }
601      i++;
602      for(j=i-1;j>=1;j--)
603      {
604         if(x[i]==x[j])
605         {
606            s=(e[j]-e[i]) mod (p-1);
607            t=(f[i]-f[j]) mod (p-1);
608            if(s!=0)
609            {
610               i=0;
611            }
612            else
613            {
614               e[1]=0;
615               f[1]=random(0,2147483629) mod (p-1);
616               x[1]=powerN(b,f[1],p);
617               i=1;
618            }
619            break;
620         }
621      }
622   }
623
624   list w=exgcdN(s,p-1);
625   number u=w[1];
626   number d=w[3];
627
628   number a=(t*u/d) mod (p-1);
629
630   while(powerN(b,a,p)!=y)
631   {
632      a=(a+(p-1)/d) mod (p-1);
633   }
634   return(a);
635}
636example
637{ "EXAMPLE:"; echo = 2;
638   ring R = 0,x,dp;
639   number b=2;
640   number y=10;
641   number p=101;
642   rho(b,y,p);
643}
644//====================================================================
645//====================== Primality Tests =============================
646//====================================================================
647
648//================================= Miller-Rabin =====================
649
650proc MillerRabin(number n, int k)
651"USAGE:  MillerRabin(n,k);
652RETURN: 1 if n is prime, 0 else
653NOTE: probabilistic test of Miller-Rabin with k loops to test if n is prime.
654       Using the theorem: If n is prime, n-1=2^s*r, r odd, then
655       powerN(a,r,n)=1 or powerN(a,r*2^i,n)=-1 for some i
656EXAMPLE:example MillerRabin; shows an example
657"
658{
659   if(n<0){n=-n;}
660   if((n==2)||(n==3)){return(1);}
661   if((n mod 2)==0){return(0);}
662
663   int i;
664   number a,b,j,r,s;
665   r=n-1;
666   s=0;
667   while((r mod 2)==0)
668   {
669      s=s+1;
670      r=r/2;
671   }
672   while(i<k)
673   {
674      i++;
675      a=random(2,2147483629) mod n; if(a==0){a=3;}
676      if(exgcdN(a,n)[3]!=1){return(0);}
677      b=powerN(a,r,n);
678      if(b!=1)
679      {
680         j=0;
681         while(j<s)
682         {
683           if(((b+1) mod n)==0) break;
684           b=powerN(b,2,n);
685           j=j+1;
686         }
687         if(j==s){return(0);}
688      }
689   }
690   return(1);
691}
692example
693{ "EXAMPLE:"; echo = 2;
694   ring R = 0,z,dp;
695   number x=2;
696   x=x^787-1;
697   MillerRabin(x,3);
698}
699
700//======================= Soloway-Strassen  ==========================
701
702proc SolowayStrassen(number n, int k)
703"USAGE:  SolowayStrassen(n,k);
704RETURN: 1 if n is prime, 0 else
705NOTE: probabilistic test of Soloway-Strassen with k loops to test if n is
706       prime using the theorem: If n is prime then
707       powerN(a,(n-1)/2,n)=Jacobi(a,n) mod n
708EXAMPLE:example SolowayStrassen; shows an example
709"
710{
711   if(n<0){n=-n;}
712   if((n==2)||(n==3)){return(1);}
713   if((n mod 2)==0){return(0);}
714
715   number a;
716   int i;
717   while(i<k)
718   {
719      i++;
720      a=random(2,2147483629) mod n; if(a==0){a=3;}
721      if(gcdN(a,n)!=1){return(0);}
722      if(powerN(a,(n-1)/2,n)!=(Jacobi(a,n) mod n)){return(0);}
723   }
724   return(1);
725}
726example
727{ "EXAMPLE:"; echo = 2;
728   ring R = 0,z,dp;
729   number h=10;
730   number p=h^100+267;
731   //p=h^100+43723;
732   //p=h^200+632347;
733   SolowayStrassen(h,3);
734}
735
736
737/*
738ring R=0,z,dp;
739number p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
740number q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
741number n=p*q;
742SolowayStrassen(n,3);
743*/
744
745//===================== Pocklington-Lehmer ==============================
746
747proc PocklingtonLehmer(number N, list #)
748"USAGE: PocklingtonLehmer(N); optional: PocklingtonLehmer(N,L);
749        L a list of the first k primes
750RETURN:message N is not prime or {A,{p},{a_p}} as certificate for N being prime
751NOTE:assumes that it is possible to factorize N-1=A*B such that gcd(A,B)=1,
752      the factorization of A is completely known and A^2>N .
753      N is prime if and only if for each prime factor p of A we can find
754      a_p such that a_p^(N-1)=1 mod N and gcd(a_p^((N-1)/p)-1,N)=1
755EXAMPLE:example PocklingtonLehmer; shows an example
756"
757{
758   number m=intRoot(N);
759   if(size(#)>0)
760   {
761      list S=PollardRho(N-1,10000,1,#);
762   }
763   else
764   {
765      list S=PollardRho(N-1,10000,1);
766   }
767   int i,j;
768   number A=1;
769   number p,a,g;
770   list PA;
771   list re;
772
773   while(i<size(S))
774   {
775      p=S[i+1];
776      A=A*p;
777      PA[i+1]=p;
778      if(A>m){break;}
779
780      while(1)
781      {
782        p=p*S[i+1];
783        if(((N-1) mod p)==0)
784        {
785           A=A*p;
786        }
787        else
788        {
789           break;
790        }
791      }
792      i++;
793   }
794   if(A<=m)
795   {
796     A=N/A;
797     PA=list(S[size(S)]);
798   }
799   for(i=1;i<=size(PA);i++)
800   {
801      a=1;
802      while(a<N-1)
803      {
804         a=a+1;
805         if(powerN(a,N-1,N)!=1){return("not prime");}
806         g=gcdN(powerN(a,(N-1)/PA[i],N),N);
807         if(g==1)
808         {
809           re[size(re)+1]=list(PA[i],a);
810           break;
811         }
812         if(g<N){"not prime";return(g);}
813      }
814   }
815   return(list(A,re));
816}
817example
818{ "EXAMPLE:"; echo = 2;
819   ring R = 0,z,dp;
820   number N=105554676553297;
821   PocklingtonLehmer(N);
822   list L=primList(1000);
823   PocklingtonLehmer(N,L);
824}
825
826//=======================================================================
827//======================= Factorization =================================
828//=======================================================================
829
830//======================= Pollards rho  =================================
831
832proc PollardRho(number n, int k, int allFactors, list #)
833"USAGE:  PollardRho(n,k,allFactors); optional: PollardRho(n,k,allFactors,L);
834         L a list of the first k primes
835RETURN: a list of factors of n (which could be just n),if allFactors=0@*
836        a list of all factors of n ,if allFactors=1
837NOTE: probabilistic rho-algorithm of Pollard to find a factor of n in k loops.
838      Creates a sequence x_i such that (x_i)^2=(x_2i)^2 mod n for some i,
839      computes gcd(x_i-x_2i,n) to find a divisor. To define the sequence
840      choose x,a and define x_n+1=x_n^2+a mod n, x_1=x.
841      If allFactors is 1, it tries to find recursively all prime factors
842      using the Soloway-Strassen test.
844EXAMPLE:example PollardRho; shows an example
845"
846{
847    int i,j;
848    list L=primList(100);
849    list re,se;
850    if(n<0){n=-n;}
851    if(n==1){return(re);}
852
853//this is optional: test whether a prime of the list # devides n
854    if(size(#)>0)
855    {
856       L=#;
857    }
858    for(i=1;i<=size(L);i++)
859    {
860       if((n mod L[i])==0)
861       {
862          re[size(re)+1]=L[i];
863          while((n mod L[i])==0)
864          {
865             n=n/L[i];
866          }
867       }
868       if(n==1){return(re);}
869    }
870    int e=size(re);
871//here the rho-algorithm starts
872    number a,d,x,y;
873    while(n>1)
874    {
875       a=random(2,2147483629);
876       x=random(2,2147483629);
877       y=x;
878       d=1;
879       i=0;
880       while(i<k)
881       {
882          i++;
883          x=powerN(x,2,n); x=(x+a) mod n;
884          y=powerN(y,2,n); y=(y+a) mod n;
885          y=powerN(y,2,n); y=(y+a) mod n;
886          d=gcdN(x-y,n);
887          if(d>1)
888          {
889             re[size(re)+1]=d;
890             while((n mod d)==0)
891             {
892               n=n/d;
893             }
894             break;
895          }
896          if(i==k)
897          {
898             re[size(re)+1]=n;
899             n=1;
900          }
901       }
902
903    }
904    if(allFactors)      //want to obtain all prime factors
905    {
906       i=e;
907       while(i<size(re))
908       {
909          i++;
910
911          if(!SolowayStrassen(re[i],5))
912          {
913             se=PollardRho(re[i],2*k,1);
914             re[i]=se[size(se)];
915             for(j=1;j<=size(se)-1;j++)
916             {
917                re[size(re)+1]=se[j];
918             }
919             i--;
920          }
921       }
922    }
923    return(re);
924}
925example
926{ "EXAMPLE:"; echo = 2;
927   ring R = 0,z,dp;
928   number h=10;
929   number p=h^30+4;
930   PollardRho(p,5000,0);
931}
932
933//======================== Pollards p-factorization ================
934proc pFactor(number n,int B, list P)
935"USAGE:  pFactor(n,B.P); n to be factorized, B a bound , P a list of primes
936RETURN: a list of factors of n or n if no factor found
937NOTE: Pollard's p-factorization
938       creates the product k of powers of primes (bounded by B)  from
939       the list P with the idea that for a prime divisor p of n we have
940       p-1|k, and then p devides gcd(a^k-1,n) for some random a
941EXAMPLE:example pFactor; shows an example
942"
943{
944   int i;
945   number k=1;
946   number w;
947   while(i<size(P))
948   {
949      i++;
950      w=P[i];
951      if(w>B) break;
952      while(w*P[i]<=B)
953      {
954         w=w*P[i];
955      }
956      k=k*w;
957   }
958   number a=random(2,2147483629);
959   number d=gcdN(powerN(a,k,n)-1,n);
960   if((d>1)&&(d<n)){return(d);}
961   return(n);
962}
963example
964{ "EXAMPLE:"; echo = 2;
965   ring R = 0,z,dp;
966   list L=primList(1000);
967   pFactor(1241143,13,L);
968   number h=10;
969   h=h^30+25;
970   pFactor(h,20,L);
971}
972
974
975proc quadraticSieve(number n, int c, list B, int k)
976"USAGE:  quadraticSieve(n,c,B,k); n to be factorized, [-c,c] the
977         sieve-intervall, B a list of primes,
978         k for using the first k elements in B
979RETURN: a list of factors of n or the message: no divisor found
980NOTE: The idea being used is to find x,y such that x^2=y^2 mod n then
981      gcd(x-y,n) can be a proper divisor of n
983"
984{
985   number f,d;
986   int i,j,l,s,p;
987   list S,tmp;
988   intvec v;
989   v[k]=0;
990
991//compute the integral part of the square root of n
992   number m=intRoot(n);
993
994//consider the function f(X)=(X+m)^2-n and compute for s in [-c,c] the values
995   while(i<=2*c)
996   {
997      f=(i-c+m)^2-n;
998      tmp[1]=i-c+m;
999      tmp[2]=f;
1000      tmp[3]=v;
1001      S[i+1]=tmp;
1002      i++;
1003   }
1004
1005//the sieve with p in B
1006//find all s in [-c,c] such that f(s) has all prime divisors in the first
1007//k elements of B and the decomposition of f(s). They are characterized
1008//by 1 or -1 at the second place of S[j]:
1009//S[j]=j-c+m,f(j-c)/p_1^v_1*...*p_k^v_k, v_1,...,v_k maximal
1010   for(i=1;i<=k;i++)
1011   {
1012      p=B[i];
1013      if((p>2)&&(Jacobi(n,p)==-1)){i++;continue;}//n is no quadratic rest mod p
1014      j=1;
1015      while(j<=p)
1016      {
1017         if(j>2*c+1) break;
1018         f=S[j][2];
1019         v=S[j][3];
1020         s=0;
1021         while((f mod p)==0)
1022         {
1023            s++;
1024            f=f/p;
1025         }
1026         if(s)
1027         {
1028           S[j][2]=f;
1029           v[i]=s;
1030           S[j][3]=v;
1031           l=j;
1032           while(l+p<=2*c+1)
1033           {
1034              l=l+p;
1035              f=S[l][2];
1036              v=S[l][3];
1037              s=0;
1038              while((f mod p)==0)
1039              {
1040                s++;
1041                f=f/p;
1042              }
1043              S[l][2]=f;
1044              v[i]=s;
1045              S[l][3]=v;
1046           }
1047         }
1048         j++;
1049      }
1050   }
1051   list T;
1052   for(j=1;j<=2*c+1;j++)
1053   {
1054      if((S[j][2]==1)||(S[j][2]==-1))
1055      {
1056         T[size(T)+1]=S[j];
1057      }
1058   }
1059
1060//the system of equations for the exponents {l_s} for the f(s) such
1061//product f(s)^l_s is a square (l_s are 1 or 0)
1062   matrix M[k+1][size(T)];
1063   for(j=1;j<=size(T);j++)
1064   {
1065      if(T[j][2]==-1){M[1,j]=1;}
1066      for(i=1;i<=k;i++)
1067      {
1068         M[i+1,j]=T[j][3][i];
1069      }
1070   }
1071   intmat G=solutionsMod2(M);
1072
1073//construction of x and y such that x^2=y^2 mod n and d=gcd(x-y,n)
1074//y=square root of product f(s)^l_s
1075//x=product s+m
1076   number x=1;
1077   number y=1;
1078
1079   for(i=1;i<=ncols(G);i++)
1080   {
1081      kill v;
1082      intvec v;
1083      v[k]=0;
1084      for(j=1;j<=size(T);j++)
1085      {
1086         x=x*T[j][1]^G[j,i] mod n;
1087         if((T[j][2]==-1)&&(G[j,i]==1)){y=-y;}
1088         v=v+G[j,i]*T[j][3];
1089
1090      }
1091      for(l=1;l<=k;l++)
1092      {
1093         y=y*B[l]^(v[l] div 2) mod n;
1094      }
1095      d=gcdN(x-y,n);
1096      if((d>1)&&(d<n)){return(d);}
1097   }
1098   return("no divisor found");
1099}
1100example
1101{ "EXAMPLE:"; echo = 2;
1102   ring R = 0,z,dp;
1103   list L=primList(5000);
1106}
1107
1108//======================================================================
1109//==================== elliptic curves  ================================
1110//======================================================================
1111
1112//================= elementary operations ==============================
1113
1114proc isOnCurve(number N, number a, number b, list P)
1115"USAGE:  isOnCurve(N,a,b,P);
1116RETURN: 1 or 0 (depending on whether P is on the curve or not)
1117NOTE: checks whether P=(P[1]:P[2]:P[3]) is a point on the elliptic
1118      curve defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1119EXAMPLE:example isOnCurve; shows an example
1120"
1121{
1122   if(((P[2]^2*P[3]-P[1]^3-a*P[1]*P[3]^2-b*P[3]^3) mod N)!=0){return(0);}
1123   return(1);
1124}
1125example
1126{ "EXAMPLE:"; echo = 2;
1127   ring R = 0,z,dp;
1128   isOnCurve(32003,5,7,list(10,16,1));
1129}
1130
1131proc ellipticAdd(number N, number a, number b, list P, list Q)
1133RETURN: list L, representing the point P+Q
1134NOTE: P=(P[1]:P[2]:P[3]), Q=(Q[1]:Q[2]:Q[3]) points on the elliptic curve
1135      defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1137"
1138{
1139   if(N==2){ERROR("not implemented for 2");}
1140   int i;
1141   for(i=1;i<=3;i++)
1142   {
1143      P[i]=P[i] mod N;
1144      Q[i]=Q[i] mod N;
1145   }
1146   list Resu;
1147   Resu[1]=number(0);
1148   Resu[2]=number(1);
1149   Resu[3]=number(0);
1150   list Error;
1151   Error[1]=0;
1152   //test for ellictic curve
1153   number D=4*a^3+27*b^2;
1154   number g=gcdN(D,N);
1155   if(g==N){return(Error);}
1156   if(g!=1)
1157   {
1158      P[4]=g;
1159      return(P);
1160   }
1161   if(((P[1]==0)&&(P[2]==0)&&(P[3]==0))||((Q[1]==0)&&(Q[2]==0)&&(Q[3]==0)))
1162   {
1163      Error[1]=-2;
1164      return(Error);
1165   }
1166   if(!isOnCurve(N,a,b,P)||!isOnCurve(N,a,b,Q))
1167   {
1168      Error[1]=-1;
1169      return(Error);
1170   }
1171   if(P[3]==0){return(Q);}
1172   if(Q[3]==0){return(P);}
1173   list I=exgcdN(P[3],N);
1174   if(I[3]!=1)
1175   {
1176      P[4]=I[3];
1177      return(P);
1178   }
1179   P[1]=P[1]*I[1] mod N;
1180   P[2]=P[2]*I[1] mod N;
1181   I=exgcdN(Q[3],N);
1182   if(I[3]!=1)
1183   {
1184      P[4]=I[3];
1185      return(P);
1186   }
1187   Q[1]=Q[1]*I[1] mod N;
1188   Q[2]=Q[2]*I[1] mod N;
1189   if((P[1]==Q[1])&&(((P[2]+Q[2]) mod N)==0)){return(Resu);}
1190   number L;
1191   if((P[1]==Q[1])&&(P[2]==Q[2]))
1192   {
1193      I=exgcdN(2*Q[2],N);
1194      if(I[3]!=1)
1195      {
1196         P[4]=I[3];
1197         return(P);
1198      }
1199      L=I[1]*(3*Q[1]^2+a) mod N;
1200   }
1201   else
1202   {
1203      I=exgcdN(Q[1]-P[1],N);
1204      if(I[3]!=1)
1205      {
1206         P[4]=I[3];
1207         return(P);
1208      }
1209      L=(Q[2]-P[2])*I[1] mod N;
1210   }
1211   Resu[1]=(L^2-P[1]-Q[1]) mod N;
1212   Resu[2]=(L*(P[1]-Resu[1])-P[2]) mod N;
1213   Resu[3]=number(1);
1214   return(Resu);
1215}
1216example
1217{ "EXAMPLE:"; echo = 2;
1218   ring R = 0,z,dp;
1219   number N=11;
1220   number a=1;
1221   number b=6;
1222   list P,Q;
1223   P[1]=2;
1224   P[2]=4;
1225   P[3]=1;
1226   Q[1]=3;
1227   Q[2]=5;
1228   Q[3]=1;
1230}
1231
1232proc ellipticMult(number N, number a, number b, list P, number k)
1233"USAGE:  ellipticMult(N,a,b,P,k);
1234RETURN: a list L representing the point k*P
1235NOTE:  P=(P[1]:P[2]:P[3]) a point on the elliptic curve defined by
1236       y^2z=x^3+a*xz^2+b*z^3  over Z/N
1237EXAMPLE:example ellipticMult; shows an example
1238"
1239{
1240   if(P[3]==0){return(P);}
1241   list resu;
1242   resu[1]=number(0);
1243   resu[2]=number(1);
1244   resu[3]=number(0);
1245
1246   if(k==0){return(resu);}
1247   if(k==1){return(P);}
1249   if(k==-1)
1250   {
1251      resu=P;
1252      resu[2]=N-P[2];
1253      return(resu);
1254   }
1255   if(k<0)
1256   {
1257      resu=ellipticMult(N,a,b,P,-k);
1258      return(ellipticMult(N,a,b,resu,-1));
1259   }
1260   if((k mod 2)==0)
1261   {
1262      resu=ellipticMult(N,a,b,P,k/2);
1264   }
1265   resu=ellipticMult(N,a,b,P,k-1);
1267}
1268example
1269{ "EXAMPLE:"; echo = 2;
1270   ring R = 0,z,dp;
1271   number N=11;
1272   number a=1;
1273   number b=6;
1274   list P;
1275   P[1]=2;
1276   P[2]=4;
1277   P[3]=1;
1278   ellipticMult(N,a,b,P,3);
1279}
1280
1281//================== Random for elliptic curves =====================
1282
1283proc ellipticRandomCurve(number N)
1284"USAGE:  ellipticRandomCurve(N);
1285RETURN: a list of two random numbers a,b and 4a^3+27b^2 mod N
1286NOTE:   y^2z=x^3+a*xz^2+b^2*z^3 defines an elliptic curve over Z/N
1287EXAMPLE:example ellipticRandomCurve; shows an example
1288"
1289{
1290   int k;
1291   while(k<=10)
1292   {
1293     k++;
1294     number a=random(1,2147483647) mod N;
1295     number b=random(1,2147483647) mod N;
1296     //test for ellictic curve
1297     number D=4*a^3+27*b^4; //the constant term is b^2
1298     number g=gcdN(D,N);
1299     if(g<N){return(list(a,b,g));}
1300   }
1301   ERROR("no random curve found");
1302}
1303example
1304{ "EXAMPLE:"; echo = 2;
1305   ring R = 0,z,dp;
1306   ellipticRandomCurve(32003);
1307}
1308
1309proc ellipticRandomPoint(number N, number a, number b)
1310"USAGE:  ellipticRandomPoint(N,a,b);
1311RETURN: a list representing  a random point (x:y:z) of the elliptic curve
1312        defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1313EXAMPLE:example ellipticRandomPoint; shows an example
1314"
1315{
1316   number x=random(1,2147483647) mod N;
1317   number h=x^3+a*x+b;
1318   h=h mod N;
1319   list resu;
1320   resu[1]=x;
1321   resu[2]=0;
1322   resu[3]=1;
1323   if(h==0){return(resu);}
1324
1325   number n=Jacobi(h,N);
1326   if(n==0)
1327   {
1328      resu=-5;
1329      "N is not prime";
1330      return(resu);
1331   }
1332   if(n==1)
1333   {
1334      resu[2]=squareRoot(h,N);
1335      return(resu);
1336   }
1337   return(ellipticRandomPoint(N,a,b));
1338}
1339example
1340{ "EXAMPLE:"; echo = 2;
1341   ring R = 0,z,dp;
1342   ellipticRandomPoint(32003,3,181);
1343}
1344
1345
1346
1347//====================================================================
1348//======== counting the points of an elliptic curve  =================
1349//====================================================================
1350
1351//==================   the trivial approaches  =======================
1352proc countPoints(number N, number a, number b)
1353"USAGE:  countPoints(N,a,b);
1354RETURN: the number of points of the elliptic curve defined by
1355        y^2=x^3+a*x+b  over Z/N
1356NOTE: trivial approach
1357EXAMPLE:example countPoints; shows an example
1358"
1359{
1360  number x;
1361  number r=N+1;
1362  while(x<N)
1363  {
1364     r=r+Jacobi((x^3+a*x+b) mod N,N);
1365     x=x+1;
1366  }
1367  return(r);
1368}
1369example
1370{ "EXAMPLE:"; echo = 2;
1371   ring R = 0,z,dp;
1372   countPoints(181,71,150);
1373}
1374
1375proc ellipticAllPoints(number N, number a, number b)
1376"USAGE:  ellipticAllPoints(N,a,b);
1377RETURN: list of points (x:y:z) of the elliptic curve defined by
1378        y^2z=x^3+a*xz^2+b*z^3  over Z/N
1379EXAMPLE:example ellipticAllPoints; shows an example
1380"
1381{
1382   list resu,point;
1383   point[1]=0;
1384   point[2]=1;
1385   point[3]=0;
1386   resu[1]=point;
1387   point[3]=1;
1388   number x,h,n;
1389   while(x<N)
1390   {
1391      h=(x^3+a*x+b) mod N;
1392      if(h==0)
1393      {
1394         point[1]=x;
1395         point[2]=0;
1396         resu[size(resu)+1]=point;
1397      }
1398      else
1399      {
1400         n=Jacobi(h,N);
1401         if(n==1)
1402         {
1403            n=squareRoot(h,N);
1404            point[1]=x;
1405            point[2]=n;
1406            resu[size(resu)+1]=point;
1407            point[2]=N-n;
1408            resu[size(resu)+1]=point;
1409         }
1410      }
1411      x=x+1;
1412   }
1413   return(resu);
1414}
1415example
1416{ "EXAMPLE:"; echo = 2;
1417   ring R = 0,z,dp;
1418   list L=ellipticAllPoints(181,71,150);
1419   size(L);
1420   L[size(L)];
1421}
1422
1423//================ the algorithm of Shanks and Mestre =================
1424
1425proc ShanksMestre(number q, number a, number b, list #)
1426"USAGE:  ShanksMestre(q,a,b); optional: ShanksMestre(q,a,b,s); s the number
1427         of loops in the algorithm (default s=1)
1428RETURN: the number of points of the elliptic curve defined by
1429         y^2=x^3+a*x+b  over Z/N
1430NOTE: algorithm of Shanks and Mestre (baby-step-giant-step)
1431EXAMPLE:example ShanksMestre; shows an example
1432"
1433{
1434   number n=intRoot(4*q);
1435   number m=intRoot(intRoot(16*q))+1;
1436   number d;
1437   int i,j,k,s;
1438   list B,K,T,P,Q,R,mP;
1439   B[1]=list(0,1,0);
1440   if(size(#)>0)
1441   {
1442      s=#[1];
1443   }
1444   else
1445   {
1446      s=1;
1447   }
1448   while(k<s)
1449   {
1450      P =ellipticRandomPoint(q,a,b);
1451      Q =ellipticMult(q,a,b,P,n+q+1);
1452
1453      while(j<m)
1454      {
1455         j++;
1457      }
1459      mP[2]=q-mP[2];
1460      while(i<m)                            //giant-step
1461      {
1462         j=0;
1463         while(j<m)
1464         {
1465            j=j+1;
1466            if((Q[1]==B[j][1])&&(Q[2]==B[j][2])&&(Q[3]==B[j][3]))
1467            {
1468
1469               T[1]=P;
1470               T[2]=q+1+n-(i*m+j-1);
1471               K[size(K)+1]=T;
1472               if(size(K)>1)
1473               {
1474                  if(K[size(K)][2]!=K[size(K)-1][2])
1475                  {
1476                     d=gcdN(K[size(K)][2],K[size(K)-1][2]);
1477                     if(ellipticMult(q,a,b,K[size(K)],d)[3]==0)
1478                     {
1479                       K[size(K)][2]=K[size(K)-1][2];
1480                     }
1481                  }
1482               }
1483               i=int(m);
1484               break;
1485            }
1486         }
1487         i=i+1;
1489      }
1490      k++;
1491   }
1492   if(size(K)>0)
1493   {
1494      int te=1;
1495      for(i=1;i<=size(K)-1;i++)
1496      {
1497         if(K[size(K)][2]!=K[i][2])
1498         {
1499           if(ellipticMult(q,a,b,K[i],K[size(K)][2])[3]!=0)
1500           {
1501             te=0;
1502             break;
1503           }
1504         }
1505      }
1506      if(te)
1507      {
1508         return(K[size(K)][2]);
1509      }
1510   }
1511   return(ShanksMestre(q,a,b,s));
1512}
1513example
1514{ "EXAMPLE:"; echo = 2;
1515   ring R = 0,z,dp;
1516   ShanksMestre(32003,71,602);
1517}
1518
1519//==================== Schoof's algorithm =============================
1520
1521proc Schoof(number N,number a, number b)
1522"USAGE:  Schoof(N,a,b);
1523RETURN: the number of points of the elliptic curve defined by
1524        y^2=x^3+a*x+b  over Z/N
1525NOTE:  algorithm of Schoof
1526EXAMPLE:example Schoof; shows an example
1527"
1528{
1529   int pr=printlevel;
1530//test for ellictic curve
1531   number D=4*a^3+27*b^2;
1532   number G=gcdN(D,N);
1533   if(G==N){ERROR("not an elliptic curve");}
1534   if(G!=1){ERROR("not a prime");}
1535
1536//=== small N
1537  // if((N<=500)&&(pr<5)){return(countPoints(int(N),a,b));}
1538
1539//=== the general case
1540   number q=intRoot(4*N);
1541   list L=primL(2*q);
1542   int r=size(L);
1543   list T;
1544   int i,j;
1545   for(j=1;j<=r;j++)
1546   {
1547      T[j]=(testElliptic(int(N),a,b,L[j])+int(q)) mod L[j];
1548   }
1549   if(pr>=5)
1550   {
1551      "===================================================================";
1552      "Chinese remainder :";
1553      for(i=1;i<=size(T);i++)
1554      {
1555         " x =",T[i]," mod ",L[i];
1556      }
1557      "gives t+ integral part of the square root of q (to be positive)";
1558      chineseRem(T,L);
1559      "we obtain t = ",chineseRem(T,L)-q;
1560      "===================================================================";
1561   }
1562   number t=chineseRem(T,L)-q;
1563   return(N+1-t);
1564}
1565example
1566{ "EXAMPLE:"; echo = 2;
1567   ring R = 0,z,dp;
1568   Schoof(32003,71,602);
1569}
1570
1571/*
1572needs 518 sec
1573Schoof(2147483629,17,3567);
15742147168895
1575*/
1576
1577
1578proc generateG(number a,number b, int m)
1579"USAGE:  generateG(a,b,m);
1580RETURN: m-th division polynomial
1581NOTE: generate the so-called division polynomials, i.e., the recursively defined
1582polynomials p_m=generateG(a,b,m) in Z[x, y] such that, for a point (x:y:1) on the
1583elliptic curve defined by y^2=x^3+a*x+b  over Z/N the point@*
1584m*P=(x-(p_(m-1)*p_(m+1))/p_m^2 :(p_(m+2)*p_(m-1)^2-p_(m-2)*p_(m+1)^2)/4y*p_m^3 :1)
1585and m*P=0 if and only if p_m(P)=0
1586EXAMPLE:example generateG; shows an example
1587"
1588{
1589   if(m==0){return(poly(0));}
1590   if(m==1){return(poly(1));}
1591   if(m==2){return(2*var(1));}
1592   if(m==3){return(3*var(2)^4+6*a*var(2)^2+12*b*var(2)-a^2);}
1593   if(m==4)
1594   {
1595      return(4*var(1)*(var(2)^6+5*a*var(2)^4+20*b*var(2)^3-5*a^2*var(2)^2
1596        -4*a*b*var(2)-8*b^2-a^3));
1597   }
1598   if((m mod 2)==0)
1599   {
1600      return((generateG(a,b,m div 2+2)*generateG(a,b,m div 2-1)^2
1601        -generateG(a,b,m div 2-2)*generateG(a,b,m div 2+1)^2)
1602        *generateG(a,b,m div 2)/(2*var(1)));
1603   }
1604   return(generateG(a,b,(m-1) div 2+2)*generateG(a,b,(m-1) div 2)^3
1605     -generateG(a,b,(m-1) div 2-1)*generateG(a,b,(m-1) div 2+1)^3);
1606}
1607example
1608{ "EXAMPLE:"; echo = 2;
1609   ring R = 0,(x,y),dp;
1610   generateG(7,15,4);
1611}
1612
1613
1614proc testElliptic(int q,number a,number b,int l)
1615"USAGE:  testElliptic(q,a,b,l);
1616RETURN: an integer t, the trace of the Frobenius
1617NOTE: the kernel for the Schoof algorithm: looks for the t such that for all
1618      points (x:y:1) in C[l]={P in C | l*P=0},C the elliptic curve defined by
1619      y^2=x^3+a*x+b  over Z/q with group structure induced by 0=(0:1:0),
1620      (x:y:1)^(q^2)-t*(x:y:1)^q -ql*(x:y:1)=(0:1:0), ql= q mod l, trace of
1621      Frobenius.
1622EXAMPLE:example testElliptic; shows an example
1623"
1624{
1625   int pr=printlevel;
1626   def R=basering;
1627   ring S=q,(y,x),(L(100000),lp);
1628   number a=imap(R,a);
1629   number b=imap(R,b);
1630   poly F=y2-x3-a*x-b;       // the curve C
1631   poly G=generateG(a,b,l);
1632   ideal I=std(ideal(F,G));  // the points C[l]
1633   poly xq=powerX(q,2,I);
1634   poly yq=powerX(q,1,I);
1635   poly xq2=reduce(subst(xq,x,xq,y,yq),I);
1636   poly yq2=reduce(subst(yq,x,xq,y,yq),I);
1637   ideal J;
1638   int ql=q mod l;
1639   if(ql==0){ERROR("q is not prime");}
1640   int t;
1641   poly F1,F2,G1,G2,P1,P2,Q1,Q2,H1,H2,L1,L2;
1642
1643   if(pr>=5)
1644   {
1645      "===================================================================";
1646      "q=",q;
1647      "l=",l;
1648      "q mod l=",ql;
1649      "the Groebner basis for C[l]:";I;
1650      "x^q mod I = ",xq;
1651      "x^(q^2) mod I = ",xq2;
1652      "y^q mod I = ",yq;
1653      "y^(q^2) mod I = ",yq2;
1654      pause();
1655   }
1656   //==== l=2 =============================================================
1657   if(l==2)
1658   {
1659      xq=powerX(q,2,std(x3+a*x+b));
1660      J=std(ideal(xq-x,x3+a*x+b));
1661      if(deg(J[1])==0){t=1;}
1662      if(pr>=5)
1663      {
1664         "===================================================================";
1665         "the case l=2";
1666         "the gcd(x^q-x,x^3+ax+b)=",J[1];
1667         pause();
1668      }
1669      setring R;
1670      return(t);
1671   }
1672   //=== (F1/G1,F2/G2)=[ql](x,y) ==========================================
1673   if(ql==1)
1674   {
1675      F1=x;G1=1;F2=y;G2=1;
1676   }
1677   else
1678   {
1679      G1=reduce(generateG(a,b,ql)^2,I);
1680      F1=reduce(x*G1-generateG(a,b,ql-1)*generateG(a,b,ql+1),I);
1681      G2=reduce(4*y*generateG(a,b,ql)^3,I);
1682      F2=reduce(generateG(a,b,ql+2)*generateG(a,b,ql-1)^2
1683               -generateG(a,b,ql-2)*generateG(a,b,ql+1)^2,I);
1684
1685   }
1686   if(pr>=5)
1687   {
1688      "===================================================================";
1689      "the point ql*(x,y)=(F1/G1,F2/G2)";
1690      "F1=",F1;
1691      "G1=",G1;
1692      "F2=",F2;
1693      "G2=",G2;
1694      pause();
1695   }
1696   //==== the case t=0 :  the equations for (x,y)^(q^2)=-[ql](x,y) ===
1697   J[1]=xq2*G1-F1;
1698   J[2]=yq2*G2+F2;
1699   if(pr>=5)
1700   {
1701      "===================================================================";
1702      "the case t=0 mod l";
1703      "the equations for (x,y)^(q^2)=-[ql](x,y) :";
1704      J;
1705      "the test, if they vanish for all points in C[l]:";
1706      reduce(J,I);
1707      pause();
1708   }
1709   //=== test if all points of C[l] satisfy  (x,y)^(q^2)=-[ql](x,y)
1710   //=== if so: t mod l =0 is returned
1711   if(size(reduce(J,I))==0){setring R;return(0);}
1712
1713   //==== test for (x,y)^(q^2)=[ql](x,y) for some point
1714
1715   J=xq2*G1-F1,yq2*G2-F2;
1716   J=std(J+I);
1717   if(pr>=5)
1718   {
1719      "===================================================================";
1720      "test if (x,y)^(q^2)=[ql](x,y) for one point";
1721      "if so, the Frobenius has an eigenvalue 2ql/t: (x,y)^q=(2ql/t)*(x,y)";
1722      "it follows that t^2=4q mod l";
1723      "if w is one square root of q mod l";
1724      "t =2w mod l or -2w mod l ";
1725      "-------------------------------------------------------------------";
1726      "the equations for (x,y)^(q^2)=[ql](x,y) :";
1727      xq2*G1-F1,yq2*G2-F2;
1728      "the test if one point satisfies them";
1729      J;
1730      pause();
1731   }
1732   if(deg(J[1])>0)
1733   {
1734      setring R;
1735      int w=int(squareRoot(q,l));
1736      setring S;
1737      //=== +/-2w mod l zurueckgeben, wenn (x,y)^q=+/-[w](x,y)
1738      //==== the case t>0 :  (Q1/P1,Q2/P2)=[w](x,y) ==============
1739      if(w==1)
1740      {
1741         Q1=x;P1=1;Q2=y;P2=1;
1742      }
1743      else
1744      {
1745         P1=reduce(generateG(a,b,w)^2,I);
1746         Q1=reduce(x*G1-generateG(a,b,w-1)*generateG(a,b,w+1),I);
1747         P2=reduce(4*y*generateG(a,b,w)^3,I);
1748         Q2=reduce(generateG(a,b,w+2)*generateG(a,b,w-1)^2
1749               -generateG(a,b,w-2)*generateG(a,b,w+1)^2,I);
1750      }
1751      J=xq*P1-Q1,yq*P2-Q2;
1752      J=std(I+J);
1753      if(pr>=5)
1754      {
1755      "===================================================================";
1756      "the Frobenius has an eigenvalue, one of the roots of  w^2=q mod l:";
1757      "one root is:";w;
1758      "test, if it is the eigenvalue (if not it must be -w):";
1759      "the equations for (x,y)^q=w*(x,y)";I;xq*P1-Q1,yq*P2-Q2;
1760      "the Groebner basis";
1761       J;
1762         pause();
1763      }
1764      if(deg(J[1])>0){return(2*w mod l);}
1765      return(-2*w mod l);
1766   }
1767
1768   //==== the case t>0 :  (Q1/P1,Q2/P2)=(x,y)^(q^2)+[ql](x,y) =====
1769   P1=reduce(G1*G2^2*(F1-xq2*G1)^2,I);
1770   Q1=reduce((F2-yq2*G2)^2*G1^3-F1*G2^2*(F1-xq2*G1)^2-xq2*P1,I);
1771   P2=reduce(P1*G2*(F1-xq2*G1),I);
1772   Q2=reduce((xq2*P1-Q1)*(F2-yq2*G2)*G1-yq2*P2,I);
1773
1774   if(pr>=5)
1775   {
1776      "we are in the general case:";
1777      "(x,y)^(q^2)!=ql*(x,y) and (x,y)^(q^2)!=-ql*(x,y) ";
1778      "the point (Q1/P1,Q2/P2)=(x,y)^(q^2)+[ql](x,y)";
1779      "Q1=",Q1;
1780      "P1=",P1;
1781      "Q2=",Q2;
1782      "P2=",P2;
1783      pause();
1784   }
1785   while(t<(l-1) div 2)
1786   {
1787      t++;
1788      //====  (H1/L1,H2/L2)=[t](x,y)^q ===============================
1789      if(t==1)
1790      {
1791         H1=xq;L1=1;
1792         H2=yq;L2=1;
1793      }
1794      else
1795      {
1796         H1=x*generateG(a,b,t)^2-generateG(a,b,t-1)*generateG(a,b,t+1);
1797         H1=subst(H1,x,xq,y,yq);
1798         H1=reduce(H1,I);
1799         L1=generateG(a,b,t)^2;
1800         L1=subst(L1,x,xq,y,yq);
1801         L1=reduce(L1,I);
1802         H2=generateG(a,b,t+2)*generateG(a,b,t-1)^2
1803           -generateG(a,b,t-2)*generateG(a,b,t+1)^2;
1804         H2=subst(H2,x,xq,y,yq);
1805         H2=reduce(H2,I);
1806         L2=4*y*generateG(a,b,t)^3;
1807         L2=subst(L2,x,xq,y,yq);
1808         L2=reduce(L2,I);
1809      }
1810      J=Q1*L1-P1*H1,Q2*L2-P2*H2;
1811      if(pr>=5)
1812      {
1813      "we test now the different t, 0<t<=(l-1)/2:";
1814      "the point (H1/L1,H2/L2)=[t](x,y)^q :";
1815      "H1=",H1;
1816      "L1=",L1;
1817      "H2=",H2;
1818      "L2=",L2;
1819      "the equations for (x,y)^(q^2)+[ql](x,y)=[t](x,y)^q :";J;
1820      "the test";reduce(J,I);
1821      "the test for l-t (the x-cordinate is the same):";
1822       Q1*L1-P1*H1,Q2*L2+P2*H2;
1823       reduce(ideal(Q1*L1-P1*H1,Q2*L2+P2*H2),I);
1824         pause();
1825      }
1826      if(size(reduce(J,I))==0){setring R;return(t);}
1827      J=Q1*L1-P1*H1,Q2*L2+P2*H2;
1828      if(size(reduce(J,I))==0){setring R;return(l-t);}
1829   }
1830   ERROR("something is wrong in testElliptic");
1831}
1832example
1833{ "EXAMPLE:"; echo = 2;
1834   ring R = 0,z,dp;
1835   testElliptic(1267985441,338474977,64740730,3);
1836}
1837
1838//============================================================================
1839//================== Factorization and Primality Test ========================
1840//============================================================================
1841
1842//============= Lenstra's ECM Factorization ==================================
1843
1844proc factorLenstraECM(number N, list S, int B, list #)
1845"USAGE:  factorLenstraECM(N,S,B); optional: factorLenstraECM(N,S,B,d);
1846         d+1 the number of loops in the algorithm (default d=0)
1847RETURN: a factor of N or the message no factor found
1848NOTE: - computes a factor of N using Lenstra's ECM factorization@*
1849      - the idea is that the fact that N is not prime is dedected using
1850        the operations on the elliptic curve
1851      - is similarly to Pollard's p-1-factorization
1852EXAMPLE:example factorLenstraECM; shows an example
1853"
1854{
1855   list L,P;
1856   number g,M,w;
1857   int i,j,k,d;
1858   int l=size(S);
1859   if(size(#)>0)
1860   {
1861      d=#[1];
1862   }
1863
1864   while(i<=d)
1865   {
1866      L=ellipticRandomCurve(N);
1867      if(L[3]>1){return(L[3]);} //the discriminant was not invertible
1868      P=list(0,L[2],1);
1869      j=0;
1870      M=1;
1871      while(j<l)
1872      {
1873         j++;
1874         w=S[j];
1875         if(w>B) break;
1876         while(w*S[j]<B)
1877         {
1878           w=w*S[j];
1879         }
1880         M=M*w;
1881         P=ellipticMult(N,L[1],L[2]^2,P,w);
1882         if(size(P)==4){return(P[4]);}  //some inverse did not exsist
1883         if(P[3]==0){break;}            //the case M*P=0
1884      }
1885      i++;
1886   }
1887   return("no factor found");
1888}
1889example
1890{ "EXAMPLE:"; echo = 2;
1891   ring R = 0,z,dp;
1892   list L=primList(1000);
1893   factorLenstraECM(181*32003,L,10,5);
1894   number h=10;
1895   h=h^30+25;
1896   factorLenstraECM(h,L,4,3);
1897}
1898
1899//================= ECPP (Goldwasser-Kilian) a primaly-test =============
1900
1901proc ECPP(number N)
1902"USAGE:  ECPP(N);
1903RETURN: message:N is not prime or {L,P,m,q} as certificate for N being prime@*
1904         L a list (y^2=x^3+L[1]*x+L[2] defines an elliptic curve C)@*
1905         P a list ((P[1]:P[2]:P[3]) is a point of C)@*
1906         m,q integers
1907ASSUME: gcd(N,6)=1
1908NOTE:   The basis of the algorithm is the following theorem:
1909         Given C, an elliptic curve over Z/N, P a point of C(Z/N),
1910         m an integer, q a prime with the following properties:
1911         - q|m
1912         - q>(4-th root(N) +1)^2
1913         - m*P=0=(0:1:0)
1914         - (m/q)*P=(x:y:z) and z a unit in Z/N
1915         Then N is prime.
1916EXAMPLE:example ECPP; shows an example
1917"
1918{
1919   list L,S,P;
1920   number m,q;
1921   int i;
1922
1923   number n=intRoot(intRoot(N));
1924   n=(n+1)^2;                         //lower bound for q
1925   while(1)
1926   {
1927      L=ellipticRandomCurve(N);       //a random elliptic curve C
1928      m=ShanksMestre(N,L[1],L[2],3);  //number of points of the curve C
1929      S=PollardRho(m,10000,1);        //factorization of m
1930      for(i=1;i<=size(S);i++)         //search for q between the primes
1931      {
1932         q=S[i];
1933         if(n<q){break;}
1934      }
1935      if(n<q){break;}
1936   }
1937   number u=m/q;
1938   while(1)
1939   {
1940      P=ellipticRandomPoint(N,L[1],L[2]);  //a random point on C
1941      if(ellipticMult(N,L[1],L[2],P,m)[3]!=0){"N is not prime";return(-5);}
1942      if(ellipticMult(N,L[1],L[2],P,u)[3]!=0)
1943      {
1944         L=delete(L,3);
1945         return(list(L,P,m,q));
1946      }
1947   }
1948}
1949example
1950{ "EXAMPLE:"; echo = 2;
1951   ring R = 0,z,dp;
1952   number N=1267985441;
1953   ECPP(N);
1954}
1955
1956static proc wordToNumber(string s)
1957{
1958   int i;
1959   intvec v;
1960   number n;
1961   number t=27;
1962   for(i=size(s);i>0;i--)
1963   {
1964      if(s[i]=="a"){v[i]=0;}
1965      if(s[i]=="b"){v[i]=1;}
1966      if(s[i]=="c"){v[i]=2;}
1967      if(s[i]=="d"){v[i]=3;}
1968      if(s[i]=="e"){v[i]=4;}
1969      if(s[i]=="f"){v[i]=5;}
1970      if(s[i]=="g"){v[i]=6;}
1971      if(s[i]=="h"){v[i]=7;}
1972      if(s[i]=="i"){v[i]=8;}
1973      if(s[i]=="j"){v[i]=9;}
1974      if(s[i]=="k"){v[i]=10;}
1975      if(s[i]=="l"){v[i]=11;}
1976      if(s[i]=="m"){v[i]=12;}
1977      if(s[i]=="n"){v[i]=13;}
1978      if(s[i]=="o"){v[i]=14;}
1979      if(s[i]=="p"){v[i]=15;}
1980      if(s[i]=="q"){v[i]=16;}
1981      if(s[i]=="r"){v[i]=17;}
1982      if(s[i]=="s"){v[i]=18;}
1983      if(s[i]=="t"){v[i]=19;}
1984      if(s[i]=="u"){v[i]=20;}
1985      if(s[i]=="v"){v[i]=21;}
1986      if(s[i]=="w"){v[i]=22;}
1987      if(s[i]=="x"){v[i]=23;}
1988      if(s[i]=="y"){v[i]=24;}
1989      if(s[i]=="z"){v[i]=25;}
1990      if(s[i]==" "){v[i]=26;}
1991   }
1992   for(i=1;i<=size(s);i++)
1993   {
1994      n=n+v[i]*t^(i-1);
1995   }
1996   return(n);
1997}
1998
1999static proc numberToWord(number n)
2000{
2001   int i,j;
2002   string v;
2003   list s;
2004   number t=27;
2005   number mm;
2006   number nn=n;
2007   while(nn>t)
2008   {
2009      j++;
2010      mm=nn mod t;
2011      s[j]=mm;
2012      nn=(nn-mm)/t;
2013   }
2014   j++;
2015   s[j]=nn;
2016   for(i=1;i<=j;i++)
2017   {
2018      if(s[i]==0){v=v+"a";}
2019      if(s[i]==1){v=v+"b";}
2020      if(s[i]==2){v=v+"c";}
2021      if(s[i]==3){v=v+"d";}
2022      if(s[i]==4){v=v+"e";}
2023      if(s[i]==5){v=v+"f";}
2024      if(s[i]==6){v=v+"g";}
2025      if(s[i]==7){v=v+"h";}
2026      if(s[i]==8){v=v+"i";}
2027      if(s[i]==9){v=v+"j";}
2028      if(s[i]==10){v=v+"k";}
2029      if(s[i]==11){v=v+"l";}
2030      if(s[i]==12){v=v+"m";}
2031      if(s[i]==13){v=v+"n";}
2032      if(s[i]==14){v=v+"o";}
2033      if(s[i]==15){v=v+"p";}
2034      if(s[i]==16){v=v+"q";}
2035      if(s[i]==17){v=v+"r";}
2036      if(s[i]==18){v=v+"s";}
2037      if(s[i]==19){v=v+"t";}
2038      if(s[i]==20){v=v+"u";}
2039      if(s[i]==21){v=v+"v";}
2040      if(s[i]==22){v=v+"w";}
2041      if(s[i]==23){v=v+"x";}
2042      if(s[i]==24){v=v+"y";}
2043      if(s[i]==25){v=v+"z";}
2044      if(s[i]==26){v=v+" ";}
2045   }
2046   return(v);
2047}
2048
2049proc code(string s)
2050"USAGE:  code(s); s a string
2051ASSUME:  s contains only small letters and space
2052COMPUTE: a number, RSA-coding of the string s
2053RETURN:  return RSA-coding of the string s as string
2054EXAMPLE: code;  shows an example
2055"
2056{
2057   ring r=0,x,dp;
2058   number
2059p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
2060   number
2061q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
2062   number n=p*q;
2063   number phi=(p-1)*(q-1);
2064   number e=1234567891;
2065   list L=exgcdN(e,phi);
2066   number d=L[1];
2067   number m=wordToNumber(s);
2068   number c=powerN(m,e,n);
2069   string cc=string(c);
2070   return(cc);
2071}
2072example
2073{"EXAMPLE:";  echo = 2;
2074  string s="i go to school";
2075  code(s);
2076}
2077
2078proc decodeString(string g)
2079"USAGE:  decodeString(s); s a string
2080ASSUME:  s is a string of a number, the output of code
2081COMPUTE: a string, RSA-decoding of the string s
2082RETURN:  return RSA-decoding of the string s as string
2083EXAMPLE: decodeString;  shows an example
2084"
2085{
2086   ring r=0,x,dp;
2087   number
2088p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
2089   number
2090q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
2091   number n=p*q;
2092   number phi=(p-1)*(q-1);
2093   number e=1234567891;
2094   list L=exgcdN(e,phi);
2095   number d=L[1];
2096   execute("number c="+g+";");
2097   number f=powerN(c,d,n);
2098   string s=numberToWord(f);
2099   return(s);
2100}
2101example
2102{"EXAMPLE:";  echo = 2;
2103  string
2104s="78638618599886548153321853785991541374544958648147340831959482696082179852616053583234149080198937632782579537867262780982185252913122030800897193851413140758915381848932565";
2105  string t=decodeString(s);
2106  t;
2107}
2108
2109/*
2110//===============================================================
2111//=======  Example for DSA  =====================================
2112//===============================================================
2113Suppose a file test is given.It contains "Oscar".
2114
2115//Hash-function MD5 under Linux
2116
2117md5sum test           8edfe37dae96cfd2466d77d3884d4196
2118
2119//================================================================
2120
2121ring R=0,x,dp;
2122
2123number q=2^19+21;          //524309
2124number o=2*3*23*number(7883)*number(16170811);
2125
2126number p=o*q+1;            //9223372036869000547
2127number b=2;
2128number g=power(2,o,p);     //8308467587808723131
2129
2130number a=111111;
2131number A=power(g,a,p);     //8566038811843553785
2132
2133number h =decimal("8edfe37dae96cfd2466d77d3884d4196");
2134
2135                           //189912871665444375716340628395668619670
2136h= h mod q;                //259847
2137
2138number k=123456;
2139
2140number ki=exgcd(k,q)[1];    //50804
2141                            //inverse von k mod q
2142
2143number r= power(g,k,p) mod q;    //76646
2144
2145number s=ki*(h+a*r) mod q;       //2065
2146
2147//========== signatur is (r,s)=(76646,2065) =====================
2148//==================== verification  ============================
2149
2150number si=exgcd(s,q)[1];  //inverse von s mod q
2151number e1=si*h mod q;
2152number e2=si*r mod q;
2153number rr=((power(g,e1,p)*power(A,e2,p)) mod p) mod q;    //76646
2154
2155//===============================================================
2156//=======  Example for knapsack  ================================
2157//===============================================================
2158ring R=(5^5,t),x,dp;
2159R;
2160//   # ground field : 3125
2161//   primitive element : t
2162//   minpoly        : 1*t^5+4*t^1+2*t^0
2163//   number of vars : 1
2164//        block   1 : ordering dp
2165//                  : names    x
2166//        block   2 : ordering C
2167
2168proc findEx(number n, number g)
2169{
2170   int i;
2171   for(i=0;i<=size(basering)-1;i++)
2172   {
2173      if(g^i==n){return(i);}
2174   }
2175}
2176
2177number g=t^3;                       //choice of the primitive root
2178
2179findEx(t+1,g);
2180//2091
2181findEx(t+2,g);
2182//2291
2183findEx(t+3,g);
2184//1043
2185
2186intvec b=1,2091,2291,1043;          // k=4
2187int z=199;
2188intvec v=1043+z,1+z,2091+z,2291+z;  //permutation   pi=(0123)
2189v;
21901242,200,2290,2490
2191
2192//(1101)=(e_3,e_2,e_1,e_0)
2193//encoding 2490+2290+1242=6022 und 1+1+0+1=3
2194
2195//(6022,3) decoding: c-z*c'=6022-199*3=5425
2196
2197ring S=5,x,dp;
2198poly F=x5+4x+2;
2199poly G=reduce((x^3)^5425,std(F));
2200G;
2201//x3+x2+x+1
2202
2203factorize(G);
2204//[1]:
2205//   _[1]=1
2206//   _[2]=x+1
2207//   _[3]=x-2
2208//   _[4]=x+2
2209//[2]:
2210//   1,1,1,1
2211
2212//factors x+1,x+2,x+3, i.e.  (1110)=(e_pi(3),e_pi(2),e_pi(1),e_pi(0))
2213
2214//pi(0)=1,pi(1)=2,pi(2)=3,pi(3)=0 gives:  (1101)
2215
2216*/
2217
Note: See TracBrowser for help on using the repository browser.