source:git/Singular/LIB/crypto.lib@0bd2cb

spielwiese
Last change on this file since 0bd2cb was 0bd2cb, checked in by Hans Schoenemann <hannes@…>, 9 years ago
• Property mode set to `100644`
File size: 90.4 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="version crypto.lib 4.0.0.1 Aug_2014 "; // \$Id\$
3category="Teaching";
4info="
5LIBRARY:  crypto.lib     Procedures for teaching cryptography
6AUTHORS:                 Gerhard Pfister, pfister@mathematik.uni-kl.de
7@*                       David Brittinger, dativ@gmx.net
8
9OVERVIEW:
10     The library contains procedures to compute the discrete logarithm,
11      primality-tests, factorization included elliptic curves.
12      The library is intended to be used for teaching purposes but not
13      for serious computations. Sufficiently high printlevel allows to
14      control each step, thus illustrating the algorithms at work.
15
16
17PROCEDURES:
18 decimal(s);                number corresponding to the hexadecimal number s
19 eexgcdN(L)                 T with sum_i L[i]*T[i]=T[n+1]=gcd(L[1],...,L[n])
20 lcmN(a,b)                  compute lcm(a,b)
21 powerN(m,d,n)              compute m^d mod n
22 chineseRem(T,L)            compute x such that x = T[i] mod L[i]
23 Jacobi(a,n)                the generalized Legendre symbol of a and n
24 primList(n)                the list of all primes <=n
25 primL(q)                   first primes p_1,...,p_r such that q<p_1*...*p_r
26 intPart(x)                 the integral part of a rational number
27 intRoot(m)                 the integral part of the square root of m
28 squareRoot(a,p)            the square root of a in Z/p, p prime
29 solutionsMod2(M)           basis solutions of Mx=0 over Z/2
30 powerX(q,i,I)              q-th power of the i-th variable modulo I
31 babyGiant(b,y,p)           discrete logarithm x: b^x=y mod p
32 rho(b,y,p)                 discrete logarithm x: b^x=y mod p
33 MillerRabin(n,k)           probabilistic primaly-test of Miller-Rabin
34 SolowayStrassen(n,k)       probabilistic primaly-test of Soloway-Strassen
35 PocklingtonLehmer(N,[])    primaly-test of Pocklington-Lehmer
36 PollardRho(n,k,a,[])       Pollard's rho factorization
37 pFactor(n,B,P)             Pollard's p-factorization
39 isOnCurve(N,a,b,P)         P is on the curve y^2z=x^3+a*xz^2+b*z^3  over Z/N
41 ellipticMult(N,a,b,P,k)    k*P on elliptic curves
42 ellipticRandomCurve(N)     generates y^2z=x^3+a*xz^2+b*z^3  over Z/N randomly
43 ellipticRandomPoint(N,a,b) random point on y^2z=x^3+a*xz^2+b*z^3  over Z/N
44 countPoints(N,a,b)         number of points of y^2=x^3+a*x+b  over Z/N
45 ellipticAllPoints(N,a,b)   points of y^2=x^3+a*x+b  over Z/N
46 ShanksMestre(q,a,b,[])     number of points of y^2=x^3+a*x+b  over Z/N
47 Schoof(N,a,b)              number of points of y^2=x^3+a*x+b  over Z/N
48 generateG(a,b,m)           m-th division polynomial of y^2=x^3+a*x+b  over Z/N
49 factorLenstraECM(N,S,B,[]) Lenstra's factorization
50 ECPP(N)                    primaly-test of Goldwasser-Kilian
51 calculate_ordering(num1, primitive, mod1)    Calculates x so that primitive^x == num1 mod mod1
52  is_primitive_root(primitive, mod1) Checks if primitive is a primitive root modulo mod1
53  find_first_primitive_root(mod1) Returns the first primitive root modulo mod1, starting with 1
55  inverse_modulus(num,mod1)       Finds a t so that t*num = 1 mod mod1
56  is_prime(n)                     Checks if n is prime
57  proc find_biggest_index(a)      Returns the index of the biggest element of a
58  find_index(a,e)                 Returns the list index of element e in list a. Returns 0 if e is not in a
59  subset_sum01(list knapsack, int solution)                          solves the subset-sum-knapsack-problem by calculating all subsets and choosing the right solution
60  subset_sum02(list knapsack, int sol)                            solves the subset-sum-knapsack-problem with a naive greedy algorithm
61  unbounded_knapsack(list knapsack, list profit, int capacity)                solves the unbounded_knapsack-problem, needing a list of knapsack weights, a list of profits and a capacity
62  multidimensional_knapsack(matrix m, list capacities, list profits)              solves the multidimensional_knapsack-problem by using the PECH algorithm, needing a weight matrix m, a list of capacities and a list of profits
63  naccache_stern_generation(int key, int primenum)                      generates a hard knapsack for the Naccache-Stern Kryptosystem for given key and prime modulus
64  naccache_stern_encryption(list knapsack, list message, int primenum)            encrypts a message with the Naccache-Stern Kryptosystem, using a hard knapsack, a message encoded as binary list and a prime modulus
65  naccache_stern_decryption(list knapsack, int key, int primenum, int message)        decrypts a message with the Naccache-Stern Kryptosystem, using the easy knapsack, the key, the prime modulus and the message encoded as integer
66  m_merkle_hellman_transformation(list knapsack, int primitive, int mod1)            generates a hard knapsack for the multiplicative Merkle-Hellman Kryptosystem for a given easy knapsack and a primitive root for a modulus mod1
67  m_merkle_hellman_encryption(list knapsack, list message)                  encrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using a hard knapsack and a message encoded as binary list
68  m_merkle_hellman_decryption(list knapsack, bigint primitive, bigint mod1, int message)    decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the easy knapsack, the key given by the primitive root, the modulus mod1 and the message encoded as integer
69  merkle_hellman_transformation(list knapsack, int key, int mod1                generates a hard knapsack for the  Merkle-Hellman Kryptosystem for a given easy knapsack , a multiplicator key and a modulus mod1
70  merkle_hellman_encryption(list knapsack, list message)                    encrypts a message with the Merkle-Hellman Kryptosystem, using a hard knapsack and a message encoded as binary list
71  merkle_hellman_decryption(list knapsack, int key, int mod1, int message)          decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the hard knapsack, the key, the modulus mod1 and the message encoded as integer
72  super_increasing_knapsack(int ksize)                            Creates the smallest super-increasing knapsack of given size ksize
73  h_increasing_knapsack(int ksize, int h)                            Creates the smallest h-increasing knapsack of given size ksize and h
74  injective_knapsack(int ksize, int kmaxelement)                        Creates all list of all injective knapsacks of given size ksize and maximal element kmaxelement
75  calculate_max_sum(list a)                                  Calculates the maximal sum of a given knapsack a
76  is_injective(list a)                                    Checks if knapsack a is injective
77  is_h_injective(list a, int h)                                Checks if knapsack a is h-injective
78  is_fix_injective(list a)                                  Checks if knapsack a is fix-injective
79  three_elements(list out, int iterations)                          Creates the smallest injective knapsack with a given injective_knapsack by using the three-elements-algorithm with a given number of iterations
80
81              [parameters in square brackets are optional]
82";
83
84LIB "poly.lib";
85LIB "atkins.lib";
86
87///////////////////////////////////////////////////////////////////////////////
88
89
90//=============================================================================
91//=========================== basic prozedures ================================
92//=============================================================================
93
94proc decimal(string s)
95"USAGE:  decimal(s); s = string
96RETURN: the (decimal) number corresponding to the hexadecimal number s
97EXAMPLE:example decimal; shows an example
98"
99{
100   int n=size(s);
101   int i;
102   bigint k;
103   bigint t=16;
104   bigint m=0;
105   for(i=1;i<=n;i++)
106   {
107      k=0;
108      if(s[i]=="1"){k=1;}
109      if(s[i]=="2"){k=2;}
110      if(s[i]=="3"){k=3;}
111      if(s[i]=="4"){k=4;}
112      if(s[i]=="5"){k=5;}
113      if(s[i]=="6"){k=6;}
114      if(s[i]=="7"){k=7;}
115      if(s[i]=="8"){k=8;}
116      if(s[i]=="9"){k=9;}
117      if(s[i]=="a"){k=10;}
118      if(s[i]=="b"){k=11;}
119      if(s[i]=="c"){k=12;}
120      if(s[i]=="d"){k=13;}
121      if(s[i]=="e"){k=14;}
122      if(s[i]=="f"){k=15;}
123      m=m*t+k;
124   }
125   return(m);
126}
127example
128{ "EXAMPLE:"; echo = 2;
129   string s  ="8edfe37dae96cfd2466d77d3884d4196";
130   decimal(s);
131}
132
133proc eexgcdN(list L)
134"USAGE:  eexgcdN(L);
135RETURN: list T such that sum_i L[i]*T[i]=T[n+1]=gcd(L[1],...,L[n])
136EXAMPLE:example eexgcdN; shows an example
137"
138{
139   if(size(L)==2)
140   {
141     list LL=extgcd(L[1],L[2]);return(list(LL[2],LL[3],LL[1]));
142   }
143   bigint p=L[size(L)];
144   L=delete(L,size(L));
145   list T=eexgcdN(L);
146   list S=extgcd(T[size(T)],p);
147   int i;
148   for(i=1;i<=size(T)-1;i++)
149   {
150      T[i]=T[i]*S[2];
151   }
152   p=T[size(T)];
153   T[size(T)]=S[3];
154   T[size(T)+1]=S[1];
155   return(T);
156}
157example
158{ "EXAMPLE:"; echo = 2;
159   eexgcdN(list(24,15,21));
160}
161
162proc lcmN(bigint a, bigint b)
163"USAGE:  lcmN(a,b);
164RETURN: lcm(a,b);
165EXAMPLE:example lcmN; shows an example
166"
167{
168   return (a*b/gcd(a,b));
169}
170example
171{ "EXAMPLE:"; echo = 2;
172   lcmN(24,15);
173}
174
175proc powerN(bigint m, bigint d, bigint n)
176"USAGE:  powerN(m,d,n);
177RETURN: m^d mod n
178EXAMPLE:example powerN; shows an example
179"
180{
181   if(d==0){return(bigint(1));}
182   int i;
183   if(n==0)
184   {
185      for(i=12;i>=2;i--)
186      {
187         if((d mod i)==0){return(powerN(m,d div i,n)^i);}
188      }
189      return(m*powerN(m,d-1,n));
190   }
191   for(i=12;i>=2;i--)
192   {
193      if((d mod i)==0)
194      {
195        bigint rr=powerN(m,d div i,n)^i mod n;
196        if (rr<0) { rr=rr+n;}
197        return(rr);
198      }
199   }
200   return(m*powerN(m,d-1,n) mod n);
201}
202example
203{ "EXAMPLE:"; echo = 2;
204   powerN(24,15,7);
205}
206
207proc chineseRem(list T,list L)
208"USAGE:  chineseRem(T,L);
209RETURN: x such that x = T[i] mod L[i]
210NOTE:   chinese remainder theorem
211EXAMPLE:example chineseRem; shows an example
212"
213{
214   int i;
215   int n=size(L);
216   bigint N=1;
217   for(i=1;i<=n;i++)
218   {
219      N=N*L[i];
220   }
221   list M;
222   for(i=1;i<=n;i++)
223   {
224      M[i]=N div L[i];
225   }
226   list S=eexgcdN(M);
227   bigint x;
228   for(i=1;i<=n;i++)
229   {
230      x=x+S[i]*M[i]*T[i];
231   }
232   x=x mod N;
233   return(x);
234}
235example
236{ "EXAMPLE:"; echo = 2;
237   chineseRem(list(24,15,7),list(2,3,5));
238}
239
240proc Jacobi(bigint a, bigint n)
241"USAGE:  Jacobi(a,n);
242RETURN: the generalized Legendre symbol
243NOTE: if n is an odd prime then Jacobi(a,n)=0,1,-1 if n|a, a=x^2 mod n,else
244EXAMPLE:example Jacobi; shows an example
245"
246{
247   int i;
248   int z=1;
249   bigint t=1;
250   bigint k;
251
252   if((((n-1) div 2) mod 2)!=0){z=-1;}
253   if(a<0){return(z*Jacobi(-a,n));}
254   a=a mod n;
255   if(n==1){return(1);}
256   if(a==0){return(0);}
257
258   while(a!=0)
259   {
260      while((a mod 2)==0)
261      {
262         a=a div 2;
263         if(((n mod 8)==3)||((n mod 8)==5)){t=-t;}
264      }
265      k=a;a=n;n=k;
266      if(((a mod 4)==3)&&((n mod 4)==3)){t=-t;}
267      a=a mod n;
268   }
269   if (n==1){return(t);}
270   return(0);
271}
272example
273{ "EXAMPLE:"; echo = 2;
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(bigint 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   bigint 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   primL(20);
339}
340
341proc intPart(number x)
342"USAGE:  intPart(x);
343RETURN: the integral part of a rational number
344EXAMPLE:example intPart; shows an example
345"
346{
347  if (x>=0)
348  {
349    return(bigint((numerator(x)-(bigint(numerator(x)) mod bigint(denominator(x)))))
350         div bigint(denominator(x)));
351  }
352  else
353  {
354    return(bigint((numerator(x)-(bigint(numerator(x)) mod bigint(denominator(x)+denominator(x)))))
355         div bigint(denominator(x)));
356  }
357}
358example
359{ "EXAMPLE:"; echo = 2;
360   ring r=0,x,dp;
361   intPart(7/3);
362}
363
364proc intRoot(bigint mm)
365"USAGE:  intRoot(m);
366RETURN: the integral part of the square root of m
367EXAMPLE:example intRoot; shows an example
368"
369{
370   ring R = 0,@x,dp;
371   number m=mm;
372   number x=1;
373   number t=x^2;
374   number s=(x+1)^2;
375   while(((m>t)&&(m>s))||((m<t)&&(m<s)))
376   {
377      if (x==0) { t=0; }
378      else
379      {
380        x=intPart(x/2+m/(2*x));  //Newton step
381        t=x^2;
382      }
383      if(t>m)
384      {
385         s=(x-1)^2;
386      }
387      else
388      {
389         s=(x+1)^2;
390      }
391   }
392   if(t>m){return(bigint(x-1));}
393   if(s==m){return(bigint(x+1));}
394   return(bigint(x));
395}
396example
397{ "EXAMPLE:"; echo = 2;
398   intRoot(20);
399}
400
401proc squareRoot(bigint a, bigint p)
402"USAGE:  squareRoot(a,p);
403RETURN: the square root of a in Z/p, p prime
404NOTE:   assumes the Jacobi symbol is 1 or p=2.
405EXAMPLE:example squareRoot; shows an example
406"
407{
408   if(p==2){return(a);}
409   if((a mod p)==0){return(0);}
410   if (a<0) { a=a+p; }
411   if(powerN(a,p-1,p)!=1)
412   {
413      "p is not prime";
414      return(bigint(-5));
415   }
416   bigint n=random(1,2147483647) mod p;
417   if(n==0){n=n+1;}
418   bigint j=Jacobi(n,p);
419   if(j==0)
420   {
421      "p is not prime";
422      return(bigint(-5));
423   }
424   if(j==1)
425   {
426      return(squareRoot(a,p));
427   }
428   bigint q=p-1;
429   bigint e=0;
430   bigint two=2;
431   bigint z,m,t;
432   while((q mod 2)==0)
433   {
434      e=e+1;
435      q=q div 2;
436   }
437   bigint y=powerN(n,q,p);
438   bigint r=e;
439   bigint x=powerN(a,(q-1) div 2,p);
440   bigint b=a*x^2 mod p;
441   x=a*x mod p;
442
443   while(((b-1) mod p)!=0)
444   {
445      m=0;z=b;
446      while(((z-1) mod p)!=0)
447      {
448         z=z^2 mod p;
449         m=m+1;
450      }
451      t=powerN(y,powerN(two,r-m-1,p),p);
452      y=t^2 mod p;
453      r=m;
454      x=x*t mod p;
455      b=b*y mod p;
456   }
457   return(x);
458}
459example
460{ "EXAMPLE:"; echo = 2;
461   squareRoot(8315890421938608,32003);
462}
463
464
465proc solutionsMod2(bigintmat MM)
466"USAGE:  solutionsMod2(M);
467RETURN: an intmat containing a basis of the vector space of solutions of the
468        linear system of equations defined by M over the prime field of
469        characteristic 2
470EXAMPLE:example solutionsMod2; shows an example
471"
472{
473   ring Rhelp=2,z,(c,dp);
474   int i,j;
475   matrix M[nrows(MM)][ncols(MM)];
476   for(i=1;i<=nrows(MM);i++)
477   {
478      for(j=1;j<=ncols(MM);j++)
479      {
480         M[i,j]=MM[i,j];
481      }
482   }
483   matrix S=syz(M);
484   intmat v[nrows(S)][ncols(S)];
485   for(i=1;i<=nrows(S);i++)
486   {
487      for(j=1;j<=ncols(S);j++)
488      {
489         if(S[i,j]==1){v[i,j]=1;}
490      }
491   }
492   return(v);
493}
494example
495{ "EXAMPLE:"; echo = 2;
496   bigintmat M[3][3]=1,2,3,4,5,6,7,6,5;
497   solutionsMod2(M);
498}
499
500proc powerX(int q, int i, ideal I)
501"USAGE:  powerX(q,i,I);
502RETURN: the q-th power of the i-th variable modulo I
503ASSUME: I is a standard basis
504EXAMPLE:example powerX; shows an example
505"
506{
507   if(q<=181){return(reduce(var(i)^int(q),I));}
508   if((q mod 5)==0){return(reduce(powerX(q div 5,i,I)^5,I));}
509   if((q mod 4)==0){return(reduce(powerX(q div 4,i,I)^4,I));}
510   if((q mod 3)==0){return(reduce(powerX(q div 3,i,I)^3,I));}
511   if((q mod 2)==0){return(reduce(powerX(q div 2,i,I)^2,I));}
512   return(reduce(var(i)*powerX(q-1,i,I),I));
513}
514example
515{ "EXAMPLE:"; echo = 2;
516   ring R = 0,(x,y),dp;
517   powerX(100,2,std(ideal(x3-1,y2-x)));
518}
519
520//======================================================================
521//=========================== Discrete Logarithm =======================
522//======================================================================
523
524//============== Shank's baby step - giant step ========================
525
526proc babyGiant(bigint b, bigint y, bigint p)
527"USAGE:  babyGiant(b,y,p);
528RETURN: the discrete logarithm x: b^x=y mod p
529NOTE:   This procedure works based on Shank's baby step - giant step method.
530EXAMPLE:example babyGiant; shows an example
531"
532{
533   int i,j,m;
534   list l;
535   bigint h=1;
536   bigint x;
537
538//choose m minimal such that m^2>p
539   for(i=1;i<=p;i++){if(i^2>p) break;}
540   m=i;
541
542//baby-step:  compute the table y*b^i for 1<=i<=m
543   for(i=1;i<=m;i++){l[i]=y*b^i mod p;}
544
545//giant-step: compute b^(m+j), 1<=j<=m and search in the baby-step table
546//for an i with y*b^i=b^(m*j). If found then x=m*j-i
547   bigint g=b^m mod p;
548   while(j<m)
549   {
550      j++;
551      h=h*g mod p;
552      for(i=1;i<=m;i++)
553      {
554         if(h==l[i])
555         {
556            x=m*j-i;
557            j=m;
558            break;
559         }
560     }
561   }
562   return(x);
563}
564example
565{ "EXAMPLE:"; echo = 2;
566   bigint b=2;
567   bigint y=10;
568   bigint p=101;
569   babyGiant(b,y,p);
570}
571
572//==============  Pollards rho  =================================
573
574proc rho(bigint b, bigint y, bigint p)
575"USAGE:  rho(b,y,p);
576RETURN: the discrete logarithm x=log_b(y): b^x=y mod p
577NOTE: Pollard's rho:
578       choose random f_0 in 0,...,p-2 ,e_0=0, define x_0=b^f_0, define
579       x_i=y^e_i*b^f_i as below. For i large enough there is i with
580       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,
581       d=gcd(s,p-1)=u*s+v*(p-1) then x=tu/d +j*(p-1)/d for some j (to be
582       found by trying)
583EXAMPLE:example rho; shows an example
584"
585{
586   int i=1;
587   int j;
588   bigint s,t;
589   list e,f,x;
590
591   e[1]=0;
592   f[1]=random(0,2147483629) mod (p-1);
593   x[1]=powerN(b,f[1],p);
594   while(i)
595   {
596      if((x[i] mod 3)==1)
597      {
598         x[i+1]=y*x[i] mod p;
599         e[i+1]=e[i]+1 mod (p-1);
600         f[i+1]=f[i];
601      }
602      if((x[i] mod 3)==2)
603      {
604         x[i+1]=x[i]^2 mod p;
605         e[i+1]=e[i]*2 mod (p-1);
606         f[i+1]=f[i]*2 mod (p-1);
607      }
608      if((x[i] mod 3)==0)
609      {
610         x[i+1]=x[i]*b mod p;
611         e[i+1]=e[i];
612         f[i+1]=f[i]+1 mod (p-1);
613      }
614      i++;
615      for(j=i-1;j>=1;j--)
616      {
617         if(x[i]==x[j])
618         {
619            s=(e[j]-e[i]) mod (p-1);
620            t=(f[i]-f[j]) mod (p-1);
621            if(s!=0)
622            {
623               i=0;
624            }
625            else
626            {
627               e[1]=0;
628               f[1]=random(0,2147483629) mod (p-1);
629               x[1]=powerN(b,f[1],p);
630               i=1;
631            }
632            break;
633         }
634      }
635   }
636
637   list w=extgcd(s,p-1);
638   bigint u=w[2];
639   bigint d=w[1];
640
641   bigint a=(t*u div d) mod (p-1);
642
643   bigint pn=powerN(b,a,p);
644   if (pn<0) { pn=pn+p;}
645   while(powerN(b,a,p)!=y)
646   {
647      a=(a+(p-1) div d) mod (p-1);
648      if (a<0) { a=a+p-1; }
649   }
650   return(a);
651}
652example
653{ "EXAMPLE:"; echo = 2;
654   bigint b=2;
655   bigint y=10;
656   bigint p=101;
657   rho(b,y,p);
658}
659//====================================================================
660//====================== Primality Tests =============================
661//====================================================================
662
663//================================= Miller-Rabin =====================
664
665proc MillerRabin(bigint n, int k)
666"USAGE:  MillerRabin(n,k);
667RETURN: 1 if n is prime, 0 else
668NOTE: probabilistic test of Miller-Rabin with k loops to test if n is prime.
669       Using the theorem: If n is prime, n-1=2^s*r, r odd, then
670       powerN(a,r,n)=1 or powerN(a,r*2^i,n)=-1 for some i
671EXAMPLE:example MillerRabin; shows an example
672"
673{
674   if(n<0){n=-n;}
675   if((n==2)||(n==3)){return(1);}
676   if((n mod 2)==0){return(0);}
677
678   int i;
679   bigint a,b,j,r,s;
680   r=n-1;
681   s=0;
682   while((r mod 2)==0)
683   {
684      s=s+1;
685      r=r div 2;
686   }
687   while(i<k)
688   {
689      i++;
690      a=random(2,2147483629) mod n; if(a==0){a=3;}
691      if(gcd(a,n)!=1){return(0);}
692      b=powerN(a,r,n);
693      if(b!=1)
694      {
695         j=0;
696         while(j<s)
697         {
698           if(((b+1) mod n)==0) break;
699           b=powerN(b,2,n);
700           j=j+1;
701         }
702         if(j==s){return(0);}
703      }
704   }
705   return(1);
706}
707example
708{ "EXAMPLE:"; echo = 2;
709   bigint x=2;
710   x=x^787-1;
711   MillerRabin(x,3);
712}
713
714//======================= Soloway-Strassen  ==========================
715
716proc SolowayStrassen(bigint n, int k)
717"USAGE:  SolowayStrassen(n,k);
718RETURN: 1 if n is prime, 0 else
719NOTE: probabilistic test of Soloway-Strassen with k loops to test if n is
720       prime using the theorem: If n is prime then
721       powerN(a,(n-1)/2,n)=Jacobi(a,n) mod n
722EXAMPLE:example SolowayStrassen; shows an example
723"
724{
725   if(n<0){n=-n;}
726   if((n==2)||(n==3)){return(1);}
727   if((n mod 2)==0){return(0);}
728
729   bigint a,pn,jn;
730   int i;
731   while(i<k)
732   {
733      i++;
734      a=random(2,2147483629) mod n; if(a==0){a=3;}
735      if(gcd(a,n)!=1){return(0);}
736      pn=powerN(a,(n-1) div 2,n);
737      if (pn<0) { pn=pn+n;}
738      jn=Jacobi(a,n) mod n;
739      if (jn<0) { jn=jn+n;}
740      if(pn!=jn){return(0);}
741   }
742   return(1);
743}
744example
745{ "EXAMPLE:"; echo = 2;
746   bigint h=10;
747   bigint p=h^100+267;
748   //p=h^100+43723;
749   //p=h^200+632347;
750   SolowayStrassen(h,3);
751}
752
753
754/*
755ring R=0,z,dp;
756number p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
757number q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
758number n=p*q;
759SolowayStrassen(n,3);
760*/
761
762//===================== Pocklington-Lehmer ==============================
763
764proc PocklingtonLehmer(bigint N, list #)
765"USAGE: PocklingtonLehmer(N); optional: PocklingtonLehmer(N,L);
766        L a list of the first k primes
767RETURN:message N is not prime or {A,{p},{a_p}} as certificate for N being prime
768NOTE:assumes that it is possible to factorize N-1=A*B such that gcd(A,B)=1,
769      the factorization of A is completely known and A^2>N .
770      N is prime if and only if for each prime factor p of A we can find
771      a_p such that a_p^(N-1)=1 mod N and gcd(a_p^((N-1)/p)-1,N)=1
772EXAMPLE:example PocklingtonLehmer; shows an example
773"
774{
775   bigint m=intRoot(N);
776   if(size(#)>0)
777   {
778      list S=PollardRho(N-1,10000,1,#);
779   }
780   else
781   {
782      list S=PollardRho(N-1,10000,1);
783   }
784   int i,j;
785   bigint A=1;
786   bigint p,a,g;
787   list PA;
788   list re;
789
790   while(i<size(S))
791   {
792      p=S[i+1];
793      A=A*p;
794      PA[i+1]=p;
795      if(A>m){break;}
796
797      while(1)
798      {
799        p=p*S[i+1];
800        if(((N-1) mod p)==0)
801        {
802           A=A*p;
803        }
804        else
805        {
806           break;
807        }
808      }
809      i++;
810   }
811   if(A<=m)
812   {
813     A=N div A;
814     PA=list(S[size(S)]);
815   }
816   for(i=1;i<=size(PA);i++)
817   {
818      a=1;
819      while(a<N-1)
820      {
821         a=a+1;
822         if(powerN(a,N-1,N)!=1){return("not prime");}
823         g=gcd(powerN(a,(N-1) div PA[i],N),N);
824         if(g==1)
825         {
826           re[size(re)+1]=list(PA[i],a);
827           break;
828         }
829         if(g<N){"not prime";return(g);}
830      }
831   }
832   return(list(A,re));
833}
834example
835{ "EXAMPLE:"; echo = 2;
836   bigint N=105554676553297;
837   PocklingtonLehmer(N);
838   list L=primList(1000);
839   PocklingtonLehmer(N,L);
840}
841
842//=======================================================================
843//======================= Factorization =================================
844//=======================================================================
845
846//======================= Pollards rho  =================================
847
848proc PollardRho(bigint n, int k, int allFactors, list #)
849"USAGE:  PollardRho(n,k,allFactors); optional: PollardRho(n,k,allFactors,L);
850         L a list of the first k primes
851RETURN: a list of factors of n (which could be just n),if allFactors=0@*
852        a list of all factors of n ,if allFactors=1
853NOTE: probabilistic rho-algorithm of Pollard to find a factor of n in k loops.
854      Creates a sequence x_i such that (x_i)^2=(x_2i)^2 mod n for some i,
855      computes gcd(x_i-x_2i,n) to find a divisor. To define the sequence
856      choose x,a and define x_n+1=x_n^2+a mod n, x_1=x.
857      If allFactors is 1, it tries to find recursively all prime factors
858      using the Soloway-Strassen test.
860EXAMPLE:example PollardRho; shows an example
861"
862{
863    int i,j;
864    list L=primList(100);
865    list re,se;
866    if(n<0){n=-n;}
867    if(n==1){return(re);}
868
869//this is optional: test whether a prime of the list # devides n
870    if(size(#)>0)
871    {
872       L=#;
873    }
874    for(i=1;i<=size(L);i++)
875    {
876       if((n mod L[i])==0)
877       {
878          re[size(re)+1]=L[i];
879          while((n mod L[i])==0)
880          {
881             n=n div L[i];
882          }
883       }
884       if(n==1){return(re);}
885    }
886    int e=size(re);
887//here the rho-algorithm starts
888    bigint a,d,x,y;
889    while(n>1)
890    {
891       a=random(2,2147483629);
892       x=random(2,2147483629);
893       y=x;
894       d=1;
895       i=0;
896       while(i<k)
897       {
898          i++;
899          x=powerN(x,2,n); x=(x+a) mod n;
900          y=powerN(y,2,n); y=(y+a) mod n;
901          y=powerN(y,2,n); y=(y+a) mod n;
902          d=gcd(x-y,n);
903          if(d>1)
904          {
905             re[size(re)+1]=d;
906             while((n mod d)==0)
907             {
908               n=n div d;
909             }
910             break;
911          }
912          if(i==k)
913          {
914             re[size(re)+1]=n;
915             n=1;
916          }
917       }
918    }
919    if(allFactors)      //want to obtain all prime factors
920    {
921       i=e;
922       while(i<size(re))
923       {
924          i++;
925
926          if(!SolowayStrassen(re[i],5))
927          {
928             se=PollardRho(re[i],2*k,1);
929             re[i]=se[size(se)];
930             for(j=1;j<=size(se)-1;j++)
931             {
932                re[size(re)+1]=se[j];
933             }
934             i--;
935          }
936       }
937    }
938    return(re);
939}
940example
941{ "EXAMPLE:"; echo = 2;
942   bigint h=10;
943   bigint p=h^30+4;
944   PollardRho(p,5000,0);
945}
946
947//======================== Pollards p-factorization ================
948proc pFactor(bigint n,int B, list P)
949"USAGE:  pFactor(n,B,P); n to be factorized, B a bound , P a list of primes
950RETURN: a list of factors of n or n if no factor found
951NOTE: Pollard's p-factorization
952       creates the product k of powers of primes (bounded by B)  from
953       the list P with the idea that for a prime divisor p of n we have
954       p-1|k, and then p divides gcd(a^k-1,n) for some random a
955EXAMPLE:example pFactor; shows an example
956"
957{
958   int i;
959   bigint k=1;
960   bigint w;
961   while(i<size(P))
962   {
963      i++;
964      w=P[i];
965      if(w>B) break;
966      while(w*P[i]<=B)
967      {
968         w=w*P[i];
969      }
970      k=k*w;
971   }
972   bigint a=random(2,2147483629);
973   bigint d=gcd(powerN(a,k,n)-1,n);
974   if((d>1)&&(d<n)){return(d);}
975   return(n);
976}
977example
978{ "EXAMPLE:"; echo = 2;
979   list L=primList(1000);
980   pFactor(1241143,13,L);
981   bigint h=10;
982   h=h^30+25;
983   pFactor(h,20,L);
984}
985
987
988proc quadraticSieve(bigint n, int c, list B, int k)
989"USAGE:  quadraticSieve(n,c,B,k); n to be factorized, [-c,c] the
990         sieve-intervall, B a list of primes,
991         k for using the first k elements in B
992RETURN: a list of factors of n or the message: no divisor found
993NOTE: The idea being used is to find x,y such that x^2=y^2 mod n then
994      gcd(x-y,n) can be a proper divisor of n
996"
997{
998   bigint f,d;
999   int i,j,l,s,p;
1000   list S,tmp;
1001   intvec v;
1002   v[k]=0;
1003
1004//compute the integral part of the square root of n
1005   bigint m=intRoot(n);
1006
1007//consider the function f(X)=(X+m)^2-n and compute for s in [-c,c] the values
1008   while(i<=2*c)
1009   {
1010      f=(i-c+m)^2-n;
1011      tmp[1]=i-c+m;
1012      tmp[2]=f;
1013      tmp[3]=v;
1014      S[i+1]=tmp;
1015      i++;
1016   }
1017
1018//the sieve with p in B
1019//find all s in [-c,c] such that f(s) has all prime divisors in the first
1020//k elements of B and the decomposition of f(s). They are characterized
1021//by 1 or -1 at the second place of S[j]:
1022//S[j]=j-c+m,f(j-c)/p_1^v_1*...*p_k^v_k, v_1,...,v_k maximal
1023   for(i=1;i<=k;i++)
1024   {
1025      p=B[i];
1026      if((p>2)&&(Jacobi(n,p)==-1)){i++;continue;}//n is no quadratic rest mod p
1027      j=1;
1028      while(j<=p)
1029      {
1030         if(j>2*c+1) break;
1031         f=S[j][2];
1032         v=S[j][3];
1033         s=0;
1034         while((f mod p)==0)
1035         {
1036            s++;
1037            f=f div p;
1038         }
1039         if(s)
1040         {
1041           S[j][2]=f;
1042           v[i]=s;
1043           S[j][3]=v;
1044           l=j;
1045           while(l+p<=2*c+1)
1046           {
1047              l=l+p;
1048              f=S[l][2];
1049              v=S[l][3];
1050              s=0;
1051              while((f mod p)==0)
1052              {
1053                s++;
1054                f=f div p;
1055              }
1056              S[l][2]=f;
1057              v[i]=s;
1058              S[l][3]=v;
1059           }
1060         }
1061         j++;
1062      }
1063   }
1064   list T;
1065   for(j=1;j<=2*c+1;j++)
1066   {
1067      if((S[j][2]==1)||(S[j][2]==-1))
1068      {
1069         T[size(T)+1]=S[j];
1070      }
1071   }
1072
1073//the system of equations for the exponents {l_s} for the f(s) such
1074//product f(s)^l_s is a square (l_s are 1 or 0)
1075   bigintmat M[k+1][size(T)];
1076   for(j=1;j<=size(T);j++)
1077   {
1078      if(T[j][2]==-1){M[1,j]=1;}
1079      for(i=1;i<=k;i++)
1080      {
1081         M[i+1,j]=T[j][3][i];
1082      }
1083   }
1084   intmat G=solutionsMod2(M);
1085
1086//construction of x and y such that x^2=y^2 mod n and d=gcd(x-y,n)
1087//y=square root of product f(s)^l_s
1088//x=product s+m
1089   bigint x=1;
1090   bigint y=1;
1091
1092   for(i=1;i<=ncols(G);i++)
1093   {
1094      kill v;
1095      intvec v;
1096      v[k]=0;
1097      for(j=1;j<=size(T);j++)
1098      {
1099         x=x*T[j][1]^G[j,i] mod n;
1100         if((T[j][2]==-1)&&(G[j,i]==1)){y=-y;}
1101         v=v+G[j,i]*T[j][3];
1102
1103      }
1104      for(l=1;l<=k;l++)
1105      {
1106         y=y*B[l]^(v[l] div 2) mod n;
1107      }
1108      d=gcd(x-y,n);
1109      if((d>1)&&(d<n)){return(d);}
1110   }
1111   return("no divisor found");
1112}
1113example
1114{ "EXAMPLE:"; echo = 2;
1115   list L=primList(5000);
1118}
1119
1120//======================================================================
1121//==================== elliptic curves  ================================
1122//======================================================================
1123
1124//================= elementary operations ==============================
1125
1126proc isOnCurve(bigint N, bigint a, bigint b, list P)
1127"USAGE:  isOnCurve(N,a,b,P);
1128RETURN: 1 or 0 (depending on whether P is on the curve or not)
1129NOTE: checks whether P=(P[1]:P[2]:P[3]) is a point on the elliptic
1130      curve defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1131EXAMPLE:example isOnCurve; shows an example
1132"
1133{
1134   if(((P[2]^2*P[3]-P[1]^3-a*P[1]*P[3]^2-b*P[3]^3) mod N)!=0){return(0);}
1135   return(1);
1136}
1137example
1138{ "EXAMPLE:"; echo = 2;
1139   isOnCurve(32003,5,7,list(10,16,1));
1140}
1141
1142proc ellipticAdd(bigint N, bigint a, bigint b, list P, list Q)
1144RETURN: list L, representing the point P+Q
1145NOTE: P=(P[1]:P[2]:P[3]), Q=(Q[1]:Q[2]:Q[3]) points on the elliptic curve
1146      defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1148"
1149{
1150   if(N==2){ERROR("not implemented for 2");}
1151   int i;
1152   for(i=1;i<=3;i++)
1153   {
1154      P[i]=P[i] mod N; if (P[i]<0) { P[i]=P[i]+N:}
1155      Q[i]=Q[i] mod N; if (Q[i]<0) { Q[i]=Q[i]+N;}
1156   }
1157   list Resu;
1158   Resu[1]=bigint(0);
1159   Resu[2]=bigint(1);
1160   Resu[3]=bigint(0);
1161   list Error;
1162   Error[1]=0;
1163   //test for ellictic curve
1164   bigint D=4*a^3+27*b^2;
1165   bigint g=gcd(D,N);
1166   if(g==N){return(Error);}
1167   if(g!=1)
1168   {
1169      P[4]=g;
1170      return(P);
1171   }
1172   if(((P[1]==0)&&(P[2]==0)&&(P[3]==0))||((Q[1]==0)&&(Q[2]==0)&&(Q[3]==0)))
1173   {
1174      Error[1]=-2;
1175      return(Error);
1176   }
1177   if(!isOnCurve(N,a,b,P)||!isOnCurve(N,a,b,Q))
1178   {
1179      Error[1]=-1;
1180      return(Error);
1181   }
1182   if(P[3]==0){return(Q);}
1183   if(Q[3]==0){return(P);}
1184   list I=extgcd(P[3],N);
1185   if(I[1]!=1)
1186   {
1187      P[4]=I[1];
1188      return(P);
1189   }
1190   P[1]=P[1]*I[2] mod N;
1191   P[2]=P[2]*I[2] mod N;
1192   I=extgcd(Q[3],N);
1193   if(I[1]!=1)
1194   {
1195      P[4]=I[1];
1196      return(P);
1197   }
1198   Q[1]=Q[1]*I[2] mod N;
1199   Q[2]=Q[2]*I[2] mod N;
1200   if((P[1]==Q[1])&&(((P[2]+Q[2]) mod N)==0)){return(Resu);}
1201   bigint L;
1202   if((P[1]==Q[1])&&(P[2]==Q[2]))
1203   {
1204      I=extgcd(2*Q[2],N);
1205      if(I[1]!=1)
1206      {
1207         P[4]=I[1];
1208         return(P);
1209      }
1210      L=I[2]*(3*Q[1]^2+a) mod N;
1211   }
1212   else
1213   {
1214      I=extgcd(Q[1]-P[1],N);
1215      if(I[1]!=1)
1216      {
1217         P[4]=I[1];
1218         return(P);
1219      }
1220      L=(Q[2]-P[2])*I[2] mod N;
1221   }
1222   Resu[1]=(L^2-P[1]-Q[1]) mod N;
1223   if (Resu[1]<0) { Resu[1]=Resu[1]+N; }
1224   Resu[2]=(L*(P[1]-Resu[1])-P[2]) mod N;
1225   if (Resu[2]<0) { Resu[2]=Resu[2]+N; }
1226   Resu[3]=bigint(1);
1227   return(Resu);
1228}
1229example
1230{ "EXAMPLE:"; echo = 2;
1231   bigint N=11;
1232   bigint a=1;
1233   bigint b=6;
1234   list P,Q;
1235   P[1]=2;
1236   P[2]=4;
1237   P[3]=1;
1238   Q[1]=3;
1239   Q[2]=5;
1240   Q[3]=1;
1242}
1243
1244proc ellipticMult(bigint N, bigint a, bigint b, list P, bigint k)
1245"USAGE:  ellipticMult(N,a,b,P,k);
1246RETURN: a list L representing the point k*P
1247NOTE:  P=(P[1]:P[2]:P[3]) a point on the elliptic curve defined by
1248       y^2z=x^3+a*xz^2+b*z^3  over Z/N
1249EXAMPLE:example ellipticMult; shows an example
1250"
1251{
1252   if(P[3]==0){return(P);}
1253   list resu;
1254   resu[1]=bigint(0);
1255   resu[2]=bigint(1);
1256   resu[3]=bigint(0);
1257
1258   if(k==0){return(resu);}
1259   if(k==1){return(P);}
1261   if(k==-1)
1262   {
1263      resu=P;
1264      resu[2]=N-P[2];
1265      return(resu);
1266   }
1267   if(k<0)
1268   {
1269      resu=ellipticMult(N,a,b,P,-k);
1270      return(ellipticMult(N,a,b,resu,-1));
1271   }
1272   if((k mod 2)==0)
1273   {
1274      resu=ellipticMult(N,a,b,P,k/2);
1276   }
1277   resu=ellipticMult(N,a,b,P,k-1);
1279}
1280example
1281{ "EXAMPLE:"; echo = 2;
1282   bigint N=11;
1283   bigint a=1;
1284   bigint b=6;
1285   list P;
1286   P[1]=2;
1287   P[2]=4;
1288   P[3]=1;
1289   ellipticMult(N,a,b,P,3);
1290}
1291
1292//================== Random for elliptic curves =====================
1293
1294proc ellipticRandomCurve(bigint N)
1295"USAGE:  ellipticRandomCurve(N);
1296RETURN: a list of two random numbers a,b and 4a^3+27b^2 mod N
1297NOTE:   y^2z=x^3+a*xz^2+b^2*z^3 defines an elliptic curve over Z/N
1298EXAMPLE:example ellipticRandomCurve; shows an example
1299"
1300{
1301   int k;
1302   while(k<=10)
1303   {
1304     k++;
1305     bigint a=random(1,2147483647) mod N;
1306     bigint b=random(1,2147483647) mod N;
1307     //test for ellictic curve
1308     bigint D=4*a^3+27*b^4; //the constant term is b^2
1309     bigint g=gcd(D,N);
1310     if(g<N){return(list(a,b,g));}
1311   }
1312   ERROR("no random curve found");
1313}
1314example
1315{ "EXAMPLE:"; echo = 2;
1316   ellipticRandomCurve(32003);
1317}
1318
1319proc ellipticRandomPoint(bigint N, bigint a, bigint b)
1320"USAGE:  ellipticRandomPoint(N,a,b);
1321RETURN: a list representing  a random point (x:y:z) of the elliptic curve
1322        defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
1323EXAMPLE:example ellipticRandomPoint; shows an example
1324"
1325{
1326   bigint x=random(1,2147483647) mod N;
1327   bigint h=x^3+a*x+b;
1328   h=h mod N;
1329   list resu;
1330   resu[1]=x;
1331   resu[2]=0;
1332   resu[3]=1;
1333   if(h==0){return(resu);}
1334
1335   bigint n=Jacobi(h,N);
1336   if(n==0)
1337   {
1338      resu=-5;
1339      "N is not prime";
1340      return(resu);
1341   }
1342   if(n==1)
1343   {
1344      resu[2]=squareRoot(h,N);
1345      return(resu);
1346   }
1347   return(ellipticRandomPoint(N,a,b));
1348}
1349example
1350{ "EXAMPLE:"; echo = 2;
1351   ellipticRandomPoint(32003,3,181);
1352}
1353
1354
1355
1356//====================================================================
1357//======== counting the points of an elliptic curve  =================
1358//====================================================================
1359
1360//==================   the trivial approaches  =======================
1361proc countPoints(bigint N, bigint a, bigint b)
1362"USAGE:  countPoints(N,a,b);
1363RETURN: the number of points of the elliptic curve defined by
1364        y^2=x^3+a*x+b  over Z/N
1365NOTE: trivial approach
1366EXAMPLE:example countPoints; shows an example
1367"
1368{
1369  bigint x;
1370  bigint r=N+1;
1371  while(x<N)
1372  {
1373     r=r+Jacobi((x^3+a*x+b) mod N,N);
1374     x=x+1;
1375  }
1376  return(r);
1377}
1378example
1379{ "EXAMPLE:"; echo = 2;
1380   countPoints(181,71,150);
1381}
1382
1383proc ellipticAllPoints(bigint N, bigint a, bigint b)
1384"USAGE:  ellipticAllPoints(N,a,b);
1385RETURN: list of points (x:y:z) of the elliptic curve defined by
1386        y^2z=x^3+a*xz^2+b*z^3  over Z/N
1387EXAMPLE:example ellipticAllPoints; shows an example
1388"
1389{
1390   list resu,point;
1391   point[1]=0;
1392   point[2]=1;
1393   point[3]=0;
1394   resu[1]=point;
1395   point[3]=1;
1396   bigint x,h,n;
1397   while(x<N)
1398   {
1399      h=(x^3+a*x+b) mod N;
1400      if(h==0)
1401      {
1402         point[1]=x;
1403         point[2]=0;
1404         resu[size(resu)+1]=point;
1405      }
1406      else
1407      {
1408         n=Jacobi(h,N);
1409         if(n==1)
1410         {
1411            n=squareRoot(h,N);
1412            point[1]=x;
1413            point[2]=n;
1414            resu[size(resu)+1]=point;
1415            point[2]=N-n;
1416            resu[size(resu)+1]=point;
1417         }
1418      }
1419      x=x+1;
1420   }
1421   return(resu);
1422}
1423example
1424{ "EXAMPLE:"; echo = 2;
1425   list L=ellipticAllPoints(181,71,150);
1426   size(L);
1427   L[size(L)];
1428}
1429
1430//================ the algorithm of Shanks and Mestre =================
1431
1432proc ShanksMestre(bigint q, bigint a, bigint b, list #)
1433"USAGE:  ShanksMestre(q,a,b); optional: ShanksMestre(q,a,b,s); s the number
1434         of loops in the algorithm (default s=1)
1435RETURN: the number of points of the elliptic curve defined by
1436         y^2=x^3+a*x+b  over Z/N
1437NOTE: algorithm of Shanks and Mestre (baby-step-giant-step)
1438EXAMPLE:example ShanksMestre; shows an example
1439"
1440{
1441   bigint n=intRoot(4*q);
1442   bigint m=intRoot(intRoot(16*q))+1;
1443   bigint d;
1444   int i,j,k,s;
1445   list B,K,T,P,Q,R,mP;
1446   B[1]=list(0,1,0);
1447   if(size(#)>0)
1448   {
1449      s=#[1];
1450   }
1451   else
1452   {
1453      s=1;
1454   }
1455   while(k<s)
1456   {
1457      P =ellipticRandomPoint(q,a,b);
1458      Q =ellipticMult(q,a,b,P,n+q+1);
1459
1460      while(j<m)
1461      {
1462         j++;
1464      }
1466      mP[2]=q-mP[2];
1467      while(i<m)                            //giant-step
1468      {
1469         j=0;
1470         while(j<m)
1471         {
1472            j++;
1473            if((Q[1]==B[j][1])&&(Q[2]==B[j][2])&&(Q[3]==B[j][3]))
1474            {
1475               T[1]=P;
1476               T[2]=q+1+n-(i*m+j-1);
1477               K[size(K)+1]=T;
1478               if(size(K)>1)
1479               {
1480                  if(K[size(K)][2]!=K[size(K)-1][2])
1481                  {
1482                     d=gcd(K[size(K)][2],K[size(K)-1][2]);
1483                     if(ellipticMult(q,a,b,K[size(K)],d)[3]==0)
1484                     {
1485                       K[size(K)][2]=K[size(K)-1][2];
1486                     }
1487                  }
1488               }
1489               i=int(m);
1490               break;
1491            }
1492         }
1493         i=i+1;
1495      }
1496      k++;
1497   }
1498   if(size(K)>0)
1499   {
1500      int te=1;
1501      for(i=1;i<=size(K)-1;i++)
1502      {
1503         if(K[size(K)][2]!=K[i][2])
1504         {
1505           if(ellipticMult(q,a,b,K[i],K[size(K)][2])[3]!=0)
1506           {
1507             te=0;
1508             break;
1509           }
1510         }
1511      }
1512      if(te)
1513      {
1514         return(K[size(K)][2]);
1515      }
1516   }
1517   return(ShanksMestre(q,a,b,s));
1518}
1519example
1520{ "EXAMPLE:"; echo = 2;
1521   ShanksMestre(32003,71,602);
1522}
1523
1524//==================== Schoof's algorithm =============================
1525
1526proc Schoof(bigint N,bigint a, bigint b)
1527"USAGE:  Schoof(N,a,b);
1528RETURN: the number of points of the elliptic curve defined by
1529        y^2=x^3+a*x+b  over Z/N
1530NOTE:  algorithm of Schoof
1531EXAMPLE:example Schoof; shows an example
1532"
1533{
1534   int pr=printlevel;
1535//test for ellictic curve
1536   bigint D=4*a^3+27*b^2;
1537   bigint G=gcd(D,N);
1538   if(G==N){ERROR("not an elliptic curve");}
1539   if(G!=1){ERROR("not a prime");}
1540
1541//=== small N
1542  // if((N<=500)&&(pr<5)){return(countPoints(int(N),a,b));}
1543
1544//=== the general case
1545   bigint q=intRoot(4*N);
1546   list L=primL(2*q);
1547   int r=size(L);
1548   list T;
1549   int i,j;
1550   for(j=1;j<=r;j++)
1551   {
1552      T[j]=(testElliptic(int(N),a,b,L[j])+int(q)) mod L[j];
1553   }
1554   if(pr>=5)
1555   {
1556      "===================================================================";
1557      "Chinese remainder :";
1558      for(i=1;i<=size(T);i++)
1559      {
1560         " x =",T[i]," mod ",L[i];
1561      }
1562      "gives t+ integral part of the square root of q (to be positive)";
1563      chineseRem(T,L);
1564      "we obtain t = ",chineseRem(T,L)-q;
1565      "===================================================================";
1566   }
1567   bigint t=chineseRem(T,L)-q;
1568   return(N+1-t);
1569}
1570example
1571{ "EXAMPLE:"; echo = 2;
1572   Schoof(32003,71,602);
1573}
1574
1575/*
1576needs 518 sec
1577Schoof(2147483629,17,3567);
15782147168895
1579*/
1580
1581
1582proc generateG(number a,number b, int m)
1583"USAGE:  generateG(a,b,m);
1584RETURN: m-th division polynomial
1585NOTE: generate the so-called division polynomials, i.e., the recursively defined
1586polynomials p_m=generateG(a,b,m) in Z[x, y] such that, for a point (x:y:1) on the
1587elliptic curve defined by y^2=x^3+a*x+b  over Z/N the point@*
1588m*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)
1589and m*P=0 if and only if p_m(P)=0
1590EXAMPLE:example generateG; shows an example
1591"
1592{
1593   if(m==0){return(poly(0));}
1594   if(m==1){return(poly(1));}
1595   if(m==2){return(2*var(1));}
1596   if(m==3){return(3*var(2)^4+6*a*var(2)^2+12*b*var(2)-a^2);}
1597   if(m==4)
1598   {
1599      return(4*var(1)*(var(2)^6+5*a*var(2)^4+20*b*var(2)^3-5*a^2*var(2)^2
1600        -4*a*b*var(2)-8*b^2-a^3));
1601   }
1602   if((m mod 2)==0)
1603   {
1604      return((generateG(a,b,m div 2+2)*generateG(a,b,m div 2-1)^2
1605        -generateG(a,b,m div 2-2)*generateG(a,b,m div 2+1)^2)
1606        *generateG(a,b,m div 2)/(2*var(1)));
1607   }
1608   return(generateG(a,b,(m-1) div 2+2)*generateG(a,b,(m-1) div 2)^3
1609     -generateG(a,b,(m-1) div 2-1)*generateG(a,b,(m-1) div 2+1)^3);
1610}
1611example
1612{ "EXAMPLE:"; echo = 2;
1613   ring R = 0,(x,y),dp;
1614   generateG(7,15,4);
1615}
1616
1617
1618static proc testElliptic(int q,bigint aa,bigint bb,int l)
1619"USAGE:  testElliptic(q,a,b,l);
1620RETURN: an integer t, the trace of the Frobenius
1621NOTE: the kernel for the Schoof algorithm: looks for the t such that for all
1622      points (x:y:1) in C[l]={P in C | l*P=0},C the elliptic curve defined by
1623      y^2=x^3+a*x+b  over Z/q with group structure induced by 0=(0:1:0),
1624      (x:y:1)^(q^2)-t*(x:y:1)^q -ql*(x:y:1)=(0:1:0), ql= q mod l, trace of
1625      Frobenius.
1626EXAMPLE:example testElliptic; shows an example
1627"
1628{
1629   int pr=printlevel;
1630   ring S=q,(y,x),(L(100000),lp);
1631   number a=aa;
1632   number b=bb;
1633   poly F=y2-x3-a*x-b;       // the curve C
1634   poly G=generateG(a,b,l);
1635   ideal I=std(ideal(F,G));  // the points C[l]
1636   poly xq=powerX(q,2,I);
1637   poly yq=powerX(q,1,I);
1638   poly xq2=reduce(subst(xq,x,xq,y,yq),I);
1639   poly yq2=reduce(subst(yq,x,xq,y,yq),I);
1640   ideal J;
1641   int ql=q mod l;
1642   if(ql==0){ERROR("q is not prime");}
1643   int t;
1644   poly F1,F2,G1,G2,P1,P2,Q1,Q2,H1,H2,L1,L2;
1645
1646   if(pr>=5)
1647   {
1648      "===================================================================";
1649      "q=",q;
1650      "l=",l;
1651      "q mod l=",ql;
1652      "the Groebner basis for C[l]:";I;
1653      "x^q mod I = ",xq;
1654      "x^(q^2) mod I = ",xq2;
1655      "y^q mod I = ",yq;
1656      "y^(q^2) mod I = ",yq2;
1657      pause();
1658   }
1659   //==== l=2 =============================================================
1660   if(l==2)
1661   {
1662      xq=powerX(q,2,std(x3+a*x+b));
1663      J=std(ideal(xq-x,x3+a*x+b));
1664      if(deg(J[1])==0){t=1;}
1665      if(pr>=5)
1666      {
1667         "===================================================================";
1668         "the case l=2";
1669         "the gcd(x^q-x,x^3+ax+b)=",J[1];
1670         pause();
1671      }
1672      return(t);
1673   }
1674   //=== (F1/G1,F2/G2)=[ql](x,y) ==========================================
1675   if(ql==1)
1676   {
1677      F1=x;G1=1;F2=y;G2=1;
1678   }
1679   else
1680   {
1681      G1=reduce(generateG(a,b,ql)^2,I);
1682      F1=reduce(x*G1-generateG(a,b,ql-1)*generateG(a,b,ql+1),I);
1683      G2=reduce(4*y*generateG(a,b,ql)^3,I);
1684      F2=reduce(generateG(a,b,ql+2)*generateG(a,b,ql-1)^2
1685               -generateG(a,b,ql-2)*generateG(a,b,ql+1)^2,I);
1686
1687   }
1688   if(pr>=5)
1689   {
1690      "===================================================================";
1691      "the point ql*(x,y)=(F1/G1,F2/G2)";
1692      "F1=",F1;
1693      "G1=",G1;
1694      "F2=",F2;
1695      "G2=",G2;
1696      pause();
1697   }
1698   //==== the case t=0 :  the equations for (x,y)^(q^2)=-[ql](x,y) ===
1699   J[1]=xq2*G1-F1;
1700   J[2]=yq2*G2+F2;
1701   if(pr>=5)
1702   {
1703      "===================================================================";
1704      "the case t=0 mod l";
1705      "the equations for (x,y)^(q^2)=-[ql](x,y) :";
1706      J;
1707      "the test, if they vanish for all points in C[l]:";
1708      reduce(J,I);
1709      pause();
1710   }
1711   //=== test if all points of C[l] satisfy  (x,y)^(q^2)=-[ql](x,y)
1712   //=== if so: t mod l =0 is returned
1713   if(size(reduce(J,I))==0){return(0);}
1714
1715   //==== test for (x,y)^(q^2)=[ql](x,y) for some point
1716
1717   J=xq2*G1-F1,yq2*G2-F2;
1718   J=std(J+I);
1719   if(pr>=5)
1720   {
1721      "===================================================================";
1722      "test if (x,y)^(q^2)=[ql](x,y) for one point";
1723      "if so, the Frobenius has an eigenvalue 2ql/t: (x,y)^q=(2ql/t)*(x,y)";
1724      "it follows that t^2=4q mod l";
1725      "if w is one square root of q mod l";
1726      "t =2w mod l or -2w mod l ";
1727      "-------------------------------------------------------------------";
1728      "the equations for (x,y)^(q^2)=[ql](x,y) :";
1729      xq2*G1-F1,yq2*G2-F2;
1730      "the test if one point satisfies them";
1731      J;
1732      pause();
1733   }
1734   if(deg(J[1])>0)
1735   {
1736      int w=int(squareRoot(q,l));
1737      //=== +/-2w mod l zurueckgeben, wenn (x,y)^q=+/-[w](x,y)
1738      //==== the case t>0 :  (Q1/P1,Q2/P2)=[w](x,y) ==============
1739      if(w==1)
1740      {
1741         Q1=x;P1=1;Q2=y;P2=1;
1742      }
1743      else
1744      {
1745         P1=reduce(generateG(a,b,w)^2,I);
1746         Q1=reduce(x*G1-generateG(a,b,w-1)*generateG(a,b,w+1),I);
1747         P2=reduce(4*y*generateG(a,b,w)^3,I);
1748         Q2=reduce(generateG(a,b,w+2)*generateG(a,b,w-1)^2
1749               -generateG(a,b,w-2)*generateG(a,b,w+1)^2,I);
1750      }
1751      J=xq*P1-Q1,yq*P2-Q2;
1752      J=std(I+J);
1753      if(pr>=5)
1754      {
1755      "===================================================================";
1756      "the Frobenius has an eigenvalue, one of the roots of  w^2=q mod l:";
1757      "one root is:";w;
1758      "test, if it is the eigenvalue (if not it must be -w):";
1759      "the equations for (x,y)^q=w*(x,y)";I;xq*P1-Q1,yq*P2-Q2;
1760      "the Groebner basis";
1761       J;
1762         pause();
1763      }
1764      if(deg(J[1])>0){return(2*w mod l);}
1765      return(-2*w mod l);
1766   }
1767
1768   //==== the case t>0 :  (Q1/P1,Q2/P2)=(x,y)^(q^2)+[ql](x,y) =====
1769   P1=reduce(G1*G2^2*(F1-xq2*G1)^2,I);
1770   Q1=reduce((F2-yq2*G2)^2*G1^3-F1*G2^2*(F1-xq2*G1)^2-xq2*P1,I);
1771   P2=reduce(P1*G2*(F1-xq2*G1),I);
1772   Q2=reduce((xq2*P1-Q1)*(F2-yq2*G2)*G1-yq2*P2,I);
1773
1774   if(pr>=5)
1775   {
1776      "we are in the general case:";
1777      "(x,y)^(q^2)!=ql*(x,y) and (x,y)^(q^2)!=-ql*(x,y) ";
1778      "the point (Q1/P1,Q2/P2)=(x,y)^(q^2)+[ql](x,y)";
1779      "Q1=",Q1;
1780      "P1=",P1;
1781      "Q2=",Q2;
1782      "P2=",P2;
1783      pause();
1784   }
1785   while(t<(l-1) div 2)
1786   {
1787      t++;
1788      //====  (H1/L1,H2/L2)=[t](x,y)^q ===============================
1789      if(t==1)
1790      {
1791         H1=xq;L1=1;
1792         H2=yq;L2=1;
1793      }
1794      else
1795      {
1796         H1=x*generateG(a,b,t)^2-generateG(a,b,t-1)*generateG(a,b,t+1);
1797         H1=subst(H1,x,xq,y,yq);
1798         H1=reduce(H1,I);
1799         L1=generateG(a,b,t)^2;
1800         L1=subst(L1,x,xq,y,yq);
1801         L1=reduce(L1,I);
1802         H2=generateG(a,b,t+2)*generateG(a,b,t-1)^2
1803           -generateG(a,b,t-2)*generateG(a,b,t+1)^2;
1804         H2=subst(H2,x,xq,y,yq);
1805         H2=reduce(H2,I);
1806         L2=4*y*generateG(a,b,t)^3;
1807         L2=subst(L2,x,xq,y,yq);
1808         L2=reduce(L2,I);
1809      }
1810      J=Q1*L1-P1*H1,Q2*L2-P2*H2;
1811      if(pr>=5)
1812      {
1813      "we test now the different t, 0<t<=(l-1)/2:";
1814      "the point (H1/L1,H2/L2)=[t](x,y)^q :";
1815      "H1=",H1;
1816      "L1=",L1;
1817      "H2=",H2;
1818      "L2=",L2;
1819      "the equations for (x,y)^(q^2)+[ql](x,y)=[t](x,y)^q :";J;
1820      "the test";reduce(J,I);
1821      "the test for l-t (the x-cordinate is the same):";
1822       Q1*L1-P1*H1,Q2*L2+P2*H2;
1823       reduce(ideal(Q1*L1-P1*H1,Q2*L2+P2*H2),I);
1824         pause();
1825      }
1826      if(size(reduce(J,I))==0){return(t);}
1827      J=Q1*L1-P1*H1,Q2*L2+P2*H2;
1828      if(size(reduce(J,I))==0){return(l-t);}
1829   }
1830   ERROR("something is wrong in testElliptic");
1831}
1832example
1833{ "EXAMPLE:"; echo = 2;
1834   testElliptic(1267985441,338474977,64740730,3);
1835}
1836
1837//============================================================================
1838//================== Factorization and Primality Test ========================
1839//============================================================================
1840
1841//============= Lenstra's ECM Factorization ==================================
1842
1843proc factorLenstraECM(bigint N, list S, int B, list #)
1844"USAGE:  factorLenstraECM(N,S,B); optional: factorLenstraECM(N,S,B,d);
1845         d+1 the number of loops in the algorithm (default d=0)
1846RETURN: a factor of N or the message no factor found
1847NOTE: - computes a factor of N using Lenstra's ECM factorization@*
1848      - the idea is that the fact that N is not prime is dedected using
1849        the operations on the elliptic curve
1850      - is similarly to Pollard's p-1-factorization
1851EXAMPLE:example factorLenstraECM; shows an example
1852"
1853{
1854   list L,P;
1855   bigint g,M,w;
1856   int i,j,k,d;
1857   int l=size(S);
1858   if(size(#)>0)
1859   {
1860      d=#[1];
1861   }
1862
1863   while(i<=d)
1864   {
1865      L=ellipticRandomCurve(N);
1866      if(L[3]>1){return(L[3]);} //the discriminant was not invertible
1867      P=list(0,L[2],1);
1868      j=0;
1869      M=1;
1870      while(j<l)
1871      {
1872         j++;
1873         w=S[j];
1874         if(w>B) break;
1875         while(w*S[j]<B)
1876         {
1877           w=w*S[j];
1878         }
1879         M=M*w;
1880         P=ellipticMult(N,L[1],L[2]^2,P,w);
1881         if(size(P)==4){return(P[4]);}  //some inverse did not exsist
1882         if(P[3]==0){break;}            //the case M*P=0
1883      }
1884      i++;
1885   }
1886   return("no factor found");
1887}
1888example
1889{ "EXAMPLE:"; echo = 2;
1890   list L=primList(1000);
1891   factorLenstraECM(181*32003,L,10,5);
1892   bigint h=10;
1893   h=h^30+25;
1894   factorLenstraECM(h,L,4,3);
1895}
1896
1897//================= ECPP (Goldwasser-Kilian) a primaly-test =============
1898
1899proc ECPP(bigint N)
1900"USAGE:  ECPP(N);
1901RETURN: message:N is not prime or {L,P,m,q} as certificate for N being prime@*
1902         L a list (y^2=x^3+L[1]*x+L[2] defines an elliptic curve C)@*
1903         P a list ((P[1]:P[2]:P[3]) is a point of C)@*
1904         m,q integers
1905ASSUME: gcd(N,6)=1
1906NOTE:   The basis of the algorithm is the following theorem:
1907         Given C, an elliptic curve over Z/N, P a point of C(Z/N),
1908         m an integer, q a prime with the following properties:
1909         - q|m
1910         - q>(4-th root(N) +1)^2
1911         - m*P=0=(0:1:0)
1912         - (m/q)*P=(x:y:z) and z a unit in Z/N
1913         Then N is prime.
1914EXAMPLE:example ECPP; shows an example
1915"
1916{
1917   list L,S,P;
1918   bigint m,q;
1919   int i;
1920
1921   bigint n=intRoot(intRoot(N));
1922   n=(n+1)^2;                         //lower bound for q
1923   while(1)
1924   {
1925      L=ellipticRandomCurve(N);       //a random elliptic curve C
1926      m=ShanksMestre(N,L[1],L[2],3);  //number of points of the curve C
1927      S=PollardRho(m,10000,1);        //factorization of m
1928      for(i=1;i<=size(S);i++)         //search for q between the primes
1929      {
1930         q=S[i];
1931         if(n<q){break;}
1932      }
1933      if(n<q){break;}
1934   }
1935   bigint u=m/q;
1936   while(1)
1937   {
1938      P=ellipticRandomPoint(N,L[1],L[2]);  //a random point on C
1939      "P=",P;
1940      if(ellipticMult(N,L[1],L[2],P,m)[3]!=0){"N is not prime";return(-5);}
1941      if(ellipticMult(N,L[1],L[2],P,u)[3]!=0)
1942      {
1943         L=delete(L,3);
1944         return(list(L,P,m,q));
1945      }
1946   }
1947}
1948example
1949{ "EXAMPLE:"; echo = 2;
1950   bigint N=1267985441;
1951   ECPP(N);
1952}
1953
1954static proc wordToNumber(string s)
1955{
1956   int i;
1957   intvec v;
1958   bigint n;
1959   bigint t=27;
1960   for(i=size(s);i>0;i--)
1961   {
1962      if(s[i]=="a"){v[i]=0;}
1963      if(s[i]=="b"){v[i]=1;}
1964      if(s[i]=="c"){v[i]=2;}
1965      if(s[i]=="d"){v[i]=3;}
1966      if(s[i]=="e"){v[i]=4;}
1967      if(s[i]=="f"){v[i]=5;}
1968      if(s[i]=="g"){v[i]=6;}
1969      if(s[i]=="h"){v[i]=7;}
1970      if(s[i]=="i"){v[i]=8;}
1971      if(s[i]=="j"){v[i]=9;}
1972      if(s[i]=="k"){v[i]=10;}
1973      if(s[i]=="l"){v[i]=11;}
1974      if(s[i]=="m"){v[i]=12;}
1975      if(s[i]=="n"){v[i]=13;}
1976      if(s[i]=="o"){v[i]=14;}
1977      if(s[i]=="p"){v[i]=15;}
1978      if(s[i]=="q"){v[i]=16;}
1979      if(s[i]=="r"){v[i]=17;}
1980      if(s[i]=="s"){v[i]=18;}
1981      if(s[i]=="t"){v[i]=19;}
1982      if(s[i]=="u"){v[i]=20;}
1983      if(s[i]=="v"){v[i]=21;}
1984      if(s[i]=="w"){v[i]=22;}
1985      if(s[i]=="x"){v[i]=23;}
1986      if(s[i]=="y"){v[i]=24;}
1987      if(s[i]=="z"){v[i]=25;}
1988      if(s[i]==" "){v[i]=26;}
1989   }
1990   for(i=1;i<=size(s);i++)
1991   {
1992      n=n+v[i]*t^(i-1);
1993   }
1994   return(n);
1995}
1996
1997static proc numberToWord(bigint n)
1998{
1999   int i,j;
2000   string v;
2001   list s;
2002   bigint t=27;
2003   bigint mm;
2004   bigint nn=n;
2005   while(nn>t)
2006   {
2007      j++;
2008      mm=nn mod t;
2009      s[j]=mm;
2010      nn=(nn-mm) div t;
2011   }
2012   j++;
2013   s[j]=nn;
2014   for(i=1;i<=j;i++)
2015   {
2016      if(s[i]==0){v=v+"a";}
2017      if(s[i]==1){v=v+"b";}
2018      if(s[i]==2){v=v+"c";}
2019      if(s[i]==3){v=v+"d";}
2020      if(s[i]==4){v=v+"e";}
2021      if(s[i]==5){v=v+"f";}
2022      if(s[i]==6){v=v+"g";}
2023      if(s[i]==7){v=v+"h";}
2024      if(s[i]==8){v=v+"i";}
2025      if(s[i]==9){v=v+"j";}
2026      if(s[i]==10){v=v+"k";}
2027      if(s[i]==11){v=v+"l";}
2028      if(s[i]==12){v=v+"m";}
2029      if(s[i]==13){v=v+"n";}
2030      if(s[i]==14){v=v+"o";}
2031      if(s[i]==15){v=v+"p";}
2032      if(s[i]==16){v=v+"q";}
2033      if(s[i]==17){v=v+"r";}
2034      if(s[i]==18){v=v+"s";}
2035      if(s[i]==19){v=v+"t";}
2036      if(s[i]==20){v=v+"u";}
2037      if(s[i]==21){v=v+"v";}
2038      if(s[i]==22){v=v+"w";}
2039      if(s[i]==23){v=v+"x";}
2040      if(s[i]==24){v=v+"y";}
2041      if(s[i]==25){v=v+"z";}
2042      if(s[i]==26){v=v+" ";}
2043   }
2044   return(v);
2045}
2046
2047proc code(string s)
2048"USAGE:  code(s); s a string
2049ASSUME:  s contains only small letters and space
2050COMPUTE: a bigint, RSA-coding of the string s
2051RETURN:  return RSA-coding of the string s as string
2052EXAMPLE: code;  shows an example
2053"
2054{
2055   ring r=0,x,dp;
2056   bigint
2057p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
2058   bigint
2059q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
2060   bigint n=p*q;
2061   bigint phi=(p-1)*(q-1);
2062   bigint e=1234567891;
2063   //bigint d=extgcd(e,phi)[2];
2064   bigint m=wordToNumber(s);
2065   bigint c=powerN(m,e,n);
2066   string cc=string(c);
2067   return(cc);
2068}
2069example
2070{"EXAMPLE:";  echo = 2;
2071  string s="i go to school";
2072  code(s);
2073}
2074
2075proc decodeString(string g)
2076"USAGE:  decodeString(s); s a string
2077ASSUME:  s is a string of a bigint, the output of code
2078COMPUTE: a string, RSA-decoding of the string s
2079RETURN:  return RSA-decoding of the string s as string
2080EXAMPLE: decodeString;  shows an example
2081"
2082{
2083   bigint
2084p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
2085   bigint
2086q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
2087   bigint n=p*q;
2088   bigint phi=(p-1)*(q-1);
2089   bigint e=1234567891;
2090   bigint d=extgcd(e,phi)[2];
2091   execute("bigint c="+g+";");
2092   bigint f=powerN(c,d,n);
2093   string s=numberToWord(f);
2094   return(s);
2095}
2096example
2097{"EXAMPLE:";  echo = 2;
2098  string
2099s="78638618599886548153321853785991541374544958648147340831959482696082179852616053583234149080198937632782579537867262780982185252913122030800897193851413140758915381848932565";
2100  string t=decodeString(s);
2101  t;
2102}
2103/*----------------------------------------------------------------------------
2104 * set stuff
2105 * -------------------------------------------------------------------------*/
2106static proc set_multiply_list_content(list h)
2107"USAGE:   set_multiply_list_content(h)
2108RETURN:  An integer c als product of all elements in h
2109EXAMPLE: example set_multiply_list_content; shows an example;
2110"
2111{
2112  int c = 1;
2113  for (int i=1;i<=size(h);i++)
2114  {
2115    c = c*h[i];
2116  }
2117  return(c);
2118}
2119example
2120{
2121  "EXAMPLE:"; echo = 2;
2122  list h=2,4,5;
2123  set_multiply_list_content(h);
2124}
2125
2126static proc set_delete_certain_element(list a, int e)
2127"USAGE:   set_delete_certain_element(a,e)
2128RETURN:  A list a without element e. If e was not in the list before, a will not be changed
2129EXAMPLE: example set_delete_certain_element; shows an example.
2130"
2131{
2132  list output_list = a;
2133  for (int i=1;i<=size(a);i++)
2134  {
2135    if (a[i]==e)
2136    {
2137      output_list = delete(output_list,i);
2138    }
2139  }
2140  return(output_list);
2141}
2142example
2143{
2144  "EXAMPLE:"; echo = 2;
2145  list h=2,4,5;
2146  set_delete_certain_element(h,4);
2147  set_delete_certain_element(h,10);
2148}
2149
2150static proc set_bubblesort_int(list output_list)
2151"USAGE:   set_bubblesort_int(output_list)
2152RETURN:  An ascending sorted list with integer values.
2153EXAMPLE: set_bubblesort_int; shows an example.
2154"
2155{
2156  output_list = bubblesort(output_list);
2157  //Cast every value into an integer
2158  for (int i=1; i<=size(output_list);i++)
2159  {
2160    output_list[i] = int(output_list[i]);
2161  }
2162  return(output_list);
2163}
2164example
2165{
2166  "EXAMPLE:"; echo = 2;
2167  list output_list=10,4,5,24,9;
2168  set_bubblesort_int(output_list);
2169}
2170
2171static proc set_is_set(list a)
2172"USAGE:   set_is_set(a)
2173RETURN:  1 if the list is a set, 0 the list contains any duplicated elements
2174EXAMPLE: set_is_set; shows an example.
2175"
2176{
2177  int i,v;
2178  for (v=1; v<=size(a); v++)
2179  {
2180    for (i=v+1; i<=size(a); i++)
2181    {
2182      if (a[i]==a[v])
2183      {
2184        return(0);
2185      }
2186    }
2187  }
2188  return(1);
2189}
2190example
2191{
2192  "EXAMPLE:"; echo = 2;
2193  list set = 1,5,10,2;
2194  list noset = 1,5,10,2,5;
2195  set_is_set(set);
2196  set_is_set(noset);
2197}
2198
2199static proc set_contains(list a, int e)
2200"USAGE:   set_contains(a,e)
2201RETURN:  1 if the list contains e, 0 otherwise
2202EXAMPLE: set_contains; shows an example.
2203"
2204{
2205  for (int v=1; v<=size(a); v++)
2206  {
2207    if (a[v]==e)
2208    {
2209      return(1);
2210    }
2211  }
2212  return(0);
2213}
2214example
2215{
2216  "EXAMPLE:"; echo = 2;
2217  list set = 1,5,10,2;
2218  set_contains(set,5);
2219  set_contains(set,40);
2220}
2221
2222static proc set_delete_duplicates(list a)
2223"USAGE:   set_delete_duplicates(a)
2224RETURN:  a list a without any duplicated elements
2225EXAMPLE: set_delete_duplicates; shows an example.
2226"
2227{
2228  int i;
2229  list output_list = a[1];
2230  for (i=1; i<=size(a); i++)
2231  {
2232    if (set_contains(output_list,a[i])==0)
2233    {
2234      output_list = insert(output_list,a[i]);
2235    }
2236  }
2237  return(output_list);
2238}
2239example
2240{
2241  "EXAMPLE:"; echo = 2;
2242  list set = 1,5,10,2,10,5;
2243  set_delete_duplicates(set);
2244}
2245
2246static proc set_equals(list a,list b)
2247"USAGE:   set_equals(a, b)
2248RETURN:  1 if the lists are equal from a set-structure standpoint, 0 otherwise
2249EXAMPLE: set_equals; shows an example.
2250"
2251{
2252  //Checks if the lists have the same length
2253  if (size(a)!=size(b))
2254  {
2255    return(0);
2256  }
2257
2258  //Sorts the lists
2259  a = set_bubblesort_int(a);
2260  b = set_bubblesort_int(b);
2261
2262  //Checks every single element of both lists
2263  for (int i=1; i<=size(a); i++)
2264  {
2265    if (a[i]!=b[i])
2266    {
2267      return(0);
2268    }
2269  }
2270  return(1);
2271}
2272example
2273{
2274  "EXAMPLE:"; echo = 2;
2275  list set1 = 1,5,10,2;
2276  list set2 = 10,2,5,1;
2277  list set3 = 1,5,9,2;
2278  set_equals(set1,set2);
2279  set_equals(set1,set3);
2280}
2281
2282static proc set_insert(list a, int e)
2283"USAGE:   set_insert(a,e)
2284RETURN:  list a containing element e
2285EXAMPLE: set_insert; shows an example.
2286"
2287{
2288  if(set_contains(a,e))
2289  {
2290    return(a);
2291  }
2292  else
2293  {
2294    a=insert(a,e);
2295    return(a);
2296  }
2297}
2298example
2299{
2300  "EXAMPLE:"; echo = 2;
2301  list set = 1,5,10,2;
2302  set_insert(set,5);
2303  set_insert(set,22);
2304}
2305
2306static proc set_union(list a, list b)
2307"USAGE:   set_union(a, b)
2308RETURN:  list a as union of a and b
2309EXAMPLE: set_union; shows an example.
2310"
2311{
2312  for (int i=1; i<=size(b); i++)
2313  {
2314    if (set_contains(a,b[i])==0)
2315    {
2316      a = insert(a,b[i]);
2317    }
2318  }
2319  return(a);
2320}
2321example
2322{
2323  "EXAMPLE:"; echo = 2;
2324  list set1 = 1,5,10,2;
2325  list set2 = 5,10,93,58,29;
2326  set_union(set1,set2);
2327}
2328
2329static proc set_section(list a, list b)
2330"USAGE:   set_section(a, b)
2331RETURN:  list output_list as intersection of a and b
2332EXAMPLE: set_section; shows an example.
2333"
2334{
2335  list output_list;
2336  for (int i=1; i<=size(a); i++)
2337  {
2338    if (set_contains(b,a[i])==1)
2339    {
2340      output_list = insert(output_list,a[i]);
2341    }
2342  }
2343  return(output_list);
2344}
2345example
2346{
2347  "EXAMPLE:"; echo = 2;
2348  list set1 = 1,5,10,2;
2349  list set2 = 5,10,93,58,29;
2350  set_section(set1,set2);
2351}
2352
2353static proc set_list_delete_duplicates(list a)
2354"USAGE:   set_list_delete_duplicates(a)
2355RETURN:  list output_list with no duplicated lists
2356EXAMPLE: set_list_delete_duplicates; shows an example.
2357"
2358{
2359  int v;
2360  int i;
2361  int counter;
2362  int out_size;
2363  list output_list=insert(output_list,a[1]);
2364  //Create a new list and try to insert every list element from a into that list. If a list is already inserted into the
2365  //new list, do nothing.
2366  for (i=2; i<=size(a); i++)
2367  {
2368    out_size = size(output_list);
2369    counter = 0;
2370
2371    for (v=1; v<=out_size; v++)
2372    {
2373      if (set_equals(output_list[v],a[i])==1)
2374      {
2375        counter++;
2376      }
2377    }
2378    if (counter==0)
2379    {
2380      output_list = insert(output_list,a[i]);
2381    }
2382  }
2383  return(output_list);
2384}
2385example
2386{
2387  "EXAMPLE:"; echo = 2;
2388  list set1 = 1,5,10,2;
2389  list set2 = 1,10,2,5;
2390  list set3 = 1,2,3;
2391  list superset = set1,set2,set3;
2392  set_list_delete_duplicates(superset);
2393}
2394
2395static proc set_construct_increasing_set(int maxelement)
2396"USAGE:   set_construct_increasing_set(maxelement)
2397RETURN:  list output_list with increasing elements from 1 to maxelement
2398EXAMPLE: set_construct_increasing_set; shows an example.
2399"
2400{
2401  list output_list;
2402  for (int i=1; i<=maxelement; i++)
2403  {
2404    output_list = insert(output_list, i);
2405  }
2406  return(output_list);
2407}
2408example
2409{
2410  "EXAMPLE:"; echo = 2;
2411  set_construct_increasing_set(10);
2412}
2413
2414static proc set_addtoall(list a, int element)
2416RETURN:  Transformed list with e_i+element for every element e_i in list a
2418"
2419{
2420  list output_list;
2421  int c;
2422  for (int i=1; i<=size(a); i++)
2423  {
2424    c = a[i]+element;
2425    output_list = insert(output_list,c);
2426  }
2427  return(set_turn(output_list));
2428}
2429example
2430{
2431  "EXAMPLE:"; echo = 2;
2432  list set = 1,5,10,2;
2434}
2435
2436static proc set_turn(list a)
2437"USAGE:   set_turn(a)
2438RETURN:  Turned list a
2439EXAMPLE: set_turn; shows an example.
2440"
2441{
2442  list output_list;
2443  for (int i=1; i<=size(a); i++)
2444  {
2445    output_list[size(a)+1-i] = a[i];
2446  }
2447  return(output_list);
2448}
2449example
2450{
2451  "EXAMPLE:"; echo = 2;
2452  list set = 1,5,10;
2453  set_turn(set);
2454}
2455
2456static proc set_subset_set(list a)
2457"USAGE:   set_subset_set(a)
2458RETURN:  Set of subsets
2459EXAMPLE: set_subset_set; shows an example.
2460"
2461{
2462  int v;
2463  int choice;
2464  list output_list;
2465  list start = a[1];
2466  list insertion_list;
2467  int output_length;
2468  output_list = insert(output_list,start);
2469
2470  for (int i=2; i<=size(a); i++)
2471  {
2472    choice = 0;
2473    start = a[i];
2474    output_list = insert(output_list,start);
2475    output_length = size(output_list);
2476    for (v=2; v<=output_length; v++)
2477    {
2478      insertion_list = set_insert(output_list[v],a[i]);
2479      output_list = insert(output_list, insertion_list,size(output_list));
2480    }
2481  }
2482  return(output_list);
2483}
2484example
2485{
2486  "EXAMPLE:"; echo = 2;
2487  list set = 1,5,10;
2488  set_subset_set(set);
2489}
2490/* -----------------------------------------------------------------
2491 * knapsack_utilities: Utility functions needed for knapsack
2492 * -----------------------------------------------------------------*/
2493proc calculate_ordering(bigint num1, bigint primitive, bigint mod1)
2494"USAGE:   calculate_ordering(num1, primitive, mod1)
2495RETURN:  x so that primitive^x == num1 mod mod1
2496EXAMPLE: example calculate_ordering; shows an example;
2497"
2498{
2499  for (int i=1;i<=int((mod1-2));i++)
2500  {
2501    if ((primitive^i%mod1)==num1)
2502    {
2503      return(i);
2504    }
2505  }
2506  return(0);
2507}
2508example
2509{
2510  "EXAMPLE:"; echo = 2;
2511  bigint mod1 = 33;
2512  bigint primitive = 14;
2513  bigint num1 = 5;
2514  calculate_ordering(num1,primitive,mod1);
2515}
2516
2517proc is_primitive_root(bigint primitive, bigint mod1)
2518"USAGE:   is_primitive_root(primitive, mod1)
2519RETURN:  1 if primitive is a primitive root modulo mod1, 0 otherwise
2520EXAMPLE: example is_primitive_root; shows an example;
2521"
2522{
2523  list output_list;
2524  for (int i=1;i<=int((mod1-1));i++)
2525  {
2526    output_list = set_insert(output_list,int((primitive^i%mod1)));
2527  }
2528  if (bigint(size(output_list))==bigint(mod1-1))
2529  {
2530    return(1);
2531  }
2532  else
2533  {
2534    return(0);
2535  }
2536}
2537example
2538{
2539  "EXAMPLE:"; echo = 2;
2540  is_primitive_root(3,7);
2541  is_primitive_root(2,7);
2542}
2543
2544proc find_first_primitive_root(bigint mod1)
2545"USAGE:   find_first_primitive_root(mod1)
2546RETURN:  First primitive root modulo mod1, 0 if no root can be found.
2547EXAMPLE: example find_first_primitive_root; shows an example;
2548"
2549{
2550  for (int i=0;i<=int(mod1-1);i++)
2551  {
2552    if (is_primitive_root(bigint(i),mod1)==1)
2553    {
2554      return(i);
2555    }
2556  }
2557  return(0);
2558}
2559example
2560{
2561  "EXAMPLE:"; echo = 2;
2562  ring r = 0,x,lp;
2563  find_first_primitive_root(7);
2564  find_first_primitive_root(557);
2565}
2566
2569RETURN:  binary encoded list, increased by 1
2570EXAMPLE: example binary_add; shows an example;
2571"
2572{
2573  int residual=1;
2574  int position = size(binary_list);
2575  while((residual==1)&&(position!=0))
2576  {
2577    if(binary_list[position]==0)
2578    {
2579      binary_list[position]=1;
2580      residual=0;
2581    }
2582    else
2583    {
2584      binary_list[position]=0;
2585      position--;
2586    }
2587  }
2588  if (position==0)
2589  {
2590    binary_list = insert(binary_list,1);
2591  }
2592  return(binary_list);
2593}
2594example
2595{
2596  "EXAMPLE:"; echo = 2;
2597  ring r = 0,x,lp;
2598  list binary_list = 1,0,1,1,1;
2600}
2601
2602proc inverse_modulus(int num, int mod1)
2603"USAGE:   inverse_modulus(num, mod1)
2604RETURN:  inverse element of num modulo mod1
2605EXAMPLE: example inverse_modulus; shows an example;
2606"
2607{
2608  if (num>=mod1)
2609  {
2610    return(0);
2611  }
2612  else
2613  {
2614    for (int i=1;i<mod1;i++)
2615    {
2616      if ((i*num%mod1)==1)
2617      {
2618        return(i);
2619      }
2620    }
2621  }
2622}
2623example
2624{
2625  "EXAMPLE:"; echo = 2;
2626  ring r = 0,x,lp;
2627  int mod1 = 13;
2628  int num = 5;
2629  inverse_modulus(num,mod1);
2630}
2631
2632proc is_prime(int n)
2633"USAGE:   is_prime(n)
2634RETURN:  1 if n is prime, 0 otherwise
2635EXAMPLE: example is_prime; shows an example;
2636"
2637{
2638  int prime1 = 1;
2639  for (int i=n-1;i>1;i--)
2640  {
2641    if(n%i==0)
2642    {
2643        prime1 = 0;
2644    }
2645  }
2646  return(prime1);
2647}
2648example
2649{
2650  "EXAMPLE:"; echo = 2;
2651  ring r = 0,x,lp;
2652  is_prime(10);
2653  is_prime(7);
2654}
2655
2656proc find_biggest_index(list a)
2657"USAGE:   find_biggest_index( a)
2658RETURN:  Returns the index of the biggest element of a
2659EXAMPLE: example find_biggest_index; shows an example;
2660"
2661{
2662  list sortedlist = bubblesort(a);
2663  return(find_index(a,sortedlist[1]));
2664}
2665
2666proc find_index(list a, bigint e)
2667"USAGE:   find_index(a, e)
2668RETURN:  Returns the list index of element e in list a. Returns 0 if e is not in a
2669EXAMPLE: example find_index; shows an example;
2670"
2671{
2672  for(int i=1;i<=size(a);i++)
2673  {
2674    if (bigint(a[i])==e)
2675    {
2676      return(i);
2677    }
2678  }
2679  return(0);
2680}
2681example
2682{
2683  "EXAMPLE:"; echo = 2;
2684  list a = 1,5,20,6,37;
2685  find_index(a,20);
2686  find_index(a,6);
2687  find_index(a,100);
2688}
2689/* ------------------------------------------------------------------
2690 * Knapsack Algorithmus such as solving several knapsack problems,
2691 * kryptographic algorithms and algorithms for creatings suitable knapsacks
2692 * ---------------------------------------------------------------- */
2693proc subset_sum01(list knapsack, int solution)
2694"USAGE:   subset_sum01(knapsack,solution)
2695RETURN:  binary list of the positions of the elements included in the subset sum or 0 if no solution exists
2696NOTE: This will return the first solution of the ssk-problem, given be the smallest binary encoding. It wont return several solutions if they exist
2697EXAMPLE: example subset_sum01; shows an example;
2698"
2699{
2700  int i;
2701  int v;
2702  int comparable;
2703  //Check if the knapsack is a set
2704  if (set_is_set(knapsack)==1)
2705  {
2706    //Create a binary list full of zeroes
2707    list binary_list;
2708    for (i=1;i<=size(knapsack);i++)
2709    {
2710      binary_list = insert(binary_list,0);
2711    }
2713
2714    for(i=1;i<=2^(size(knapsack));i++)
2715    {
2716      comparable = 0;
2717      //Create the Subset-Sum for the actual binary coding of binary_list
2718      for (v=1;v<=size(knapsack);v++)
2719      {
2720        comparable = comparable+knapsack[v]*binary_list[v];
2721      }
2722      //Check if the sum equals the solution
2723      if (comparable==solution)
2724      {
2725        return(binary_list);
2726      }
2727      else
2728      {
2730      }
2731    }
2732    return(0);
2733  }
2734  else
2735  {
2736    return(0);
2737  }
2738}
2739example
2740{
2741  "EXAMPLE:"; echo = 2;
2742  list h=1,4,7,32;
2743  subset_sum01(h,20);
2744  subset_sum01(h,11);
2745  subset_sum01(h,33);
2746}
2747
2748proc subset_sum02(list knapsack, int sol)
2749"USAGE:   subset_sum02(knapsack,sol)
2750RETURN:  binary list of the positions of the elements included in the subset sum or 0 if no solution exists
2751EXAMPLE: example subset_sum02; shows an example;
2752"
2753{
2754  list summands;
2755  int calcu;
2756  int i;
2757  //Create a sorted copy of the knapsack, calling it worksack
2758  list worksack = set_bubblesort_int(knapsack);
2759  int counter = 1;
2760  while((counter<=size(worksack))&&(sol>0))
2761  {
2762    //Try to substract an element of the knapsack from the capacity. Create a list with all the summands used
2763    calcu = sol-worksack[counter];
2764    if (calcu>=0)
2765    {
2766      sol = sol-worksack[counter];
2767      summands = insert(summands,int(worksack[counter]));
2768    }
2769    counter++;
2770  }
2771  if(sol>0)
2772  {
2773    return(0);
2774  }
2775
2776  //Get the index of the summands of the original knapsack and change the bits in the binary list wich will be the solution
2777  list binary_list;
2778  for (i=1;i<=size(knapsack);i++)
2779  {
2780    binary_list = insert(binary_list,0);
2781  }
2782
2783  for (i=1; i<=size(knapsack);i++)
2784  {
2785    if (set_contains(summands,knapsack[i])==1)
2786    {
2787      binary_list[i]=1;
2788    }
2789  }
2790  return(binary_list);
2791}
2792example
2793{
2794  "EXAMPLE:"; echo = 2;
2795  list h=1,4,7,32;
2796  subset_sum02(h,20);
2797  subset_sum02(h,11);
2798  subset_sum02(h,33);
2799}
2800
2801proc unbounded_knapsack(list knapsack, list profit, int capacity)
2802"USAGE:   unbounded_knapsack(knapsack,profit,capacity)
2803RETURN:  list of maximum profit of each iteration. For example, output_list[2] contains the maximum profit that can be achieved if the knapsack has capacity 2.
2804EXAMPLE: unbounded_knapsack; shows an example;
2805"
2806{
2807  int i;
2808  int v;
2809  list output_list;
2810  for (i=1;i<=capacity+1;i++)
2811  {
2812    output_list = insert(output_list,0);
2813  }
2814  for (i=1;i<=capacity+1;i++)
2815  {
2816    for(v=1;v<=size(knapsack);v++)
2817    {
2818      if (knapsack[v]<i)
2819      {
2820        if(output_list[i]<(output_list[i-knapsack[v]]+profit[v]))
2821        {
2822          output_list[i] = output_list[i-knapsack[v]]+profit[v];
2823        }
2824      }
2825    }
2826  }
2827  return(output_list);
2828}
2829example
2830{
2831  "EXAMPLE:"; echo = 2;
2832  list h=1,4,7,32;
2833  list knapsack = 5,2;
2834  list profit = 10,3;
2835  int capacity = 5;
2836  unbounded_knapsack(knapsack,profit,capacity);
2837}
2838
2839proc multidimensional_knapsack(matrix m, list capacities, list profits)
2840"USAGE:   multidimensional_knapsack(m,capacities,profits)
2841RETURN:  binary list of the positions of the elements included in the optimal selection
2842EXAMPLE: example multidimensional_knapsack; shows an example;
2843"
2844{
2845  int index;
2846  list output_list;
2847  list nolist;
2848  list y_list;
2849  list minmax_list;
2850  int i;
2851  int v;
2852  int checkint;
2853  //create the output_list full of zeroes with the length of all given selections
2854  for (i=1;i<=size(profits);i++)
2855  {
2856    output_list = insert(output_list,0);
2857  }
2858
2859  //Create the List E with all indices of the output_list that haven't been used yet.
2860  list E = set_turn(set_construct_increasing_set(size(profits)));
2861
2862  //Repeat till every index in E is used.
2863  while(size(E)>0)
2864  {
2865    y_list = nolist;
2866    for (i=1;i<=size(E);i++)
2867    {
2868      //Create all possible elements of y_i (y_list). minimax_list will be replaced of an empty list in each iteration
2869      minmax_list = nolist;
2870      for (v=1; v<=size(capacities);v++)
2871      {
2872        if (set_contains(E,v)==1)
2873        {
2874        minmax_list = insert(minmax_list, bigint(capacities[v])/bigint(m[v,E[i]]));
2875        }
2876      }
2877      //Sort the elements so that it is easy to pick the smallest one
2878      minmax_list = bubblesort(minmax_list);
2879
2880      //insert Element y_i into y_list, wich is the smallest of (b_i/w_ij) like in the PECH algorithm description.
2881      y_list = insert(y_list,minmax_list[size(minmax_list)],size(y_list));
2882
2883
2884    }
2885
2886    //Check if all y_i in y_list are smaller than 1. If so, every additional selection will exceed the capacity and the algorithm stops.
2887    checkint=0;
2888    for(i=1;i<=size(y_list);i++)
2889    {
2890      if (y_list[i]>=1)
2891      {
2892        checkint=1;
2893      }
2894    }
2895    if (checkint==0)
2896    {
2897      return(output_list);
2898    }
2899
2900
2901    //Find the index of the selection and update the binary output_list
2902    minmax_list = nolist;
2903    for (i=1;i<=size(E);i++)
2904    {
2905      minmax_list = insert(minmax_list, profits[E[i]]*y_list[i],size(minmax_list));
2906    }
2907    index = find_biggest_index(minmax_list);
2908
2909    output_list[E[index]]=1;
2910
2911    //Update the capacities by substracting the weights of the selection
2912    for (i=1;i<=size(capacities);i++)
2913    {
2914      capacities[i] = capacities[i]- m[i,E[index]];
2915    }
2916    E = set_delete_certain_element(E,index);
2917
2918  }
2919  return(output_list);
2920}
2921example
2922{
2923  "EXAMPLE:"; echo = 2;
2924  ring r = 0,x,lp;
2925  matrix m[3][3] = 1,4,10,7,8,3,1,9,7;
2926  list c = 12,17,10;
2927  list p = 3,2,5;
2928  multidimensional_knapsack(m,c,p);
2929}
2930
2933RETURN:  a hard knapsack list
2934EXAMPLE: example naccache_stern_generation; shows an example;
2935"
2936{
2937  //Check if primenum is a prime and the gcd-Condition holds
2939  {
2940    return(0);
2941  }
2942  else
2943  {
2944    int i;
2945    int p;
2946    list primelist;
2947    int primecounter=2;
2948    //Generate the knapsack containing the smallest prime numbers so that primenum exceeds the product of all of them
2950    {
2951      if (is_prime(primecounter)==1)
2952      {
2953      primelist = insert(primelist,primecounter);
2954      }
2955      primecounter++;
2956    }
2957    primelist = delete(primelist,1);
2958
2959
2960    //Generate the hard knapsack of the length of the prime numbers containing zeroes.
2961    list hardknapsack;
2962    for (i=1;i<=size(primelist);i++)
2963    {
2964      hardknapsack = insert(hardknapsack,0);
2965    }
2966
2967    //Create the elements of the hard knapsack
2968    primecounter = 1;
2969    bigint calcu;
2970    i=0;
2971    while (i<size(primelist))
2972    {
2973      //Create some v_i^key%primenum and store it in calcu
2975
2976      //If calcu is one of the prime numbers in the knapsack, find the index and insert v_i in the hard knapsack at that given index
2977      if(set_contains(primelist,int(calcu))==1)
2978      {
2979        p=find_index(primelist,int(calcu));
2980        hardknapsack[p] = primecounter;
2981        i++;
2982      }
2983      primecounter++;
2984    }
2985  return(hardknapsack);
2986  }
2987}
2988example
2989{
2990  "EXAMPLE:"; echo = 2;
2991  naccache_stern_generation(5,292);
2992  naccache_stern_generation(5,293);
2993}
2994
2995proc naccache_stern_encryption(list knapsack, list message, int primenum)
2997RETURN:  an encrypted message as integer
2998EXAMPLE: example naccache_stern_encryption; shows an example;
2999"
3000{
3001  bigint solution = 1;
3002  if (size(knapsack)==size(message))
3003  {
3004    for(int i=1;i<=size(knapsack);i++)
3005    {
3007    }
3008    return(solution);
3009  }
3010  else
3011  {
3012    return(0);
3013  }
3014
3015}
3016example
3017{
3018  "EXAMPLE:"; echo = 2;
3019  //Please note that the values for primenum and hardknapsack have been obtained from the example of naccache_stern_generation!
3020  list hardknapsack = 85,164,117,44;
3022  list message = 1,0,1,0;
3024}
3025
3026proc naccache_stern_decryption(list knapsack, int key, int primenum, int message)
3028RETURN:  decrypted binary list
3029EXAMPLE: example naccache_stern_decryption; shows an example;
3030"
3031{
3032  //create a binary list with the length of the knapsack containing zeros
3034  int i;
3035  list binary_list;
3036  for (i=1;i<=size(knapsack);i++)
3037  {
3038    binary_list = insert(binary_list,0);
3039  }
3040
3041  //create primelist like in int naccache_stern_generation process
3042  list primelist;
3043  int primecounter=2;
3044  while(size(primelist)<size(knapsack))
3045  {
3046    if (is_prime(primecounter)==1)
3047    {
3048      primelist = insert(primelist,primecounter);
3049    }
3050      primecounter++;
3051  }
3052
3053  //find divisors of k and update the binarylist in a way that the positions of the divisors in primelist are marked
3054  for (i=1;i<=size(primelist);i++)
3055  {
3056    if(k%primelist[i]==0)
3057    {
3058      binary_list[i]=1;
3059    }
3060  }
3061  return(binary_list);
3062
3063}
3064example
3065{
3066  "EXAMPLE:"; echo = 2;
3067  //Please note that the values have been obtained from the example of naccache_stern_generation and naccache_stern_encryption!
3069  int message = 9945;
3070  int key = 5;
3071  list hardknapsack = 85,164,117,44;
3073}
3074
3075proc m_merkle_hellman_transformation(list knapsack, int primitive, int mod1)
3076"USAGE:   m_merkle_hellman_transformation(knapsack, primitive, mod1)
3077RETURN:  list containing a hard knapsack
3078EXAMPLE: example m_merkle_hellman_transformation; shows an example;
3079"
3080{
3081  bigint new_element;
3082  list output_list;
3083  //calculate the primitiv root of every element in knapsack and insert it into a new knapsack
3084  for (int i=size(knapsack);i>=1;i--)
3085  {
3086    new_element = calculate_ordering(knapsack[i],primitive,mod1);
3087    output_list = insert(output_list,int(new_element));
3088  }
3089  return(output_list);
3090}
3091example
3092{
3093  "EXAMPLE:"; echo = 2;
3094  //Please note that the values for primenum and hardknapsack have been obtained from the example of naccache_stern_generation and naccache_stern_encryption!
3095  list knapsack = 2,3,5,7;
3096  int mod1 = 211;
3097  int primitive = 2;
3098  m_merkle_hellman_transformation(knapsack,primitive,mod1);
3099}
3100
3101proc m_merkle_hellman_encryption(list knapsack, list message)
3102"USAGE:    m_merkle_hellman_encryption(knapsack, message)
3103RETURN:  an encrypted message as integer
3104NOTE: This works in the same way as merkle_hellman_encryption. The additional function is created to keep consistency with the needed functions for every kryptosystem.
3105EXAMPLE: example m_merkle_hellman_encryption; shows an example;
3106"
3107{
3108  return(merkle_hellman_encryption(knapsack,message));
3109}
3110example
3111{
3112  "EXAMPLE:"; echo = 2;
3113  //Please note that the values for primenum and hardknapsack have been obtained from the example of m_merkle_hellman_transformation!
3114  list knapsack = 1,43,132,139;
3115  list message = 1,0,0,1;
3116  m_merkle_hellman_encryption(knapsack,message);
3117}
3118
3119proc m_merkle_hellman_decryption(list knapsack, bigint primitive, bigint mod1, int message)
3120"USAGE:  m_merkle_hellman_decryption(knapsack, primitive, mod1, message)
3121RETURN:  decrypted binary list
3122EXAMPLE: example merkle_hellman_decryption; shows an example;
3123"
3124{
3125  //Convert message
3126  int factorizing = int((primitive^message)%mod1);
3127  int i;
3128
3129  //Create binary list of length of the knapsack, containing zeroes.
3130  list binary_list;
3131  for (i=1;i<=size(knapsack);i++)
3132  {
3133    binary_list = insert(binary_list,0);
3134  }
3135
3136  //factorize the converted message, mark the factor positions in knapsack as bits in binary_list
3137  for (i=1;i<=size(knapsack);i++)
3138  {
3139    if(factorizing%knapsack[i]==0)
3140    {
3141      binary_list[i]=1;
3142    }
3143  }
3144  return(binary_list);
3145}
3146example
3147{
3148  "EXAMPLE:"; echo = 2;
3149  //Please note that the values  have been obtained from the example of m_merkle_hellman_encryption and m_merkle_hellman_transformation!
3150  list knapsack = 2,3,5,7;
3151  int message = 140;
3152  bigint primitive = 2;
3153  bigint mod1 = 211;
3154  m_merkle_hellman_decryption(knapsack,primitive,mod1,message);
3155}
3156
3157proc merkle_hellman_transformation(list knapsack, int key, int mod1)
3158"USAGE:   merkle_hellman_transformation(knapsack, key, mod1)
3159RETURN:  hard knapsack
3160EXAMPLE: example merkle_hellman_transformation; shows an example;
3161"
3162{
3163  list output_list;
3164  int new_element;
3165  //transform every element in the knapsack with normal strong modular multiplication
3166  for (int i=size(knapsack);i>=1;i--)
3167  {
3168    new_element=knapsack[i]*key%mod1;
3169    output_list = insert(output_list,new_element);
3170  }
3171  return(output_list);
3172}
3173example
3174{
3175  "EXAMPLE:"; echo = 2;
3176  list knapsack = 1,3,5,12;
3177  int key = 3;
3178  int mod1 = 23;
3179  merkle_hellman_transformation(knapsack,key,mod1);
3180}
3181
3182proc merkle_hellman_encryption(list knapsack, list message)
3183"USAGE:   merkle_hellman_encryption(knapsack, message)
3184RETURN:  encrypted integer
3185EXAMPLE: example merkle_hellman_encryption; shows an example;
3186"
3187{
3188  int solution = 0;
3189  if (size(knapsack)!=size(message)||(set_is_set(knapsack)==0))
3190  {
3191    return(0);
3192  }
3193  else
3194  {
3195    for (int i=1;i<=size(knapsack);i++)
3196    {
3197      solution = solution+knapsack[i]*message[i];
3198    }
3199    return(solution);
3200  }
3201}
3202example
3203{
3204  "EXAMPLE:"; echo = 2;
3205  //Please note that the values  have been obtained from the example of merkle_hellman_transformation!
3206  list hardknapsack =3,9,15,13;
3207  list message = 0,1,0,1;
3208  merkle_hellman_encryption(hardknapsack,message);
3209}
3210
3211proc merkle_hellman_decryption(list knapsack, int key, int mod1, int message)
3212"USAGE:   merkle_hellman_decryption(knapsack, key, mod1, message)
3213RETURN:  decrypted binary list
3214EXAMPLE: example merkle_hellman_decryption; shows an example;
3215"
3216{
3217  int new_element;
3218  int t = inverse_modulus(key,mod1);
3219  int transformed_message;
3220  list binary_list;
3221  if ((set_is_set(knapsack)==1)&&(key<mod1))
3222  {
3223    //reconstruct easy knapsack be multiplying with the inverse modulus t
3224    list easy_knapsack;
3225    for (int i=size(knapsack);i>=1;i--)
3226    {
3227      new_element=knapsack[i]*t%mod1;
3228      easy_knapsack = insert(easy_knapsack,new_element);
3229    }
3230
3231    //solve the easy knapsack problem with subset_sum01 or subset_sum02
3232    transformed_message = (message*t)%mod1;
3233    transformed_message;
3234    binary_list = subset_sum01(easy_knapsack,transformed_message);
3235    return(binary_list)
3236  }
3237  else
3238  {
3239    return(0)
3240  }
3241}
3242example
3243{
3244  "EXAMPLE:"; echo = 2;
3245  //Please note that the values  have been obtained from the example of merkle_hellman_decryption and merkle_hellman_transformation!
3246  list hardknapsack =3,9,15,13;
3247  int key = 3;
3248  int message = 22;
3249  int mod1 = 23;
3250  merkle_hellman_decryption(hardknapsack, key, mod1, message);
3251}
3252
3253proc super_increasing_knapsack(int ksize)
3254"USAGE:   super_increasing_knapsack(ksize)
3255RETURN:  super-increasing knapsack list
3256EXAMPLE: super_increasing_knapsack; shows an example;
3257"
3258{
3259  list output_list = insert(output_list,1);
3260  int next_element;
3261
3262  for (int i=2; i<=ksize; i++)
3263  {
3264    next_element = calculate_max_sum(output_list)+1;
3265    output_list = insert(output_list,next_element);
3266  }
3267  return(output_list);
3268}
3269example
3270{
3271  "EXAMPLE:"; echo = 2;
3272  super_increasing_knapsack(10);
3273}
3274
3275proc h_increasing_knapsack(int ksize, int h)
3276"USAGE:   h_increasing_knapsack(ksize, h)
3277RETURN:  h-increasing knapsack list
3278EXAMPLE: h_increasing_knapsack; shows an example;
3279"
3280{
3281  int v;
3282  if (ksize<=h+1)
3283  {
3284    return(set_turn(super_increasing_knapsack(ksize)))
3285  }
3286  else
3287  {
3288    list out = set_turn(super_increasing_knapsack(h+1));
3289    int next_element;
3290    for (int i=h+2; i<=ksize; i++)
3291    {
3292      next_element = 0;
3293      for (v=i-h; v<=i-1; v++)
3294      {
3295        next_element = next_element+out[v];
3296      }
3297      next_element++;
3298      out = insert(out,next_element,size(out));
3299    }
3300    return(out);
3301  }
3302}
3303example
3304{
3305  "EXAMPLE:"; echo = 2;
3306  h_increasing_knapsack(10,5);
3307}
3308
3309proc injective_knapsack(int ksize, int kmaxelement)
3310"USAGE:   injective_knapsack(ksize, kmaxelement)
3311RETURN:  list of injective knapsacks with maximal element kmaxelement and size ksize
3312EXAMPLE: injective_knapsack; shows an example;
3313"
3314{
3315  //Create a List of size ksize with the greatest possible elements keeping the set structure
3316  list list_of_lists;
3317  list A = insert(A,kmaxelement);
3318  int i;
3319  for (i=2;i<=ksize;i++)
3320  {
3321    A = insert(A,kmaxelement-(i-1));
3322  }
3323  A = set_turn(A);
3324  list_of_lists = insert(list_of_lists,A);
3325
3326  //Create all possible sets containing the possible elements of A
3327  int residual;
3328  int position;
3329  while(A[1]==kmaxelement)
3330  {
3331    residual=3;
3332    position = ksize;
3333    while((residual!=0))
3334    {
3335      if(A[position]==1)
3336      {
3337        A[position]=kmaxelement-position+1;
3338        residual=1;
3339        position--;
3340      }
3341      else
3342      {
3343        A[position]=A[position]-1;
3344        residual=0;
3345      }
3346    }
3347    //Insert the list into the overall list if its a set
3348    if (set_is_set(A)==1)
3349    {
3350      list_of_lists = insert(list_of_lists,A);
3351    }
3352  }
3353  //delete the first element since it is smaller than kmaxelement
3354  list_of_lists = delete(list_of_lists,1);
3355  //delete duplicates
3356  list_of_lists = set_list_delete_duplicates(list_of_lists);
3357
3358  //Check if the remaining knapsacks are injective
3359  list output_list;
3360  for(i=1;i<=size(list_of_lists);i++)
3361  {
3362    if (is_injective(list_of_lists[i])==1)
3363    {
3364      output_list=insert(output_list,list_of_lists[i]);
3365    }
3366  }
3367  return(output_list);
3368}
3369example
3370{
3371  "EXAMPLE:"; echo = 2;
3372  injective_knapsack(3,9);
3373}
3374
3375proc calculate_max_sum(list a)
3376"USAGE:   calculate_max_sum(a)
3377RETURN:  sum of all elements in a
3378EXAMPLE: calculate_max_sum; shows an example;
3379"
3380{
3381  int sum = a[1];
3382  for (int i=2; i<=size(a);i++)
3383  {
3384    sum = sum+a[i];
3385  }
3386  return(sum);
3387}
3388example
3389{
3390  "EXAMPLE:"; echo = 2;
3391  list a = 1,5,3,2,12;
3392  calculate_max_sum(a);
3393}
3394
3395proc is_injective(list a)
3396"USAGE:   is_injective(a)
3397RETURN:  1 if a is injective, 0 otherwise
3398EXAMPLE: is_injective; shows an example;
3399"
3400{
3401  //Create all subsets of the set a
3402  list subsum = set_subset_set(a);
3403  list checklist=calculate_max_sum(subsum[1]);
3404  int calculator;
3405  for (int i=2; i<=size(subsum);i++)
3406  {
3407    //calculate the maximal subset_sum for every subset. Check if there are duplicated subset_sums. If so, a is not injective
3408    calculator = calculate_max_sum(subsum[i]);
3409    if (set_contains(checklist, calculator))
3410    {
3411      return(0);
3412    }
3413    else
3414    {
3415      checklist = insert(checklist,calculator);
3416    }
3417  }
3418  return(1);
3419}
3420example
3421{
3422  "EXAMPLE:"; echo = 2;
3423  list inj = 1,5,7,41;
3424  list non_inj = 1,2,3,4;
3425  is_injective(inj);
3426  is_injective(non_inj);
3427}
3428
3429proc is_h_injective(list a, int h)
3430"USAGE:   is_h_injective(a, h)
3431RETURN:  1 if a is h-injective, 0 otherwise
3432EXAMPLE: is_h_injective; shows an example;
3433"
3434{
3435  //Create all sets of subsets
3436  list subsetlist = set_subset_set(a);
3437  list h_subsetlist;
3438  //delete every list with elements more than h+1 since they are not needed to check h-injectivity
3439  for (int i=1; i<=size(subsetlist); i++)
3440  {
3441    if(size(subsetlist[i])<=h)
3442    {
3443      h_subsetlist = insert(h_subsetlist,subsetlist[i]);
3444    }
3445  }
3446
3447  //Check if the remaining max_sums do not occure more than once
3448  list checklist=calculate_max_sum(h_subsetlist[1]);
3449  int calculator;
3450  for (i=2; i<=size(h_subsetlist);i++)
3451  {
3452    calculator = calculate_max_sum(h_subsetlist[i]);
3453    if (set_contains(checklist, calculator)==1)
3454    {
3455      return(0);
3456    }
3457    else
3458    {
3459      checklist = insert(checklist,calculator);
3460    }
3461  }
3462  return(1);
3463
3464}
3465example
3466{
3467  "EXAMPLE:"; echo = 2;
3468  list h_inj = 1,2,4,10,17;
3469  is_h_injective(h_inj,3);
3470  //1+2+4+10=17
3471  is_h_injective(h_inj,4);
3472}
3473
3474proc is_fix_injective(list a)
3475"USAGE:   is_fix_injective(a)
3476RETURN:  1 if a is fix-injective, 0 otherwise
3477EXAMPLE: is_fix_injective; shows an example;
3478"
3479{
3480  //Generation of the list-list-list
3481  list subsetlist = set_subset_set(a);
3483  list listoflists;
3484  list emptylist1;
3485  list worklist;
3486  int i;
3487  int v;
3488  list checklist;
3489  int calculator;
3490
3491  int set_destination;
3492
3493  //create list of lists wich contain the lists of a certain length as elements
3494  for (i = 1; i<= size(subsetlist); i++)
3495  {
3496    //Determine the size of the acutal list to choose where to insert it in the listoflists
3497    set_destination = size(subsetlist[i]);
3499    {
3500      //There is already an element with the same set size, so just insert it
3501      listoflists[set_destination] = insert(listoflists[set_destination],subsetlist[i]);
3502    }
3503    else
3504    {
3505      //There is not yet an element with the same set size, so create a new one
3506      listoflists[set_destination] = insert(emptylist1,subsetlist[i]);
3508    }
3509  }
3510
3511  //Check for injectivity of each seperate list. Works as in injectivity or h-injectivity
3512  for (v=1; v<=size(listoflists); v++)
3513  {
3514    worklist = listoflists[v];
3515
3516    checklist=calculate_max_sum(worklist[1]);
3517    for (i=2; i<=size(worklist); i++)
3518    {
3519      calculator = calculate_max_sum(worklist[i]);
3520      if (set_contains(checklist, calculator)==1)
3521      {
3522        return(0);
3523      }
3524      else
3525      {
3526        checklist = insert(checklist,calculator);
3527      }
3528    }
3529  }
3530  return(1);
3531}
3532example
3533{
3534  "EXAMPLE:"; echo = 2;
3535  //this is fix-injective because 17=10+2+4+1 with different numbers of addens.
3536  list fix_inj = 1,2,4,10,17;
3537  //this is not fix-injective because 4+1=2+3.
3538  list not_fix_inj = 1,2,3,4;
3539  is_fix_injective(fix_inj);
3540  is_fix_injective(not_fix_inj);
3541}
3542
3543proc three_elements(list out, int iterations)
3544"USAGE:   three_elements(out, iterations)
3545RETURN:  Injective_knapsack created with the three elements method
3546EXAMPLE: three_elements; shows an example;
3547"
3548{
3549  int a;
3550  int b;
3551  int c;
3552  int subsum;
3554  int condition = 0;
3555  out = set_turn(out);
3556  if (is_injective(out)==0)
3557  {
3558    return(0);
3559  }
3560  else
3561  {
3562    for (int i=1; i<=iterations; i++)
3563    {
3564      while(condition==0)
3565      {
3566        subsum = calculate_max_sum(out);
3569        c = b+subsum+1;
3570        if ((a+b)>(c+subsum))
3571        {
3572          condition=1;
3573        }
3574        else
3575        {
3577        }
3578      }
3580      condition=0;
3581      out=set_insert(out, a);
3582      out=set_insert(out, b);
3583      out=set_insert(out, c);
3584    }
3585    return(out);
3586  }
3587}
3588example
3589{
3590  "EXAMPLE:"; echo = 2;
3591  //this is fix-injective because 17=10+2+4+1 with different numbers of addens.
3592  list super_increasing = 1,2,4,10,20;
3593  list a = three_elements(super_increasing,2);
3594  a;
3595  is_injective(a);
3596}
3597/*
3598//===============================================================
3599//=======  Example for DSA  =====================================
3600//===============================================================
3601Suppose a file test is given.It contains "Oscar".
3602
3603//Hash-function MD5 under Linux
3604
3605md5sum test           8edfe37dae96cfd2466d77d3884d4196
3606
3607//================================================================
3608
3609ring R=0,x,dp;
3610
3611number q=2^19+21;          //524309
3612number o=2*3*23*number(7883)*number(16170811);
3613
3614number p=o*q+1;            //9223372036869000547
3615number b=2;
3616number g=power(2,o,p);     //8308467587808723131
3617
3618number a=111111;
3619number A=power(g,a,p);     //8566038811843553785
3620
3621number h =decimal("8edfe37dae96cfd2466d77d3884d4196");
3622
3623                           //189912871665444375716340628395668619670
3624h= h mod q;                //259847
3625
3626number k=123456;
3627
3628number ki=exgcd(k,q)[1];    //50804
3629                            //inverse von k mod q
3630
3631number r= power(g,k,p) mod q;    //76646
3632
3633number s=ki*(h+a*r) mod q;       //2065
3634
3635//========== signatur is (r,s)=(76646,2065) =====================
3636//==================== verification  ============================
3637
3638number si=exgcd(s,q)[1];  //inverse von s mod q
3639number e1=si*h mod q;
3640number e2=si*r mod q;
3641number rr=((power(g,e1,p)*power(A,e2,p)) mod p) mod q;    //76646
3642
3643//===============================================================
3644//=======  Example for knapsack  ================================
3645//===============================================================
3646ring R=(5^5,t),x,dp;
3647R;
3648//   # ground field : 3125
3649//   primitive element : t
3650//   minpoly        : 1*t^5+4*t^1+2*t^0
3651//   number of vars : 1
3652//        block   1 : ordering dp
3653//                  : names    x
3654//        block   2 : ordering C
3655
3656proc findEx(number n, number g)
3657{
3658   int i;
3659   for(i=0;i<=size(basering)-1;i++)
3660   {
3661      if(g^i==n){return(i);}
3662   }
3663}
3664
3665number g=t^3;                       //choice of the primitive root
3666
3667findEx(t+1,g);
3668//2091
3669findEx(t+2,g);
3670//2291
3671findEx(t+3,g);
3672//1043
3673
3674intvec b=1,2091,2291,1043;          // k=4
3675int z=199;
3676intvec v=1043+z,1+z,2091+z,2291+z;  //permutation   pi=(0123)
3677v;
36781242,200,2290,2490
3679
3680//(1101)=(e_3,e_2,e_1,e_0)
3681//encoding 2490+2290+1242=6022 und 1+1+0+1=3
3682
3683//(6022,3) decoding: c-z*c'=6022-199*3=5425
3684
3685ring S=5,x,dp;
3686poly F=x5+4x+2;
3687poly G=reduce((x^3)^5425,std(F));
3688G;
3689//x3+x2+x+1
3690
3691factorize(G);
3692//[1]:
3693//   _[1]=1
3694//   _[2]=x+1
3695//   _[3]=x-2
3696//   _[4]=x+2
3697//[2]:
3698//   1,1,1,1
3699
3700//factors x+1,x+2,x+3, i.e.  (1110)=(e_pi(3),e_pi(2),e_pi(1),e_pi(0))
3701
3702//pi(0)=1,pi(1)=2,pi(2)=3,pi(3)=0 gives:  (1101)
3703
3704*/
3705
Note: See TracBrowser for help on using the repository browser.