source: git/Singular/LIB/crypto.lib @ bee06d

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