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

spielwiese
Last change on this file since 7d56875 was 7d56875, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: typos reported by gorzelc git-svn-id: file:///usr/local/Singular/svn/trunk@11114 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 52.3 KB
Line
2///////////////////////////////////////////////////////////////////////////////
3version="\$Id: crypto.lib,v 1.6 2008-10-09 09:31:57 Singular Exp \$";
4category="Teaching";
5info="
6LIBRARY:  crypto.lib     Procedures for teaching cryptography
7AUTHOR:                  Gerhard Pfister, pfister@mathematik.uni-kl.de
8
9NOTE: The library contains procedures to compute the discrete logarithm,
10      primaly-tests, factorization included elliptic curve methodes.
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)                   all 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, 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}
156example
157{ "EXAMPLE:"; echo = 2;
158   ring R = 0,x,dp;
159   gcdN(24,15);
160}
161
162proc lcmN(number a, number b)
163"USAGE:  lcmN(a,b);
164RETURN: lcm(a,b);
165EXAMPLE:example lcmN; shows an example
166"
167{
168   number d=gcdN(a,b);
169   return(a*b/d);
170}
171example
172{ "EXAMPLE:"; echo = 2;
173   ring R = 0,x,dp;
174   lcmN(24,15);
175}
176
177proc powerN(number m, number d, number n)
178"USAGE:  powerN(m,d,n);
179RETURN: m^d mod n
180EXAMPLE:example powerN; shows an example
181"
182{
183   if(d==0){return(1)}
184   int i;
185   if(n==0)
186   {
187      for(i=12;i>=2;i--)
188      {
189         if((d mod i)==0){return(powerN(m,d/i,n)^i);}
190      }
191      return(m*powerN(m,d-1,n));
192   }
193   for(i=12;i>=2;i--)
194   {
195      if((d mod i)==0){return(powerN(m,d/i,n)^i mod n);}
196   }
197   return(m*powerN(m,d-1,n) mod n);
198}
199example
200{ "EXAMPLE:"; echo = 2;
201   ring R = 0,x,dp;
202   powerN(24,15,7);
203}
204
205proc chineseRem(list T,list L)
206"USAGE:  chineseRem(T,L);
207RETURN: x such that x = T[i] mod L[i]
208NOTE:   chinese remainder theorem
209EXAMPLE:example chineseRem; shows an example
210"
211{
212   int i;
213   int n=size(L);
214   number N=1;
215   for(i=1;i<=n;i++)
216   {
217      N=N*L[i];
218   }
219   list M;
220   for(i=1;i<=n;i++)
221   {
222      M[i]=N/L[i];
223   }
224   list S=eexgcdN(M);
225   number x;
226   for(i=1;i<=n;i++)
227   {
228      x=x+S[i]*M[i]*T[i];
229   }
230   x=x mod N;
231   return(x);
232}
233example
234{ "EXAMPLE:"; echo = 2;
235   ring R = 0,x,dp;
236   chineseRem(list(24,15,7),list(2,3,5));
237}
238
239proc Jacobi(number a, number n)
240"USAGE:  Jacobi(a,n);
241RETURN: the generalized Legendre symbol
242NOTE: if n is an odd prime then Jacobi(a,n)=0,1,-1 if n|a, a=x^2 mod n,else
243EXAMPLE:example Jacobi; shows an example
244"
245{
246   int i;
247   int z=1;
248   number t=1;
249   number k;
250
251   if((((n-1)/2) mod 2)!=0){z=-1;}
252   if(a<0){return(z*Jacobi(-a,n));}
253   a=a mod n;
254   if(n==1){return(1);}
255   if(a==0){return(0);}
256
257   while(a!=0)
258   {
259      while((a mod 2)==0)
260      {
261         a=a/2;
262         if(((n mod 8)==3)||((n mod 8)==5)){t=-t;}
263      }
264      k=a;a=n;n=k;
265      if(((a mod 4)==3)&&((n mod 4)==3)){t=-t;}
266      a=a mod n;
267   }
268   if (n==1){return(t);}
269   return(0);
270}
271example
272{ "EXAMPLE:"; echo = 2;
273   ring R = 0,x,dp;
274   Jacobi(13580555397810650806,5792543);
275}
276
277proc primList(int n)
278"USAGE:  primList(n);
279RETURN: the list of all primes <=n
280EXAMPLE:example primList; shows an example
281"
282{
283   int i,j;
284   list re;
285   re[1]=2;
286   re[2]=3;
287   for(i=5;i<=n;i=i+2)
288   {
289     j=1;
290     while(j<=size(re))
291     {
292        if((i mod re[j])==0){break;}
293        j++;
294     }
295     if(j==size(re)+1){re[size(re)+1]=i;}
296   }
297   return(re);
298}
299example
300{ "EXAMPLE:"; echo = 2;
301   list L=primList(100);
302   size(L);
303   L[size(L)];
304}
305
306proc primL(number q)
307"USAGE:  primL(q);
308RETURN: list of all primes p_1,...,p_r such that q<p_1*...*p_r
309EXAMPLE:example primL; shows an example
310"
311{
312   int i,j;
313   list re;
314   re[1]=2;
315   re[2]=3;
316   number s=6;
317   i=3;
318   while(s<=q)
319   {
320     i++;
321     j=1;
322     while(j<=size(re))
323     {
324        if((i mod re[j])==0){break;}
325        j++;
326     }
327     if(j==size(re)+1)
328     {
329        re[size(re)+1]=i;
330        s=s*i;
331     }
332   }
333   return(re);
334}
335example
336{ "EXAMPLE:"; echo = 2;
337   ring R = 0,x,dp;
338   primL(20);
339}
340
341proc intPart(number x)
342"USAGE:  intPart(x);
343RETURN: the integral part of a rational number
344EXAMPLE:example intPart; shows an example
345"
346{
347   return((numerator(x)-(numerator(x) mod denominator(x)))/denominator(x));
348}
349example
350{ "EXAMPLE:"; echo = 2;
351   ring R = 0,x,dp;
352   intPart(7/3);
353}
354
355proc intRoot(number m)
356"USAGE:  intRoot(m);
357RETURN: the integral part of the square root of m
358EXAMPLE:example intRoot; shows an example
359"
360{
361   number x=1;
362   number t=x^2;
363   number s=(x+1)^2;
364   while(((m>t)&&(m>s))||((m<t)&&(m<s)))
365   {
366      x=intPart(x/2+m/(2*x));  //Newton step
367      t=x^2;
368      if(t>m)
369      {
370         s=(x-1)^2;
371      }
372      else
373      {
374         s=(x+1)^2;
375      }
376   }
377   if(t>m){return(x-1);}
378   if(s==m){return(x+1);}
379   return(x);
380}
381example
382{ "EXAMPLE:"; echo = 2;
383   ring R = 0,x,dp;
384   intRoot(20);
385}
386
387proc squareRoot(number a, number p)
388"USAGE:  squareRoot(a,p);
389RETURN: the square root of a in Z/p, p prime
390NOTE:   assumes the Jacobi symbol is 1 or p=2.
391EXAMPLE:example squareRoot; shows an example
392"
393{
394   if(p==2){return(a);}
395   if((a mod p)==0){return(0);}
396   if(powerN(a,p-1,p)!=1)
397   {
398      "p is not prime";
399      return(number(-5));
400   }
401   number n=random(1,2147483647) mod p;
402   if(n==0){n=n+1;}
403   number j=Jacobi(n,p);
404   if(j==0)
405   {
406      "p is not prime";
407      return(number(-5));
408   }
409   if(j==1)
410   {
411      return(squareRoot(a,p));
412   }
413   number q=p-1;
414   number e=0;
415   number two=2;
416   number z,m,t;
417   while((q mod 2)==0)
418   {
419      e=e+1;
420      q=q/2;
421   }
422   number y=powerN(n,q,p);
423   number r=e;
424   number x=powerN(a,(q-1)/2,p);
425   number b=a*x^2 mod p;
426   x=a*x mod p;
427
428   while(((b-1) mod p)!=0)
429   {
430      m=0;z=b;
431      while(((z-1) mod p)!=0)
432      {
433         z=z^2 mod p;
434         m=m+1;
435      }
436      t=powerN(y,powerN(two,r-m-1,p),p);
437      y=t^2 mod p;
438      r=m;
439      x=x*t mod p;
440      b=b*y mod p;
441   }
442   return(x);
443}
444example
445{ "EXAMPLE:"; echo = 2;
446   ring R = 0,x,dp;
447   squareRoot(8315890421938608,32003);
448}
449
450
451proc solutionsMod2(matrix M)
452"USAGE:  solutionsMod2(M);
453RETURN: an intmat containing a basis of the vector space of solutions of the
454        linear system of equations defined by M over the prime field of
455        characteristic 2
456EXAMPLE:example solutionsMod2; shows an example
457"
458{
459   def R=basering;
460   ring Rhelp=2,z,(c,dp);
461   matrix M=imap(R,M);
462   matrix S=syz(M);
463   setring(R);
464   matrix S=imap(Rhelp,S);
465   int i,j;
466   intmat v[nrows(S)][ncols(S)];
467   for(i=1;i<=nrows(S);i++)
468   {
469      for(j=1;j<=ncols(S);j++)
470      {
471         if(S[i,j]==1){v[i,j]=1;}
472      }
473   }
474   return(v);
475}
476example
477{ "EXAMPLE:"; echo = 2;
478   ring R = 0,x,dp;
479   matrix M[3][3]=1,2,3,4,5,6,7,6,5;
480   solutionsMod2(M);
481}
482
483proc powerX(int q, int i, ideal I)
484"USAGE:  powerX(q,i,I);
485RETURN: the q-th power of the i-th variable modulo I
486ASSUME: I is a standard basis
487EXAMPLE:example powerX; shows an example
488"
489{
490   if(q<=181){return(reduce(var(i)^int(q),I));}
491   if((q mod 5)==0){return(reduce(powerX(q div 5,i,I)^5,I));}
492   if((q mod 4)==0){return(reduce(powerX(q div 4,i,I)^4,I));}
493   if((q mod 3)==0){return(reduce(powerX(q div 3,i,I)^3,I));}
494   if((q mod 2)==0){return(reduce(powerX(q div 2,i,I)^2,I));}
495   return(reduce(var(i)*powerX(q-1,i,I),I));
496}
497example
498{ "EXAMPLE:"; echo = 2;
499   ring R = 0,(x,y),dp;
500   powerX(100,2,std(ideal(x3-1,y2-x)));
501}
502
503//======================================================================
504//=========================== Discrete Logarithm =======================
505//======================================================================
506
507//============== Shank's baby step - giant step ========================
508
509proc babyGiant(number b, number y, number p)
510"USAGE:  babyGiant(b,y,p);
511RETURN: the discrete logarithm x: b^x=y mod p
512NOTE:   giant-step-baby-step
513EXAMPLE:example babyGiant; shows an example
514"
515{
516   int i,j,m;
517   list l;
518   number h=1;
519   number x;
520
521//choose m minimal such that m^2>p
522   for(i=1;i<=p;i++){if(i^2>p) break;}
523   m=i;
524
525//baby-step:  compute the table y*b^i for 1<=i<=m
526   for(i=1;i<=m;i++){l[i]=y*b^i mod p;}
527
528//giant-step: compute b^(m+j), 1<=j<=m and search in the baby-step table
529//for an i with y*b^i=b^(m*j). If found then x=m*j-i
530   number g=b^m mod p;
531   while(j<m)
532   {
533      j++;
534      h=h*g mod p;
535      for(i=1;i<=m;i++)
536      {
537         if(h==l[i])
538         {
539            x=m*j-i;
540            j=m;
541            break;
542         }
543     }
544   }
545   return(x);
546}
547example
548{ "EXAMPLE:"; echo = 2;
549   ring R = 0,z,dp;
550   number b=2;
551   number y=10;
552   number p=101;
553   babyGiant(b,y,p);
554}
555
556//==============  Pollards rho  =================================
557
558proc rho(number b, number y, number p)
559"USAGE:  rho(b,y,p);
560RETURN: the discrete logarithm x=log_b(y): b^x=y mod p
561NOTE: Pollard's rho:
562       choose random f_0 in 0,...,p-2 ,e_0=0, define x_0=b^f_0, define
563       x_i=y^e_ib^f_i as below. For i large enough there is i with
564       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,
565       d=gcd(s,p-1)=u*s+v*(p-1) then x=tu/d +j*(p-1)/d for some j (to be
566       found by trying)
567EXAMPLE:example rho; shows an example
568"
569{
570   int i=1;
571   int j;
572   number s,t;
573   list e,f,x;
574
575   e[1]=0;
576   f[1]=random(0,2147483629) mod (p-1);
577   x[1]=powerN(b,f[1],p);
578   while(i)
579   {
580      if((x[i] mod 3)==1)
581      {
582         x[i+1]=y*x[i] mod p;
583         e[i+1]=e[i]+1 mod (p-1);
584         f[i+1]=f[i];
585      }
586      if((x[i] mod 3)==2)
587      {
588         x[i+1]=x[i]^2 mod p;
589         e[i+1]=e[i]*2 mod (p-1);
590         f[i+1]=f[i]*2 mod (p-1);
591      }
592      if((x[i] mod 3)==0)
593      {
594         x[i+1]=x[i]*b mod p;
595         e[i+1]=e[i];
596         f[i+1]=f[i]+1 mod (p-1);
597      }
598      i++;
599      for(j=i-1;j>=1;j--)
600      {
601         if(x[i]==x[j])
602         {
603            s=(e[j]-e[i]) mod (p-1);
604            t=(f[i]-f[j]) mod (p-1);
605            if(s!=0)
606            {
607               i=0;
608            }
609            else
610            {
611               e[1]=0;
612               f[1]=random(0,2147483629) mod (p-1);
613               x[1]=powerN(b,f[1],p);
614               i=1;
615            }
616            break;
617         }
618      }
619   }
620
621   list w=exgcdN(s,p-1);
622   number u=w[1];
623   number d=w[3];
624
625   number a=(t*u/d) mod (p-1);
626
627   while(powerN(b,a,p)!=y)
628   {
629      a=(a+(p-1)/d) mod (p-1);
630   }
631   return(a);
632}
633example
634{ "EXAMPLE:"; echo = 2;
635   ring R = 0,x,dp;
636   number b=2;
637   number y=10;
638   number p=101;
639   rho(b,y,p);
640}
641//====================================================================
642//====================== Primality Tests =============================
643//====================================================================
644
645//================================= Miller-Rabin =====================
646
647proc MillerRabin(number n, int k)
648"USAGE:  MillerRabin(n,k);
649RETURN: 1 if n is prime, 0 else
650NOTE: probabilistic test of Miller-Rabin with k loops to test if n is prime.
651       Using the theorem:If n is prime, n-1=2^s*r, r odd, then
652       powerN(a,r,n)=1 or powerN(a,r*2^i,n)=-1 for some i
653EXAMPLE:example MillerRabin; shows an example
654"
655{
656   if(n<0){n=-n;}
657   if((n==2)||(n==3)){return(1);}
658   if((n mod 2)==0){return(0);}
659
660   int i;
661   number a,b,j,r,s;
662   r=n-1;
663   s=0;
664   while((r mod 2)==0)
665   {
666      s=s+1;
667      r=r/2;
668   }
669   while(i<k)
670   {
671      i++;
672      a=random(2,2147483629) mod n; if(a==0){a=3;}
673      if(exgcdN(a,n)[3]!=1){return(0);}
674      b=powerN(a,r,n);
675      if(b!=1)
676      {
677         j=0;
678         while(j<s)
679         {
680           if(((b+1) mod n)==0) break;
681           b=powerN(b,2,n);
682           j=j+1;
683         }
684         if(j==s){return(0);}
685      }
686   }
687   return(1);
688}
689example
690{ "EXAMPLE:"; echo = 2;
691   ring R = 0,z,dp;
692   number x=2;
693   x=x^787-1;
694   MillerRabin(x,3);
695}
696
697//======================= Soloway-Strassen  ==========================
698
699proc SolowayStrassen(number n, int k)
700"USAGE:  SolowayStrassen(n,k);
701RETURN: 1 if n is prime, 0 else
702NOTE: probabilistic test of Soloway-Strassen with k loops to test if n is
703       prime using the theorem:
704       If n is prime then powerN(a,(n-1)/2,n)=Jacobi(a,n) mod n
705EXAMPLE:example SolowayStrassen; shows an example
706"
707{
708   if(n<0){n=-n;}
709   if((n==2)||(n==3)){return(1);}
710   if((n mod 2)==0){return(0);}
711
712   number a;
713   int i;
714   while(i<k)
715   {
716      i++;
717      a=random(2,2147483629) mod n; if(a==0){a=3;}
718      if(gcdN(a,n)!=1){return(0);}
719      if(powerN(a,(n-1)/2,n)!=(Jacobi(a,n) mod n)){return(0);}
720   }
721   return(1);
722}
723example
724{ "EXAMPLE:"; echo = 2;
725   ring R = 0,z,dp;
726   number h=10;
727   number p=h^100+267;
728   //p=h^100+43723;
729   //p=h^200+632347;
730   SolowayStrassen(h,3);
731}
732
733
734/*
735ring R=0,z,dp;
736number p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
737number q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
738number n=p*q;
739SolowayStrassen(n,3);
740*/
741
742//===================== Pocklington-Lehmer ==============================
743
744proc PocklingtonLehmer(number N, list #)
745"USAGE: PocklingtonLehmer(N); optional: PocklingtonLehmer(N,L);
746        L a list of the first k primes
747RETURN:message N is not prime or {A,{p},{a_p}} as certificate for N being prime
748NOTE:assumes that it is possible to factorize N-1=A*B such that gcd(A,B)=1
749      the factorization of A is completely known and A^2>N .
750      N is prime if and only if for each prime factor p of A we can find
751      a_p such that a_p^(N-1)=1 mod N and gcd(a_p^((N-1)/p)-1,N)=1
752EXAMPLE:example PocklingtonLehmer; shows an example
753"
754{
755   number m=intRoot(N);
756   if(size(#)>0)
757   {
758      list S=PollardRho(N-1,10000,1,#);
759   }
760   else
761   {
762      list S=PollardRho(N-1,10000,1);
763   }
764   int i,j;
765   number A=1;
766   number p,a,g;
767   list PA;
768   list re;
769
770   while(i<size(S))
771   {
772      p=S[i+1];
773      A=A*p;
774      PA[i+1]=p;
775      if(A>m){break;}
776
777      while(1)
778      {
779        p=p*S[i+1];
780        if(((N-1) mod p)==0)
781        {
782           A=A*p;
783        }
784        else
785        {
786           break;
787        }
788      }
789      i++;
790   }
791   if(A<=m)
792   {
793     A=N/A;
794     PA=list(S[size(S)]);
795   }
796   for(i=1;i<=size(PA);i++)
797   {
798      a=1;
799      while(a<N-1)
800      {
801         a=a+1;
802         if(powerN(a,N-1,N)!=1){return("not prime");}
803         g=gcdN(powerN(a,(N-1)/PA[i],N),N);
804         if(g==1)
805         {
806           re[size(re)+1]=list(PA[i],a);
807           break;
808         }
809         if(g<N){"not prime";return(g);}
810      }
811   }
812   return(list(A,re));
813}
814example
815{ "EXAMPLE:"; echo = 2;
816   ring R = 0,z,dp;
817   number N=105554676553297;
818   PocklingtonLehmer(N);
819   list L=primList(1000);
820   PocklingtonLehmer(N,L);
821}
822
823//=======================================================================
824//======================= Factorization =================================
825//=======================================================================
826
827//======================= Pollards rho  =================================
828
829proc PollardRho(number n, int k, int allFactors, list #)
830"USAGE:  PollardRho(n,k,allFactors); optional: PollardRho(n,k,allFactors,L);
831         L a list of the first k primes
832RETURN: a list of factors of n (which could be just n),if allFactors=0
833        a list of all factors of n ,if allFactors=1
834NOTE: probabilistic rho-algorithm of Pollard to find a factor of n in k loops.
835      Creates a sequence x_i such that (x_i)^2=(x_2i)^2 mod n for some i,
836      computes gcd(x_i-x_2i,n) to find a divisor. To define the sequence
837      choose x,a and define x_n+1=x_n^2+a mod n, x_1=x.
838      If allFactors is 1, it tries to find recursively all prime factors
839      using the Soloway-Strassen test.
840EXAMPLE:example PollardRho; shows an example
841"
842{
843    int i,j;
844    list L=primList(100);
845    list re,se;
846    if(n<0){n=-n;}
847    if(n==1){return(re);}
848
849//this is optional: test whether a prime of the list # devides n
850    if(size(#)>0)
851    {
852       L=#;
853    }
854    for(i=1;i<=size(L);i++)
855    {
856       if((n mod L[i])==0)
857       {
858          re[size(re)+1]=L[i];
859          while((n mod L[i])==0)
860          {
861             n=n/L[i];
862          }
863       }
864       if(n==1){return(re);}
865    }
866    int e=size(re);
867//here the rho-algorithm starts
868    number a,d,x,y;
869    while(n>1)
870    {
871       a=random(2,2147483629);
872       x=random(2,2147483629);
873       y=x;
874       d=1;
875       i=0;
876       while(i<k)
877       {
878          i++;
879          x=powerN(x,2,n); x=(x+a) mod n;
880          y=powerN(y,2,n); y=(y+a) mod n;
881          y=powerN(y,2,n); y=(y+a) mod n;
882          d=gcdN(x-y,n);
883          if(d>1)
884          {
885             re[size(re)+1]=d;
886             while((n mod d)==0)
887             {
888               n=n/d;
889             }
890             break;
891          }
892          if(i==k)
893          {
894             re[size(re)+1]=n;
895             n=1;
896          }
897       }
898
899    }
900    if(allFactors)      //want to obtain all prime factors
901    {
902       i=e;
903       while(i<size(re))
904       {
905          i++;
906
907          if(!SolowayStrassen(re[i],5))
908          {
909             se=PollardRho(re[i],2*k,1);
910             re[i]=se[size(se)];
911             for(j=1;j<=size(se)-1;j++)
912             {
913                re[size(re)+1]=se[j];
914             }
915             i--;
916          }
917       }
918    }
919    return(re);
920}
921example
922{ "EXAMPLE:"; echo = 2;
923   ring R = 0,z,dp;
924   number h=10;
925   number p=h^30+4;
926   PollardRho(p,5000,0);
927}
928
929//======================== Pollards p-factorization ================
930proc pFactor(number n,int B, list P)
931"USAGE:  pFactor(n,B.P); n to be factorized, B a bound , P a list of primes
932RETURN: a list of factors of n or n if no factor found
933NOTE: Pollard's p-factorization
934       creates the product k of powers of primes (bounded by B)  from
935       the list P with the idea that for a prime divisor p of n p-1|k
936       then p devides gcd(a^k-1,n) for some random a
937EXAMPLE:example pFactor; shows an example
938"
939{
940   int i;
941   number k=1;
942   number w;
943   while(i<size(P))
944   {
945      i++;
946      w=P[i];
947      if(w>B) break;
948      while(w*P[i]<=B)
949      {
950         w=w*P[i];
951      }
952      k=k*w;
953   }
954   number a=random(2,2147483629);
955   number d=gcdN(powerN(a,k,n)-1,n);
956   if((d>1)&&(d<n)){return(d);}
957   return(n);
958}
959example
960{ "EXAMPLE:"; echo = 2;
961   ring R = 0,z,dp;
962   list L=primList(1000);
963   pFactor(1241143,13,L);
964   number h=10;
965   h=h^30+25;
966   pFactor(h,20,L);
967}
968
970
971proc quadraticSieve(number n, int c, list B, int k)
972"USAGE:  quadraticSieve(n,c,B,k); n to be factorized, [-c,c] the
973         sieve-intervall, B a list of primes,
974         k for using the first k elements in B
975RETURN: a list of factors of n or the message: no divisor found
976NOTE: quadraticSieve: Idea is to find x,y such that x^2=y^2 mod n then
977      gcd(x-y,n) can be a proper divisor of n
979"
980{
981   number f,d;
982   int i,j,l,s,p;
983   list S,tmp;
984   intvec v;
985   v[k]=0;
986
987//compute the integral part of the square root of n
988   number m=intRoot(n);
989
990//consider the function f(X)=(X+m)^2-n and compute for s in [-c,c] the values
991   while(i<=2*c)
992   {
993      f=(i-c+m)^2-n;
994      tmp[1]=i-c+m;
995      tmp[2]=f;
996      tmp[3]=v;
997      S[i+1]=tmp;
998      i++;
999   }
1000
1001//the sieve with p in B
1002//find all s in [-c,c] such that f(s) has all prime divisors in the first
1003//k elements of B and the decomposition of f(s). They are characterized
1004//by 1 or -1 at the second place of S[j]:
1005//S[j]=j-c+m,f(j-c)/p_1^v_1*...*p_k^v_k, v_1,...,v_k maximal
1006   for(i=1;i<=k;i++)
1007   {
1008      p=B[i];
1009      if((p>2)&&(Jacobi(n,p)==-1)){i++;continue;}//n is no quadratic rest mod p
1010      j=1;
1011      while(j<=p)
1012      {
1013         if(j>2*c+1) break;
1014         f=S[j][2];
1015         v=S[j][3];
1016         s=0;
1017         while((f mod p)==0)
1018         {
1019            s++;
1020            f=f/p;
1021         }
1022         if(s)
1023         {
1024           S[j][2]=f;
1025           v[i]=s;
1026           S[j][3]=v;
1027           l=j;
1028           while(l+p<=2*c+1)
1029           {
1030              l=l+p;
1031              f=S[l][2];
1032              v=S[l][3];
1033              s=0;
1034              while((f mod p)==0)
1035              {
1036                s++;
1037                f=f/p;
1038              }
1039              S[l][2]=f;
1040              v[i]=s;
1041              S[l][3]=v;
1042           }
1043         }
1044         j++;
1045      }
1046   }
1047   list T;
1048   for(j=1;j<=2*c+1;j++)
1049   {
1050      if((S[j][2]==1)||(S[j][2]==-1))
1051      {
1052         T[size(T)+1]=S[j];
1053      }
1054   }
1055
1056//the system of equations for the exponents {l_s} for the f(s) such
1057//product f(s)^l_s is a square (l_s are 1 or 0)
1058   matrix M[k+1][size(T)];
1059   for(j=1;j<=size(T);j++)
1060   {
1061      if(T[j][2]==-1){M[1,j]=1;}
1062      for(i=1;i<=k;i++)
1063      {
1064         M[i+1,j]=T[j][3][i];
1065      }
1066   }
1067   intmat G=solutionsMod2(M);
1068
1069//construction of x and y such that x^2=y^2 mod n and d=gcd(x-y,n)
1070//y=square root of product f(s)^l_s
1071//x=product s+m
1072   number x=1;
1073   number y=1;
1074
1075   for(i=1;i<=ncols(G);i++)
1076   {
1077      kill v;
1078      intvec v;
1079      v[k]=0;
1080      for(j=1;j<=size(T);j++)
1081      {
1082         x=x*T[j][1]^G[j,i] mod n;
1083         if((T[j][2]==-1)&&(G[j,i]==1)){y=-y;}
1084         v=v+G[j,i]*T[j][3];
1085
1086      }
1087      for(l=1;l<=k;l++)
1088      {
1089         y=y*B[l]^(v[l]/2) mod n;
1090      }
1091      d=gcdN(x-y,n);
1092      if((d>1)&&(d<n)){return(d);}
1093   }
1094   return("no divisor found");
1095}
1096example
1097{ "EXAMPLE:"; echo = 2;
1098   ring R = 0,z,dp;
1099   list L=primList(5000);
1102}
1103
1104//======================================================================
1105//==================== elliptic curves  ================================
1106//======================================================================
1107
1108//================= elementary operations ==============================
1109
1110proc isOnCurve(number N, number a, number b, list P)
1111"USAGE:  isOnCurve(N,a,b,P);
1112RETURN: 1 or 0 (depending on whether P is on the curve or not)
1113NOTE: checks whether P=(P[1]:P[2]:P[3]) is a point on the elliptic
1114      curve defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1115EXAMPLE:example isOnCurve; shows an example
1116"
1117{
1118   if(((P[2]^2*P[3]-P[1]^3-a*P[1]*P[3]^2-b*P[3]^3) mod N)!=0){return(0);}
1119   return(1);
1120}
1121example
1122{ "EXAMPLE:"; echo = 2;
1123   ring R = 0,z,dp;
1124   isOnCurve(32003,5,7,list(10,16,1));
1125}
1126
1127proc ellipticAdd(number N, number a, number b, list P, list Q)
1129RETURN: list L, representing the point P+Q
1130NOTE: P=(P[1]:P[2]:P[3]),Q =(Q[1]:Q[2]:Q[3])points on the elliptic curve
1131      defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1133"
1134{
1135   if(N==2){ERROR("not implemented for 2");}
1136   int i;
1137   for(i=1;i<=3;i++)
1138   {
1139      P[i]=P[i] mod N;
1140      Q[i]=Q[i] mod N;
1141   }
1142   list Resu;
1143   Resu[1]=number(0);
1144   Resu[2]=number(1);
1145   Resu[3]=number(0);
1146   list Error;
1147   Error[1]=0;
1148   //test for ellictic curve
1149   number D=4*a^3+27*b^2;
1150   number g=gcdN(D,N);
1151   if(g==N){return(Error);}
1152   if(g!=1)
1153   {
1154      P[4]=g;
1155      return(P);
1156   }
1157   if(((P[1]==0)&&(P[2]==0)&&(P[3]==0))||((Q[1]==0)&&(Q[2]==0)&&(Q[3]==0)))
1158   {
1159      Error[1]=-2;
1160      return(Error);
1161   }
1162   if(!isOnCurve(N,a,b,P)||!isOnCurve(N,a,b,Q))
1163   {
1164      Error[1]=-1;
1165      return(Error);
1166   }
1167   if(P[3]==0){return(Q);}
1168   if(Q[3]==0){return(P);}
1169   list I=exgcdN(P[3],N);
1170   if(I[3]!=1)
1171   {
1172      P[4]=I[3];
1173      return(P);
1174   }
1175   P[1]=P[1]*I[1] mod N;
1176   P[2]=P[2]*I[1] mod N;
1177   I=exgcdN(Q[3],N);
1178   if(I[3]!=1)
1179   {
1180      P[4]=I[3];
1181      return(P);
1182   }
1183   Q[1]=Q[1]*I[1] mod N;
1184   Q[2]=Q[2]*I[1] mod N;
1185   if((P[1]==Q[1])&&(((P[2]+Q[2]) mod N)==0)){return(Resu);}
1186   number L;
1187   if((P[1]==Q[1])&&(P[2]==Q[2]))
1188   {
1189      I=exgcdN(2*Q[2],N);
1190      if(I[3]!=1)
1191      {
1192         P[4]=I[3];
1193         return(P);
1194      }
1195      L=I[1]*(3*Q[1]^2+a) mod N;
1196   }
1197   else
1198   {
1199      I=exgcdN(Q[1]-P[1],N);
1200      if(I[3]!=1)
1201      {
1202         P[4]=I[3];
1203         return(P);
1204      }
1205      L=(Q[2]-P[2])*I[1] mod N;
1206   }
1207   Resu[1]=(L^2-P[1]-Q[1]) mod N;
1208   Resu[2]=(L*(P[1]-Resu[1])-P[2]) mod N;
1209   Resu[3]=number(1);
1210   return(Resu);
1211}
1212example
1213{ "EXAMPLE:"; echo = 2;
1214   ring R = 0,z,dp;
1215   number N=11;
1216   number a=1;
1217   number b=6;
1218   list P,Q;
1219   P[1]=2;
1220   P[2]=4;
1221   P[3]=1;
1222   Q[1]=3;
1223   Q[2]=5;
1224   Q[3]=1;
1226}
1227
1228proc ellipticMult(number N, number a, number b, list P, number k)
1229"USAGE:  ellipticMult(N,a,b,P,k);
1230RETURN: a list L representing the point k*P
1231NOTE:  P=(P[1]:P[2]:P[3]) a point on the elliptic curve defined by
1232       y^2z=x^3+a*xz^2+b*z^3  over Z/N
1233EXAMPLE:example ellipticMult; shows an example
1234"
1235{
1236   if(P[3]==0){return(P);}
1237   list resu;
1238   resu[1]=number(0);
1239   resu[2]=number(1);
1240   resu[3]=number(0);
1241
1242   if(k==0){return(resu);}
1243   if(k==1){return(P);}
1245   if(k==-1)
1246   {
1247      resu=P;
1248      resu[2]=N-P[2];
1249      return(resu);
1250   }
1251   if(k<0)
1252   {
1253      resu=ellipticMult(N,a,b,P,-k);
1254      return(ellipticMult(N,a,b,resu,-1));
1255   }
1256   if((k mod 2)==0)
1257   {
1258      resu=ellipticMult(N,a,b,P,k/2);
1260   }
1261   resu=ellipticMult(N,a,b,P,k-1);
1263}
1264example
1265{ "EXAMPLE:"; echo = 2;
1266   ring R = 0,z,dp;
1267   number N=11;
1268   number a=1;
1269   number b=6;
1270   list P;
1271   P[1]=2;
1272   P[2]=4;
1273   P[3]=1;
1274   ellipticMult(N,a,b,P,3);
1275}
1276
1277//================== Random for elliptic curves =====================
1278
1279proc ellipticRandomCurve(number N)
1280"USAGE:  ellipticRandomCurve(N);
1281RETURN: a list of two random numbers a,b and 4a^3+27b^2 mod N
1282NOTE:   y^2z=x^3+a*xz^2+b^2*z^3 defines an elliptic curve over Z/N
1283EXAMPLE:example ellipticRandomCurve; shows an example
1284"
1285{
1286   int k;
1287   while(k<=10)
1288   {
1289     k++;
1290     number a=random(1,2147483647) mod N;
1291     number b=random(1,2147483647) mod N;
1292     //test for ellictic curve
1293     number D=4*a^3+27*b^4; //the constant term is b^2
1294     number g=gcdN(D,N);
1295     if(g<N){return(list(a,b,g));}
1296   }
1297   ERROR("no random curve found");
1298}
1299example
1300{ "EXAMPLE:"; echo = 2;
1301   ring R = 0,z,dp;
1302   ellipticRandomCurve(32003);
1303}
1304
1305proc ellipticRandomPoint(number N, number a, number b)
1306"USAGE:  ellipticRandomPoint(N,a,b);
1307RETURN: a list representing  a random point (x:y:z) of the elliptic curve
1308        defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1309EXAMPLE:example ellipticRandomPoint; shows an example
1310"
1311{
1312   number x=random(1,2147483647) mod N;
1313   number h=x^3+a*x+b;
1314   h=h mod N;
1315   list resu;
1316   resu[1]=x;
1317   resu[2]=0;
1318   resu[3]=1;
1319   if(h==0){return(resu);}
1320
1321   number n=Jacobi(h,N);
1322   if(n==0)
1323   {
1324      resu=-5;
1325      "N is not prime";
1326      return(resu);
1327   }
1328   if(n==1)
1329   {
1330      resu[2]=squareRoot(h,N);
1331      return(resu);
1332   }
1333   return(ellipticRandomPoint(N,a,b));
1334}
1335example
1336{ "EXAMPLE:"; echo = 2;
1337   ring R = 0,z,dp;
1338   ellipticRandomPoint(32003,3,181);
1339}
1340
1341
1342
1343//====================================================================
1344//======== counting the points of an elliptic curve  =================
1345//====================================================================
1346
1347//==================   the trivial approaches  =======================
1348proc countPoints(number N, number a, number b)
1349"USAGE:  countPoints(N,a,b);
1350RETURN: the number of points of the elliptic curve defined by
1351        y^2=x^3+a*x+b  over Z/N
1352NOTE: trivial aproach
1353EXAMPLE:example countPoints; shows an example
1354"
1355{
1356  number x;
1357  number r=N+1;
1358  while(x<N)
1359  {
1360     r=r+Jacobi((x^3+a*x+b) mod N,N);
1361     x=x+1;
1362  }
1363  return(r);
1364}
1365example
1366{ "EXAMPLE:"; echo = 2;
1367   ring R = 0,z,dp;
1368   countPoints(181,71,150);
1369}
1370
1371proc ellipticAllPoints(number N, number a, number b)
1372"USAGE:  ellipticAllPoints(N,a,b);
1373RETURN: list of points (x:y:z) of the elliptic curve defined by
1374        y^2z=x^3+a*xz^2+b*z^3  over Z/N
1375EXAMPLE:example ellipticAllPoints; shows an example
1376"
1377{
1378   list resu,point;
1379   point[1]=0;
1380   point[2]=1;
1381   point[3]=0;
1382   resu[1]=point;
1383   point[3]=1;
1384   number x,h,n;
1385   while(x<N)
1386   {
1387      h=(x^3+a*x+b) mod N;
1388      if(h==0)
1389      {
1390         point[1]=x;
1391         point[2]=0;
1392         resu[size(resu)+1]=point;
1393      }
1394      else
1395      {
1396         n=Jacobi(h,N);
1397         if(n==1)
1398         {
1399            n=squareRoot(h,N);
1400            point[1]=x;
1401            point[2]=n;
1402            resu[size(resu)+1]=point;
1403            point[2]=N-n;
1404            resu[size(resu)+1]=point;
1405         }
1406      }
1407      x=x+1;
1408   }
1409   return(resu);
1410}
1411example
1412{ "EXAMPLE:"; echo = 2;
1413   ring R = 0,z,dp;
1414   list L=ellipticAllPoints(181,71,150);
1415   size(L);
1416   L[size(L)];
1417}
1418
1419//================ the algorithm of Shanks and Mestre =================
1420
1421proc ShanksMestre(number q, number a, number b, list #)
1422"USAGE:  ShanksMestre(q,a,b); optional:ShanksMestre(q,a,b,s); s the number
1423         of loops in the algorithm (default s=1)
1424RETURN: the number of points of the elliptic curve defined by
1425         y^2=x^3+a*x+b  over Z/N
1426NOTE: algorithm of Shanks and Mestre (giant-step-baby-step)
1427EXAMPLE:example ShanksMestre; shows an example
1428"
1429{
1430   number n=intRoot(4*q);
1431   number m=intRoot(intRoot(16*q))+1;
1432   number d;
1433   int i,j,k,s;
1434   list B,K,T,P,Q,R,mP;
1435   B[1]=list(0,1,0);
1436   if(size(#)>0)
1437   {
1438      s=#[1];
1439   }
1440   else
1441   {
1442      s=1;
1443   }
1444   while(k<s)
1445   {
1446      P =ellipticRandomPoint(q,a,b);
1447      Q =ellipticMult(q,a,b,P,n+q+1);
1448
1449      while(j<m)
1450      {
1451         j++;
1453      }
1455      mP[2]=q-mP[2];
1456      while(i<m)                            //giant-step
1457      {
1458         j=0;
1459         while(j<m)
1460         {
1461            j=j+1;
1462            if((Q[1]==B[j][1])&&(Q[2]==B[j][2])&&(Q[3]==B[j][3]))
1463            {
1464
1465               T[1]=P;
1466               T[2]=q+1+n-(i*m+j-1);
1467               K[size(K)+1]=T;
1468               if(size(K)>1)
1469               {
1470                  if(K[size(K)][2]!=K[size(K)-1][2])
1471                  {
1472                     d=gcdN(K[size(K)][2],K[size(K)-1][2]);
1473                     if(ellipticMult(q,a,b,K[size(K)],d)[3]==0)
1474                     {
1475                       K[size(K)][2]=K[size(K)-1][2];
1476                     }
1477                  }
1478               }
1479               i=int(m);
1480               break;
1481            }
1482         }
1483         i=i+1;
1485      }
1486      k++;
1487   }
1488   if(size(K)>0)
1489   {
1490      int te=1;
1491      for(i=1;i<=size(K)-1;i++)
1492      {
1493         if(K[size(K)][2]!=K[i][2])
1494         {
1495           if(ellipticMult(q,a,b,K[i],K[size(K)][2])[3]!=0)
1496           {
1497             te=0;
1498             break;
1499           }
1500         }
1501      }
1502      if(te)
1503      {
1504         return(K[size(K)][2]);
1505      }
1506   }
1507   return(ShanksMestre(q,a,b,s));
1508}
1509example
1510{ "EXAMPLE:"; echo = 2;
1511   ring R = 0,z,dp;
1512   ShanksMestre(32003,71,602);
1513}
1514
1515//==================== Schoof's algorithm =============================
1516
1517proc Schoof(number N,number a, number b)
1518"USAGE:  Schoof(N,a,b);
1519RETURN: the number of points of the elliptic curve defined by
1520        y^2=x^3+a*x+b  over Z/N
1521NOTE:  algorithm of Schoof
1522EXAMPLE:example Schoof; shows an example
1523"
1524{
1525   int pr=printlevel;
1526//test for ellictic curve
1527   number D=4*a^3+27*b^2;
1528   number G=gcdN(D,N);
1529   if(G==N){ERROR("not an elliptic curve");}
1530   if(G!=1){ERROR("not a prime");}
1531
1532//=== small N
1533  // if((N<=500)&&(pr<5)){return(countPoints(int(N),a,b));}
1534
1535//=== the general case
1536   number q=intRoot(4*N);
1537   list L=primL(2*q);
1538   int r=size(L);
1539   list T;
1540   int i,j;
1541   for(j=1;j<=r;j++)
1542   {
1543      T[j]=(testElliptic(int(N),a,b,L[j])+int(q)) mod L[j];
1544   }
1545   if(pr>=5)
1546   {
1547      "===================================================================";
1548      "Chinese remainder :";
1549      for(i=1;i<=size(T);i++)
1550      {
1551         " x =",T[i]," mod ",L[i];
1552      }
1553      "gives t+ integral part of the square root of q (to be positive)";
1554      chineseRem(T,L);
1555      "we obtain t = ",chineseRem(T,L)-q;
1556      "===================================================================";
1557   }
1558   number t=chineseRem(T,L)-q;
1559   return(N+1-t);
1560}
1561example
1562{ "EXAMPLE:"; echo = 2;
1563   ring R = 0,z,dp;
1564   Schoof(32003,71,602);
1565}
1566
1567/*
1568needs 518 sec
1569Schoof(2147483629,17,3567);
15702147168895
1571*/
1572
1573
1574proc generateG(number a,number b, int m)
1575"USAGE:  generateG(a,b,m);
1576RETURN: m-th division polynomial
1577NOTE: generate the recursively defined polynomials in Z[x,y],so called
1578 division polynomials, p_m=generateG(a,b,m) such that on the elliptic curve
1579 defined by y^2=x^3+a*x+b  over Z/N and a point P=(x:y:1) the point m*P is
1580 (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)
1581 m*P=0  iff p_m(P)=0
1582EXAMPLE:example generateG; shows an example
1583"
1584{
1585   poly f;
1586   if(m==0){return(f);}
1587   if(m==1){return(1);}
1588   if(m==2){f=2*var(1);return(f);}
1589   if(m==3){f=3*var(2)^4+6*a*var(2)^2+12*b*var(2)-a^2;return(f);}
1590   if(m==4)
1591   {
1592      f=4*var(1)*(var(2)^6+5*a*var(2)^4+20*b*var(2)^3-5*a^2*var(2)^2
1593        -4*a*b*var(2)-8*b^2-a^3);
1594      return(f);
1595   }
1596   if((m mod 2)==0)
1597   {
1598      f=(generateG(a,b,m/2+2)*generateG(a,b,m/2-1)^2
1599        -generateG(a,b,m/2-2)*generateG(a,b,m/2+1)^2)
1600        *generateG(a,b,m/2)/(2*var(1));
1601      return(f);
1602   }
1603   f=generateG(a,b,(m-1)/2+2)*generateG(a,b,(m-1)/2)^3
1604     -generateG(a,b,(m-1)/2-1)*generateG(a,b,(m-1)/2+1)^3;
1605   return(f);
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),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)/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 decode(string g)
2079"USAGE:  decode(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: decode;  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=decode(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.