source: git/Singular/LIB/krypto.lib @ 8265bdc

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