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

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