//////////////////////////////////////////////////////////////////////////////
version="version crypto.lib 4.0.0.1 Aug_2014 "; // $Id$
category="Teaching";
info="
LIBRARY:  crypto.lib     Procedures for teaching cryptography
AUTHORS:                 Gerhard Pfister, pfister@mathematik.uni-kl.de
@*                       David Brittinger, dativ@gmx.net

OVERVIEW:
     The library contains procedures to compute the discrete logarithm,
      primality-tests, factorization included elliptic curves.
      The library is intended to be used for teaching purposes but not
      for serious computations. Sufficiently high printlevel allows to
      control each step, thus illustrating the algorithms at work.


PROCEDURES:
 decimal(s);                number corresponding to the hexadecimal number s
 eexgcdN(L)                 T with sum_i L[i]*T[i]=T[n+1]=gcd(L[1],...,L[n])
 lcmN(a,b)                  compute lcm(a,b)
 powerN(m,d,n)              compute m^d mod n
 chineseRem(T,L)            compute x such that x = T[i] mod L[i]
 Jacobi(a,n)                the generalized Legendre symbol of a and n
 primList(n)                the list of all primes <=n
 primL(q)                   first primes p_1,...,p_r such that q<p_1*...*p_r
 intPart(x)                 the integral part of a rational number
 intRoot(m)                 the integral part of the square root of m
 squareRoot(a,p)            the square root of a in Z/p, p prime
 solutionsMod2(M)           basis solutions of Mx=0 over Z/2
 powerX(q,i,I)              q-th power of the i-th variable modulo I
 babyGiant(b,y,p)           discrete logarithm x: b^x=y mod p
 rho(b,y,p)                 discrete logarithm x: b^x=y mod p
 MillerRabin(n,k)           probabilistic primaly-test of Miller-Rabin
 SolowayStrassen(n,k)       probabilistic primaly-test of Soloway-Strassen
 PocklingtonLehmer(N,[])    primaly-test of Pocklington-Lehmer
 PollardRho(n,k,a,[])       Pollard's rho factorization
 pFactor(n,B,P)             Pollard's p-factorization
 quadraticSieve(n,c,B,k)    quadratic sieve factorization
 isOnCurve(N,a,b,P)         P is on the curve y^2z=x^3+a*xz^2+b*z^3  over Z/N
 ellipticAdd(N,a,b,P,Q)     P+Q, addition on elliptic curves
 ellipticMult(N,a,b,P,k)    k*P on elliptic curves
 ellipticRandomCurve(N)     generates y^2z=x^3+a*xz^2+b*z^3  over Z/N randomly
 ellipticRandomPoint(N,a,b) random point on y^2z=x^3+a*xz^2+b*z^3  over Z/N
 countPoints(N,a,b)         number of points of y^2=x^3+a*x+b  over Z/N
 ellipticAllPoints(N,a,b)   points of y^2=x^3+a*x+b  over Z/N
 ShanksMestre(q,a,b,[])     number of points of y^2=x^3+a*x+b  over Z/N
 Schoof(N,a,b)              number of points of y^2=x^3+a*x+b  over Z/N
 generateG(a,b,m)           m-th division polynomial of y^2=x^3+a*x+b  over Z/N
 factorLenstraECM(N,S,B,[]) Lenstra's factorization
 ECPP(N)                    primaly-test of Goldwasser-Kilian
 calculate_ordering(num1, primitive, mod1)    Calculates x so that primitive^x == num1 mod mod1
  is_primitive_root(primitive, mod1) Checks if primitive is a primitive root modulo mod1
  find_first_primitive_root(mod1) Returns the first primitive root modulo mod1, starting with 1
  binary_add(binary_list)         Adds a 1 to a binary encoded list
  inverse_modulus(num,mod1)       Finds a t so that t*num = 1 mod mod1
  is_prime(n)                     Checks if n is prime
  proc find_biggest_index(a)      Returns the index of the biggest element of a
  find_index(a,e)                 Returns the list index of element e in list a. Returns 0 if e is not in a
  subset_sum01(list knapsack, int solution)                          solves the subset-sum-knapsack-problem by calculating all subsets and choosing the right solution
  subset_sum02(list knapsack, int sol)                            solves the subset-sum-knapsack-problem with a naive greedy algorithm
  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
  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
  naccache_stern_generation(int key, int primenum)                      generates a hard knapsack for the Naccache-Stern Kryptosystem for given key and prime modulus
  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
  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
  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
  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
  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
  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
  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
  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
  super_increasing_knapsack(int ksize)                            Creates the smallest super-increasing knapsack of given size ksize
  h_increasing_knapsack(int ksize, int h)                            Creates the smallest h-increasing knapsack of given size ksize and h
  injective_knapsack(int ksize, int kmaxelement)                        Creates all list of all injective knapsacks of given size ksize and maximal element kmaxelement
  calculate_max_sum(list a)                                  Calculates the maximal sum of a given knapsack a
  set_is_injective(list a)                                    Checks if knapsack a is injective
  is_h_injective(list a, int h)                                Checks if knapsack a is h-injective
  is_fix_injective(list a)                                  Checks if knapsack a is fix-injective
  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

              [parameters in square brackets are optional]
";

LIB "poly.lib";
LIB "atkins.lib";

///////////////////////////////////////////////////////////////////////////////


//=============================================================================
//=========================== basic prozedures ================================
//=============================================================================

proc decimal(string s)
"USAGE:  decimal(s); s = string
RETURN: the (decimal) number corresponding to the hexadecimal number s
EXAMPLE:example decimal; shows an example
"
{
   int n=size(s);
   int i;
   bigint k;
   bigint t=16;
   bigint m=0;
   for(i=1;i<=n;i++)
   {
      k=0;
      if(s[i]=="1"){k=1;}
      if(s[i]=="2"){k=2;}
      if(s[i]=="3"){k=3;}
      if(s[i]=="4"){k=4;}
      if(s[i]=="5"){k=5;}
      if(s[i]=="6"){k=6;}
      if(s[i]=="7"){k=7;}
      if(s[i]=="8"){k=8;}
      if(s[i]=="9"){k=9;}
      if(s[i]=="a"){k=10;}
      if(s[i]=="b"){k=11;}
      if(s[i]=="c"){k=12;}
      if(s[i]=="d"){k=13;}
      if(s[i]=="e"){k=14;}
      if(s[i]=="f"){k=15;}
      m=m*t+k;
   }
   return(m);
}
example
{ "EXAMPLE:"; echo = 2;
   string s  ="8edfe37dae96cfd2466d77d3884d4196";
   decimal(s);
}

proc eexgcdN(list L)
"USAGE:  eexgcdN(L);
RETURN: list T such that sum_i L[i]*T[i]=T[n+1]=gcd(L[1],...,L[n])
EXAMPLE:example eexgcdN; shows an example
"
{
   if(size(L)==2)
   {
     list LL=extgcd(L[1],L[2]);return(list(LL[2],LL[3],LL[1]));
   }
   bigint p=L[size(L)];
   L=delete(L,size(L));
   list T=eexgcdN(L);
   list S=extgcd(T[size(T)],p);
   int i;
   for(i=1;i<=size(T)-1;i++)
   {
      T[i]=T[i]*S[2];
   }
   p=T[size(T)];
   T[size(T)]=S[3];
   T[size(T)+1]=S[1];
   return(T);
}
example
{ "EXAMPLE:"; echo = 2;
   eexgcdN(list(24,15,21));
}

proc lcmN(bigint a, bigint b)
"USAGE:  lcmN(a,b);
RETURN: lcm(a,b);
EXAMPLE:example lcmN; shows an example
"
{
   return (a*b/gcd(a,b));
}
example
{ "EXAMPLE:"; echo = 2;
   lcmN(24,15);
}

proc powerN(bigint m, bigint d, bigint n)
"USAGE:  powerN(m,d,n);
RETURN: m^d mod n
EXAMPLE:example powerN; shows an example
"
{
   if(d==0){return(bigint(1));}
   int i;
   if(n==0)
   {
      for(i=12;i>=2;i--)
      {
         if((d mod i)==0){return(powerN(m,d div i,n)^i);}
      }
      return(m*powerN(m,d-1,n));
   }
   for(i=12;i>=2;i--)
   {
      if((d mod i)==0)
      {
        bigint rr=powerN(m,d div i,n)^i mod n;
        if (rr<0) { rr=rr+n;}
        return(rr);
      }
   }
   return(m*powerN(m,d-1,n) mod n);
}
example
{ "EXAMPLE:"; echo = 2;
   powerN(24,15,7);
}

proc chineseRem(list T,list L)
"USAGE:  chineseRem(T,L);
RETURN: x such that x = T[i] mod L[i]
NOTE:   chinese remainder theorem
EXAMPLE:example chineseRem; shows an example
"
{
   int i;
   int n=size(L);
   bigint N=1;
   for(i=1;i<=n;i++)
   {
      N=N*L[i];
   }
   list M;
   for(i=1;i<=n;i++)
   {
      M[i]=N div L[i];
   }
   list S=eexgcdN(M);
   bigint x;
   for(i=1;i<=n;i++)
   {
      x=x+S[i]*M[i]*T[i];
   }
   x=x mod N;
   return(x);
}
example
{ "EXAMPLE:"; echo = 2;
   chineseRem(list(24,15,7),list(2,3,5));
}

proc Jacobi(bigint a, bigint n)
"USAGE:  Jacobi(a,n);
RETURN: the generalized Legendre symbol
NOTE: if n is an odd prime then Jacobi(a,n)=0,1,-1 if n|a, a=x^2 mod n,else
EXAMPLE:example Jacobi; shows an example
"
{
   int i;
   int z=1;
   bigint t=1;
   bigint k;

   if((((n-1) div 2) mod 2)!=0){z=-1;}
   if(a<0){return(z*Jacobi(-a,n));}
   a=a mod n;
   if(n==1){return(1);}
   if(a==0){return(0);}

   while(a!=0)
   {
      while((a mod 2)==0)
      {
         a=a div 2;
         if(((n mod 8)==3)||((n mod 8)==5)){t=-t;}
      }
      k=a;a=n;n=k;
      if(((a mod 4)==3)&&((n mod 4)==3)){t=-t;}
      a=a mod n;
   }
   if (n==1){return(t);}
   return(0);
}
example
{ "EXAMPLE:"; echo = 2;
   Jacobi(13580555397810650806,5792543);
}

proc primList(int n)
"USAGE:  primList(n);
RETURN: the list of all primes <=n
EXAMPLE:example primList; shows an example
"
{
   int i,j;
   list re;
   re[1]=2;
   re[2]=3;
   for(i=5;i<=n;i=i+2)
   {
     j=1;
     while(j<=size(re))
     {
        if((i mod re[j])==0){break;}
        j++;
     }
     if(j==size(re)+1){re[size(re)+1]=i;}
   }
   return(re);
}
example
{ "EXAMPLE:"; echo = 2;
   list L=primList(100);
   size(L);
   L[size(L)];
}

proc primL(bigint q)
"USAGE:  primL(q);
RETURN: list of the first primes p_1,...,p_r such that q>p_1*...*p_(r-1)
        and q<p_1*...*p_r
EXAMPLE:example primL; shows an example
"
{
   int i,j;
   list re;
   re[1]=2;
   re[2]=3;
   bigint s=6;
   i=3;
   while(s<=q)
   {
     i=i+2;
     j=1;
     while(j<=size(re))
     {
        if((i mod re[j])==0){break;}
        j++;
     }
     if(j==size(re)+1)
     {
        re[size(re)+1]=i;
        s=s*i;
     }
   }
   return(re);
}
example
{ "EXAMPLE:"; echo = 2;
   primL(20);
}

proc intPart(number x)
"USAGE:  intPart(x);
RETURN: the integral part of a rational number
EXAMPLE:example intPart; shows an example
"
{
  if (x>=0)
  {
    return(bigint((numerator(x)-(bigint(numerator(x)) mod bigint(denominator(x)))))
         div bigint(denominator(x)));
  }
  else
  {
    return(bigint((numerator(x)-(bigint(numerator(x)) mod bigint(denominator(x)+denominator(x)))))
         div bigint(denominator(x)));
  }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,x,dp;
   intPart(7/3);
}

proc intRoot(bigint mm)
"USAGE:  intRoot(m);
RETURN: the integral part of the square root of m
EXAMPLE:example intRoot; shows an example
"
{
   ring R = 0,@x,dp;
   number m=mm;
   number x=1;
   number t=x^2;
   number s=(x+1)^2;
   while(((m>t)&&(m>s))||((m<t)&&(m<s)))
   {
      if (x==0) { t=0; }
      else
      {
        x=intPart(x/2+m/(2*x));  //Newton step
        t=x^2;
      }
      if(t>m)
      {
         s=(x-1)^2;
      }
      else
      {
         s=(x+1)^2;
      }
   }
   if(t>m){return(bigint(x-1));}
   if(s==m){return(bigint(x+1));}
   return(bigint(x));
}
example
{ "EXAMPLE:"; echo = 2;
   intRoot(20);
}

proc squareRoot(bigint a, bigint p)
"USAGE:  squareRoot(a,p);
RETURN: the square root of a in Z/p, p prime
NOTE:   assumes the Jacobi symbol is 1 or p=2.
EXAMPLE:example squareRoot; shows an example
"
{
   if(p==2){return(a);}
   if((a mod p)==0){return(0);}
   if (a<0) { a=a+p; }
   if(powerN(a,p-1,p)!=1)
   {
      "p is not prime";
      return(bigint(-5));
   }
   bigint n=random(1,2147483647) mod p;
   if(n==0){n=n+1;}
   bigint j=Jacobi(n,p);
   if(j==0)
   {
      "p is not prime";
      return(bigint(-5));
   }
   if(j==1)
   {
      return(squareRoot(a,p));
   }
   bigint q=p-1;
   bigint e=0;
   bigint two=2;
   bigint z,m,t;
   while((q mod 2)==0)
   {
      e=e+1;
      q=q div 2;
   }
   bigint y=powerN(n,q,p);
   bigint r=e;
   bigint x=powerN(a,(q-1) div 2,p);
   bigint b=a*x^2 mod p;
   x=a*x mod p;

   while(((b-1) mod p)!=0)
   {
      m=0;z=b;
      while(((z-1) mod p)!=0)
      {
         z=z^2 mod p;
         m=m+1;
      }
      t=powerN(y,powerN(two,r-m-1,p),p);
      y=t^2 mod p;
      r=m;
      x=x*t mod p;
      b=b*y mod p;
   }
   return(x);
}
example
{ "EXAMPLE:"; echo = 2;
   squareRoot(8315890421938608,32003);
}


proc solutionsMod2(bigintmat MM)
"USAGE:  solutionsMod2(M);
RETURN: an intmat containing a basis of the vector space of solutions of the
        linear system of equations defined by M over the prime field of
        characteristic 2
EXAMPLE:example solutionsMod2; shows an example
"
{
   ring Rhelp=2,z,(c,dp);
   int i,j;
   matrix M[nrows(MM)][ncols(MM)];
   for(i=1;i<=nrows(MM);i++)
   {
      for(j=1;j<=ncols(MM);j++)
      {
         M[i,j]=MM[i,j];
      }
   }
   matrix S=syz(M);
   intmat v[nrows(S)][ncols(S)];
   for(i=1;i<=nrows(S);i++)
   {
      for(j=1;j<=ncols(S);j++)
      {
         if(S[i,j]==1){v[i,j]=1;}
      }
   }
   return(v);
}
example
{ "EXAMPLE:"; echo = 2;
   bigintmat M[3][3]=1,2,3,4,5,6,7,6,5;
   solutionsMod2(M);
}

proc powerX(int q, int i, ideal I)
"USAGE:  powerX(q,i,I);
RETURN: the q-th power of the i-th variable modulo I
ASSUME: I is a standard basis
EXAMPLE:example powerX; shows an example
"
{
   if(q<=181){return(reduce(var(i)^int(q),I));}
   if((q mod 5)==0){return(reduce(powerX(q div 5,i,I)^5,I));}
   if((q mod 4)==0){return(reduce(powerX(q div 4,i,I)^4,I));}
   if((q mod 3)==0){return(reduce(powerX(q div 3,i,I)^3,I));}
   if((q mod 2)==0){return(reduce(powerX(q div 2,i,I)^2,I));}
   return(reduce(var(i)*powerX(q-1,i,I),I));
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,(x,y),dp;
   powerX(100,2,std(ideal(x3-1,y2-x)));
}

//======================================================================
//=========================== Discrete Logarithm =======================
//======================================================================

//============== Shank's baby step - giant step ========================

proc babyGiant(bigint b, bigint y, bigint p)
"USAGE:  babyGiant(b,y,p);
RETURN: the discrete logarithm x: b^x=y mod p
NOTE:   This procedure works based on Shank's baby step - giant step method.
EXAMPLE:example babyGiant; shows an example
"
{
   int i,j,m;
   list l;
   bigint h=1;
   bigint x;

//choose m minimal such that m^2>p
   for(i=1;i<=p;i++){if(i^2>p) break;}
   m=i;

//baby-step:  compute the table y*b^i for 1<=i<=m
   for(i=1;i<=m;i++){l[i]=y*b^i mod p;}

//giant-step: compute b^(m+j), 1<=j<=m and search in the baby-step table
//for an i with y*b^i=b^(m*j). If found then x=m*j-i
   bigint g=b^m mod p;
   while(j<m)
   {
      j++;
      h=h*g mod p;
      for(i=1;i<=m;i++)
      {
         if(h==l[i])
         {
            x=m*j-i;
            j=m;
            break;
         }
     }
   }
   return(x);
}
example
{ "EXAMPLE:"; echo = 2;
   bigint b=2;
   bigint y=10;
   bigint p=101;
   babyGiant(b,y,p);
}

//==============  Pollards rho  =================================

proc rho(bigint b, bigint y, bigint p)
"USAGE:  rho(b,y,p);
RETURN: the discrete logarithm x=log_b(y): b^x=y mod p
NOTE: Pollard's rho:
       choose random f_0 in 0,...,p-2 ,e_0=0, define x_0=b^f_0, define
       x_i=y^e_i*b^f_i as below. For i large enough there is i with
       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,
       d=gcd(s,p-1)=u*s+v*(p-1) then x=tu/d +j*(p-1)/d for some j (to be
       found by trying)
EXAMPLE:example rho; shows an example
"
{
   int i=1;
   int j;
   bigint s,t;
   list e,f,x;

   e[1]=0;
   f[1]=random(0,2147483629) mod (p-1);
   x[1]=powerN(b,f[1],p);
   while(i)
   {
      if((x[i] mod 3)==1)
      {
         x[i+1]=y*x[i] mod p;
         e[i+1]=e[i]+1 mod (p-1);
         f[i+1]=f[i];
      }
      if((x[i] mod 3)==2)
      {
         x[i+1]=x[i]^2 mod p;
         e[i+1]=e[i]*2 mod (p-1);
         f[i+1]=f[i]*2 mod (p-1);
      }
      if((x[i] mod 3)==0)
      {
         x[i+1]=x[i]*b mod p;
         e[i+1]=e[i];
         f[i+1]=f[i]+1 mod (p-1);
      }
      i++;
      for(j=i-1;j>=1;j--)
      {
         if(x[i]==x[j])
         {
            s=(e[j]-e[i]) mod (p-1);
            t=(f[i]-f[j]) mod (p-1);
            if(s!=0)
            {
               i=0;
            }
            else
            {
               e[1]=0;
               f[1]=random(0,2147483629) mod (p-1);
               x[1]=powerN(b,f[1],p);
               i=1;
            }
            break;
         }
      }
   }

   list w=extgcd(s,p-1);
   bigint u=w[2];
   bigint d=w[1];

   bigint a=(t*u div d) mod (p-1);

   bigint pn=powerN(b,a,p);
   if (pn<0) { pn=pn+p;}
   while(powerN(b,a,p)!=y)
   {
      a=(a+(p-1) div d) mod (p-1);
      if (a<0) { a=a+p-1; }
   }
   return(a);
}
example
{ "EXAMPLE:"; echo = 2;
   bigint b=2;
   bigint y=10;
   bigint p=101;
   rho(b,y,p);
}
//====================================================================
//====================== Primality Tests =============================
//====================================================================

//================================= Miller-Rabin =====================

proc MillerRabin(bigint n, int k)
"USAGE:  MillerRabin(n,k);
RETURN: 1 if n is prime, 0 else
NOTE: probabilistic test of Miller-Rabin with k loops to test if n is prime.
       Using the theorem: If n is prime, n-1=2^s*r, r odd, then
       powerN(a,r,n)=1 or powerN(a,r*2^i,n)=-1 for some i
EXAMPLE:example MillerRabin; shows an example
"
{
   if(n<0){n=-n;}
   if((n==2)||(n==3)){return(1);}
   if((n mod 2)==0){return(0);}

   int i;
   bigint a,b,j,r,s;
   r=n-1;
   s=0;
   while((r mod 2)==0)
   {
      s=s+1;
      r=r div 2;
   }
   while(i<k)
   {
      i++;
      a=random(2,2147483629) mod n; if(a==0){a=3;}
      if(gcd(a,n)!=1){return(0);}
      b=powerN(a,r,n);
      if(b!=1)
      {
         j=0;
         while(j<s)
         {
           if(((b+1) mod n)==0) break;
           b=powerN(b,2,n);
           j=j+1;
         }
         if(j==s){return(0);}
      }
   }
   return(1);
}
example
{ "EXAMPLE:"; echo = 2;
   bigint x=2;
   x=x^787-1;
   MillerRabin(x,3);
}

//======================= Soloway-Strassen  ==========================

proc SolowayStrassen(bigint n, int k)
"USAGE:  SolowayStrassen(n,k);
RETURN: 1 if n is prime, 0 else
NOTE: probabilistic test of Soloway-Strassen with k loops to test if n is
       prime using the theorem: If n is prime then
       powerN(a,(n-1)/2,n)=Jacobi(a,n) mod n
EXAMPLE:example SolowayStrassen; shows an example
"
{
   if(n<0){n=-n;}
   if((n==2)||(n==3)){return(1);}
   if((n mod 2)==0){return(0);}

   bigint a,pn,jn;
   int i;
   while(i<k)
   {
      i++;
      a=random(2,2147483629) mod n; if(a==0){a=3;}
      if(gcd(a,n)!=1){return(0);}
      pn=powerN(a,(n-1) div 2,n);
      if (pn<0) { pn=pn+n;}
      jn=Jacobi(a,n) mod n;
      if (jn<0) { jn=jn+n;}
      if(pn!=jn){return(0);}
   }
   return(1);
}
example
{ "EXAMPLE:"; echo = 2;
   bigint h=10;
   bigint p=h^100+267;
   //p=h^100+43723;
   //p=h^200+632347;
   SolowayStrassen(h,3);
}


/*
ring R=0,z,dp;
number p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
number q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
number n=p*q;
SolowayStrassen(n,3);
*/

//===================== Pocklington-Lehmer ==============================

proc PocklingtonLehmer(bigint N, list #)
"USAGE: PocklingtonLehmer(N); optional: PocklingtonLehmer(N,L);
        L a list of the first k primes
RETURN:message N is not prime or {A,{p},{a_p}} as certificate for N being prime
NOTE:assumes that it is possible to factorize N-1=A*B such that gcd(A,B)=1,
      the factorization of A is completely known and A^2>N .
      N is prime if and only if for each prime factor p of A we can find
      a_p such that a_p^(N-1)=1 mod N and gcd(a_p^((N-1)/p)-1,N)=1
EXAMPLE:example PocklingtonLehmer; shows an example
"
{
   bigint m=intRoot(N);
   if(size(#)>0)
   {
      list S=PollardRho(N-1,10000,1,#);
   }
   else
   {
      list S=PollardRho(N-1,10000,1);
   }
   int i,j;
   bigint A=1;
   bigint p,a,g;
   list PA;
   list re;

   while(i<size(S))
   {
      p=S[i+1];
      A=A*p;
      PA[i+1]=p;
      if(A>m){break;}

      while(1)
      {
        p=p*S[i+1];
        if(((N-1) mod p)==0)
        {
           A=A*p;
        }
        else
        {
           break;
        }
      }
      i++;
   }
   if(A<=m)
   {
     A=N div A;
     PA=list(S[size(S)]);
   }
   for(i=1;i<=size(PA);i++)
   {
      a=1;
      while(a<N-1)
      {
         a=a+1;
         if(powerN(a,N-1,N)!=1){return("not prime");}
         g=gcd(powerN(a,(N-1) div PA[i],N),N);
         if(g==1)
         {
           re[size(re)+1]=list(PA[i],a);
           break;
         }
         if(g<N){"not prime";return(g);}
      }
   }
   return(list(A,re));
}
example
{ "EXAMPLE:"; echo = 2;
   bigint N=105554676553297;
   PocklingtonLehmer(N);
   list L=primList(1000);
   PocklingtonLehmer(N,L);
}

//=======================================================================
//======================= Factorization =================================
//=======================================================================

//======================= Pollards rho  =================================

proc PollardRho(bigint n, int k, int allFactors, list #)
"USAGE:  PollardRho(n,k,allFactors); optional: PollardRho(n,k,allFactors,L);
         L a list of the first k primes
RETURN: a list of factors of n (which could be just n),if allFactors=0@*
        a list of all factors of n ,if allFactors=1
NOTE: probabilistic rho-algorithm of Pollard to find a factor of n in k loops.
      Creates a sequence x_i such that (x_i)^2=(x_2i)^2 mod n for some i,
      computes gcd(x_i-x_2i,n) to find a divisor. To define the sequence
      choose x,a and define x_n+1=x_n^2+a mod n, x_1=x.
      If allFactors is 1, it tries to find recursively all prime factors
      using the Soloway-Strassen test.
SEE ALSO: primefactors
EXAMPLE:example PollardRho; shows an example
"
{
    int i,j;
    list L=primList(100);
    list re,se;
    if(n<0){n=-n;}
    if(n==1){return(re);}

//this is optional: test whether a prime of the list # devides n
    if(size(#)>0)
    {
       L=#;
    }
    for(i=1;i<=size(L);i++)
    {
       if((n mod L[i])==0)
       {
          re[size(re)+1]=L[i];
          while((n mod L[i])==0)
          {
             n=n div L[i];
          }
       }
       if(n==1){return(re);}
    }
    int e=size(re);
//here the rho-algorithm starts
    bigint a,d,x,y;
    while(n>1)
    {
       a=random(2,2147483629);
       x=random(2,2147483629);
       y=x;
       d=1;
       i=0;
       while(i<k)
       {
          i++;
          x=powerN(x,2,n); x=(x+a) mod n;
          y=powerN(y,2,n); y=(y+a) mod n;
          y=powerN(y,2,n); y=(y+a) mod n;
          d=gcd(x-y,n);
          if(d>1)
          {
             re[size(re)+1]=d;
             while((n mod d)==0)
             {
               n=n div d;
             }
             break;
          }
          if(i==k)
          {
             re[size(re)+1]=n;
             n=1;
          }
       }
    }
    if(allFactors)      //want to obtain all prime factors
    {
       i=e;
       while(i<size(re))
       {
          i++;

          if(!SolowayStrassen(re[i],5))
          {
             se=PollardRho(re[i],2*k,1);
             re[i]=se[size(se)];
             for(j=1;j<=size(se)-1;j++)
             {
                re[size(re)+1]=se[j];
             }
             i--;
          }
       }
    }
    return(re);
}
example
{ "EXAMPLE:"; echo = 2;
   bigint h=10;
   bigint p=h^30+4;
   PollardRho(p,5000,0);
}

//======================== Pollards p-factorization ================
proc pFactor(bigint n,int B, list P)
"USAGE:  pFactor(n,B,P); n to be factorized, B a bound , P a list of primes
RETURN: a list of factors of n or n if no factor found
NOTE: Pollard's p-factorization
       creates the product k of powers of primes (bounded by B)  from
       the list P with the idea that for a prime divisor p of n we have
       p-1|k, and then p divides gcd(a^k-1,n) for some random a
EXAMPLE:example pFactor; shows an example
"
{
   int i;
   bigint k=1;
   bigint w;
   while(i<size(P))
   {
      i++;
      w=P[i];
      if(w>B) break;
      while(w*P[i]<=B)
      {
         w=w*P[i];
      }
      k=k*w;
   }
   bigint a=random(2,2147483629);
   bigint d=gcd(powerN(a,k,n)-1,n);
   if((d>1)&&(d<n)){return(d);}
   return(n);
}
example
{ "EXAMPLE:"; echo = 2;
   list L=primList(1000);
   pFactor(1241143,13,L);
   bigint h=10;
   h=h^30+25;
   pFactor(h,20,L);
}

//==================== quadratic sieve ==============================

proc quadraticSieve(bigint n, int c, list B, int k)
"USAGE:  quadraticSieve(n,c,B,k); n to be factorized, [-c,c] the
         sieve-intervall, B a list of primes,
         k for using the first k elements in B
RETURN: a list of factors of n or the message: no divisor found
NOTE: The idea being used is to find x,y such that x^2=y^2 mod n then
      gcd(x-y,n) can be a proper divisor of n
EXAMPLE:example quadraticSieve; shows an example
"
{
   bigint f,d;
   int i,j,l,s,p;
   list S,tmp;
   intvec v;
   v[k]=0;

//compute the integral part of the square root of n
   bigint m=intRoot(n);

//consider the function f(X)=(X+m)^2-n and compute for s in [-c,c] the values
   while(i<=2*c)
   {
      f=(i-c+m)^2-n;
      tmp[1]=i-c+m;
      tmp[2]=f;
      tmp[3]=v;
      S[i+1]=tmp;
      i++;
   }

//the sieve with p in B
//find all s in [-c,c] such that f(s) has all prime divisors in the first
//k elements of B and the decomposition of f(s). They are characterized
//by 1 or -1 at the second place of S[j]:
//S[j]=j-c+m,f(j-c)/p_1^v_1*...*p_k^v_k, v_1,...,v_k maximal
   for(i=1;i<=k;i++)
   {
      p=B[i];
      if((p>2)&&(Jacobi(n,p)==-1)){i++;continue;}//n is no quadratic rest mod p
      j=1;
      while(j<=p)
      {
         if(j>2*c+1) break;
         f=S[j][2];
         v=S[j][3];
         s=0;
         while((f mod p)==0)
         {
            s++;
            f=f div p;
         }
         if(s)
         {
           S[j][2]=f;
           v[i]=s;
           S[j][3]=v;
           l=j;
           while(l+p<=2*c+1)
           {
              l=l+p;
              f=S[l][2];
              v=S[l][3];
              s=0;
              while((f mod p)==0)
              {
                s++;
                f=f div p;
              }
              S[l][2]=f;
              v[i]=s;
              S[l][3]=v;
           }
         }
         j++;
      }
   }
   list T;
   for(j=1;j<=2*c+1;j++)
   {
      if((S[j][2]==1)||(S[j][2]==-1))
      {
         T[size(T)+1]=S[j];
      }
   }

//the system of equations for the exponents {l_s} for the f(s) such
//product f(s)^l_s is a square (l_s are 1 or 0)
   bigintmat M[k+1][size(T)];
   for(j=1;j<=size(T);j++)
   {
      if(T[j][2]==-1){M[1,j]=1;}
      for(i=1;i<=k;i++)
      {
         M[i+1,j]=T[j][3][i];
      }
   }
   intmat G=solutionsMod2(M);

//construction of x and y such that x^2=y^2 mod n and d=gcd(x-y,n)
//y=square root of product f(s)^l_s
//x=product s+m
   bigint x=1;
   bigint y=1;

   for(i=1;i<=ncols(G);i++)
   {
      kill v;
      intvec v;
      v[k]=0;
      for(j=1;j<=size(T);j++)
      {
         x=x*T[j][1]^G[j,i] mod n;
         if((T[j][2]==-1)&&(G[j,i]==1)){y=-y;}
         v=v+G[j,i]*T[j][3];

      }
      for(l=1;l<=k;l++)
      {
         y=y*B[l]^(v[l] div 2) mod n;
      }
      d=gcd(x-y,n);
      if((d>1)&&(d<n)){return(d);}
   }
   return("no divisor found");
}
example
{ "EXAMPLE:"; echo = 2;
   list L=primList(5000);
   quadraticSieve(7429,3,L,4);
   quadraticSieve(1241143,100,L,50);
}

//======================================================================
//==================== elliptic curves  ================================
//======================================================================

//================= elementary operations ==============================

proc isOnCurve(bigint N, bigint a, bigint b, list P)
"USAGE:  isOnCurve(N,a,b,P);
RETURN: 1 or 0 (depending on whether P is on the curve or not)
NOTE: checks whether P=(P[1]:P[2]:P[3]) is a point on the elliptic
      curve defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
EXAMPLE:example isOnCurve; shows an example
"
{
   if(((P[2]^2*P[3]-P[1]^3-a*P[1]*P[3]^2-b*P[3]^3) mod N)!=0){return(0);}
   return(1);
}
example
{ "EXAMPLE:"; echo = 2;
   isOnCurve(32003,5,7,list(10,16,1));
}

proc ellipticAdd(bigint N, bigint a, bigint b, list P, list Q)
"USAGE:  ellipticAdd(N,a,b,P,Q);
RETURN: list L, representing the point P+Q
NOTE: P=(P[1]:P[2]:P[3]), Q=(Q[1]:Q[2]:Q[3]) points on the elliptic curve
      defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
EXAMPLE:example ellipticAdd; shows an example
"
{
   if(N==2){ERROR("not implemented for 2");}
   int i;
   for(i=1;i<=3;i++)
   {
      P[i]=P[i] mod N; if (P[i]<0) { P[i]=P[i]+N:}
      Q[i]=Q[i] mod N; if (Q[i]<0) { Q[i]=Q[i]+N;}
   }
   list Resu;
   Resu[1]=bigint(0);
   Resu[2]=bigint(1);
   Resu[3]=bigint(0);
   list Error;
   Error[1]=0;
   //test for ellictic curve
   bigint D=4*a^3+27*b^2;
   bigint g=gcd(D,N);
   if(g==N){return(Error);}
   if(g!=1)
   {
      P[4]=g;
      return(P);
   }
   if(((P[1]==0)&&(P[2]==0)&&(P[3]==0))||((Q[1]==0)&&(Q[2]==0)&&(Q[3]==0)))
   {
      Error[1]=-2;
      return(Error);
   }
   if(!isOnCurve(N,a,b,P)||!isOnCurve(N,a,b,Q))
   {
      Error[1]=-1;
      return(Error);
   }
   if(P[3]==0){return(Q);}
   if(Q[3]==0){return(P);}
   list I=extgcd(P[3],N);
   if(I[1]!=1)
   {
      P[4]=I[1];
      return(P);
   }
   P[1]=P[1]*I[2] mod N;
   P[2]=P[2]*I[2] mod N;
   I=extgcd(Q[3],N);
   if(I[1]!=1)
   {
      P[4]=I[1];
      return(P);
   }
   Q[1]=Q[1]*I[2] mod N;
   Q[2]=Q[2]*I[2] mod N;
   if((P[1]==Q[1])&&(((P[2]+Q[2]) mod N)==0)){return(Resu);}
   bigint L;
   if((P[1]==Q[1])&&(P[2]==Q[2]))
   {
      I=extgcd(2*Q[2],N);
      if(I[1]!=1)
      {
         P[4]=I[1];
         return(P);
      }
      L=I[2]*(3*Q[1]^2+a) mod N;
   }
   else
   {
      I=extgcd(Q[1]-P[1],N);
      if(I[1]!=1)
      {
         P[4]=I[1];
         return(P);
      }
      L=(Q[2]-P[2])*I[2] mod N;
   }
   Resu[1]=(L^2-P[1]-Q[1]) mod N;
   if (Resu[1]<0) { Resu[1]=Resu[1]+N; }
   Resu[2]=(L*(P[1]-Resu[1])-P[2]) mod N;
   if (Resu[2]<0) { Resu[2]=Resu[2]+N; }
   Resu[3]=bigint(1);
   return(Resu);
}
example
{ "EXAMPLE:"; echo = 2;
   bigint N=11;
   bigint a=1;
   bigint b=6;
   list P,Q;
   P[1]=2;
   P[2]=4;
   P[3]=1;
   Q[1]=3;
   Q[2]=5;
   Q[3]=1;
   ellipticAdd(N,a,b,P,Q);
}

proc ellipticMult(bigint N, bigint a, bigint b, list P, bigint k)
"USAGE:  ellipticMult(N,a,b,P,k);
RETURN: a list L representing the point k*P
NOTE:  P=(P[1]:P[2]:P[3]) a point on the elliptic curve defined by
       y^2z=x^3+a*xz^2+b*z^3  over Z/N
EXAMPLE:example ellipticMult; shows an example
"
{
   if(P[3]==0){return(P);}
   list resu;
   resu[1]=bigint(0);
   resu[2]=bigint(1);
   resu[3]=bigint(0);

   if(k==0){return(resu);}
   if(k==1){return(P);}
   if(k==2){return(ellipticAdd(N,a,b,P,P));}
   if(k==-1)
   {
      resu=P;
      resu[2]=N-P[2];
      return(resu);
   }
   if(k<0)
   {
      resu=ellipticMult(N,a,b,P,-k);
      return(ellipticMult(N,a,b,resu,-1));
   }
   if((k mod 2)==0)
   {
      resu=ellipticMult(N,a,b,P,k div 2);
      return(ellipticAdd(N,a,b,resu,resu));
   }
   resu=ellipticMult(N,a,b,P,k-1);
   return(ellipticAdd(N,a,b,resu,P));
}
example
{ "EXAMPLE:"; echo = 2;
   bigint N=11;
   bigint a=1;
   bigint b=6;
   list P;
   P[1]=2;
   P[2]=4;
   P[3]=1;
   ellipticMult(N,a,b,P,3);
}

//================== Random for elliptic curves =====================

proc ellipticRandomCurve(bigint N)
"USAGE:  ellipticRandomCurve(N);
RETURN: a list of two random numbers a,b and 4a^3+27b^2 mod N
NOTE:   y^2z=x^3+a*xz^2+b^2*z^3 defines an elliptic curve over Z/N
EXAMPLE:example ellipticRandomCurve; shows an example
"
{
   int k;
   while(k<=10)
   {
     k++;
     bigint a=random(1,2147483647) mod N;
     bigint b=random(1,2147483647) mod N;
     //test for ellictic curve
     bigint D=4*a^3+27*b^4; //the constant term is b^2
     bigint g=gcd(D,N);
     if(g<N){return(list(a,b,g));}
   }
   ERROR("no random curve found");
}
example
{ "EXAMPLE:"; echo = 2;
   ellipticRandomCurve(32003);
}

proc ellipticRandomPoint(bigint N, bigint a, bigint b)
"USAGE:  ellipticRandomPoint(N,a,b);
RETURN: a list representing  a random point (x:y:z) of the elliptic curve
        defined by y^2z=x^3+a*xz^2+b*z^3  over Z/N
EXAMPLE:example ellipticRandomPoint; shows an example
"
{
   bigint x=random(1,2147483647) mod N;
   bigint h=x^3+a*x+b;
   h=h mod N;
   list resu;
   resu[1]=x;
   resu[2]=0;
   resu[3]=1;
   if(h==0){return(resu);}

   bigint n=Jacobi(h,N);
   if(n==0)
   {
      resu=-5;
      "N is not prime";
      return(resu);
   }
   if(n==1)
   {
      resu[2]=squareRoot(h,N);
      return(resu);
   }
   return(ellipticRandomPoint(N,a,b));
}
example
{ "EXAMPLE:"; echo = 2;
   ellipticRandomPoint(32003,3,181);
}



//====================================================================
//======== counting the points of an elliptic curve  =================
//====================================================================

//==================   the trivial approaches  =======================
proc countPoints(bigint N, bigint a, bigint b)
"USAGE:  countPoints(N,a,b);
RETURN: the number of points of the elliptic curve defined by
        y^2=x^3+a*x+b  over Z/N
NOTE: trivial approach
EXAMPLE:example countPoints; shows an example
"
{
  bigint x;
  bigint r=N+1;
  while(x<N)
  {
     r=r+Jacobi((x^3+a*x+b) mod N,N);
     x=x+1;
  }
  return(r);
}
example
{ "EXAMPLE:"; echo = 2;
   countPoints(181,71,150);
}

proc ellipticAllPoints(bigint N, bigint a, bigint b)
"USAGE:  ellipticAllPoints(N,a,b);
RETURN: list of points (x:y:z) of the elliptic curve defined by
        y^2z=x^3+a*xz^2+b*z^3  over Z/N
EXAMPLE:example ellipticAllPoints; shows an example
"
{
   list resu,point;
   point[1]=0;
   point[2]=1;
   point[3]=0;
   resu[1]=point;
   point[3]=1;
   bigint x,h,n;
   while(x<N)
   {
      h=(x^3+a*x+b) mod N;
      if(h==0)
      {
         point[1]=x;
         point[2]=0;
         resu[size(resu)+1]=point;
      }
      else
      {
         n=Jacobi(h,N);
         if(n==1)
         {
            n=squareRoot(h,N);
            point[1]=x;
            point[2]=n;
            resu[size(resu)+1]=point;
            point[2]=N-n;
            resu[size(resu)+1]=point;
         }
      }
      x=x+1;
   }
   return(resu);
}
example
{ "EXAMPLE:"; echo = 2;
   list L=ellipticAllPoints(181,71,150);
   size(L);
   L[size(L)];
}

//================ the algorithm of Shanks and Mestre =================

proc ShanksMestre(bigint q, bigint a, bigint b, list #)
"USAGE:  ShanksMestre(q,a,b); optional: ShanksMestre(q,a,b,s); s the number
         of loops in the algorithm (default s=1)
RETURN: the number of points of the elliptic curve defined by
         y^2=x^3+a*x+b  over Z/N
NOTE: algorithm of Shanks and Mestre (baby-step-giant-step)
EXAMPLE:example ShanksMestre; shows an example
"
{
   bigint n=intRoot(4*q);
   bigint m=intRoot(intRoot(16*q))+1;
   bigint d;
   int i,j,k,s;
   list B,K,T,P,Q,R,mP;
   B[1]=list(0,1,0);
   if(size(#)>0)
   {
      s=#[1];
   }
   else
   {
      s=1;
   }
   while(k<s)
   {
      P =ellipticRandomPoint(q,a,b);
      Q =ellipticMult(q,a,b,P,n+q+1);

      while(j<m)
      {
         j++;
         B[j+1]=ellipticAdd(q,a,b,P,B[j]);  //baby-step list
      }
      mP=ellipticAdd(q,a,b,P,B[j]);
      mP[2]=q-mP[2];
      while(i<m)                            //giant-step
      {
         j=0;
         while(j<m)
         {
            j++;
            if((Q[1]==B[j][1])&&(Q[2]==B[j][2])&&(Q[3]==B[j][3]))
            {
               T[1]=P;
               T[2]=q+1+n-(i*m+j-1);
               K[size(K)+1]=T;
               if(size(K)>1)
               {
                  if(K[size(K)][2]!=K[size(K)-1][2])
                  {
                     d=gcd(K[size(K)][2],K[size(K)-1][2]);
                     if(ellipticMult(q,a,b,K[size(K)],d)[3]==0)
                     {
                       K[size(K)][2]=K[size(K)-1][2];
                     }
                  }
               }
               i=int(m);
               break;
            }
         }
         i=i+1;
         Q=ellipticAdd(q,a,b,mP,Q);
      }
      k++;
   }
   if(size(K)>0)
   {
      int te=1;
      for(i=1;i<=size(K)-1;i++)
      {
         if(K[size(K)][2]!=K[i][2])
         {
           if(ellipticMult(q,a,b,K[i],K[size(K)][2])[3]!=0)
           {
             te=0;
             break;
           }
         }
      }
      if(te)
      {
         return(K[size(K)][2]);
      }
   }
   return(ShanksMestre(q,a,b,s));
}
example
{ "EXAMPLE:"; echo = 2;
   ShanksMestre(32003,71,602);
}

//==================== Schoof's algorithm =============================

proc Schoof(bigint N,bigint a, bigint b)
"USAGE:  Schoof(N,a,b);
RETURN: the number of points of the elliptic curve defined by
        y^2=x^3+a*x+b  over Z/N
NOTE:  algorithm of Schoof
EXAMPLE:example Schoof; shows an example
"
{
   int pr=printlevel;
//test for ellictic curve
   bigint D=4*a^3+27*b^2;
   bigint G=gcd(D,N);
   if(G==N){ERROR("not an elliptic curve");}
   if(G!=1){ERROR("not a prime");}

//=== small N
  // if((N<=500)&&(pr<5)){return(countPoints(int(N),a,b));}

//=== the general case
   bigint q=intRoot(4*N);
   list L=primL(2*q);
   int r=size(L);
   list T;
   int i,j;
   for(j=1;j<=r;j++)
   {
      T[j]=(testElliptic(int(N),a,b,L[j])+int(q)) mod L[j];
   }
   if(pr>=5)
   {
      "===================================================================";
      "Chinese remainder :";
      for(i=1;i<=size(T);i++)
      {
         " x =",T[i]," mod ",L[i];
      }
      "gives t+ integral part of the square root of q (to be positive)";
      chineseRem(T,L);
      "we obtain t = ",chineseRem(T,L)-q;
      "===================================================================";
   }
   bigint t=chineseRem(T,L)-q;
   return(N+1-t);
}
example
{ "EXAMPLE:"; echo = 2;
   Schoof(32003,71,602);
}

/*
needs 518 sec
Schoof(2147483629,17,3567);
2147168895
*/


proc generateG(number a,number b, int m)
"USAGE:  generateG(a,b,m);
RETURN: m-th division polynomial
NOTE: generate the so-called division polynomials, i.e., the recursively defined
polynomials p_m=generateG(a,b,m) in Z[x, y] such that, for a point (x:y:1) on the
elliptic curve defined by y^2=x^3+a*x+b  over Z/N the point@*
m*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)
and m*P=0 if and only if p_m(P)=0
EXAMPLE:example generateG; shows an example
"
{
   if(m==0){return(poly(0));}
   if(m==1){return(poly(1));}
   if(m==2){return(2*var(1));}
   if(m==3){return(3*var(2)^4+6*a*var(2)^2+12*b*var(2)-a^2);}
   if(m==4)
   {
      return(4*var(1)*(var(2)^6+5*a*var(2)^4+20*b*var(2)^3-5*a^2*var(2)^2
        -4*a*b*var(2)-8*b^2-a^3));
   }
   if((m mod 2)==0)
   {
      return((generateG(a,b,m div 2+2)*generateG(a,b,m div 2-1)^2
        -generateG(a,b,m div 2-2)*generateG(a,b,m div 2+1)^2)
        *generateG(a,b,m div 2)/(2*var(1)));
   }
   return(generateG(a,b,(m-1) div 2+2)*generateG(a,b,(m-1) div 2)^3
     -generateG(a,b,(m-1) div 2-1)*generateG(a,b,(m-1) div 2+1)^3);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,(x,y),dp;
   generateG(7,15,4);
}


static proc testElliptic(int q,bigint aa,bigint bb,int l)
"USAGE:  testElliptic(q,a,b,l);
RETURN: an integer t, the trace of the Frobenius
NOTE: the kernel for the Schoof algorithm: looks for the t such that for all
      points (x:y:1) in C[l]={P in C | l*P=0},C the elliptic curve defined by
      y^2=x^3+a*x+b  over Z/q with group structure induced by 0=(0:1:0),
      (x:y:1)^(q^2)-t*(x:y:1)^q -ql*(x:y:1)=(0:1:0), ql= q mod l, trace of
      Frobenius.
EXAMPLE:example testElliptic; shows an example
"
{
   int pr=printlevel;
   ring S=q,(y,x),(L(100000),lp);
   number a=aa;
   number b=bb;
   poly F=y2-x3-a*x-b;       // the curve C
   poly G=generateG(a,b,l);
   ideal I=std(ideal(F,G));  // the points C[l]
   poly xq=powerX(q,2,I);
   poly yq=powerX(q,1,I);
   poly xq2=reduce(subst(xq,x,xq,y,yq),I);
   poly yq2=reduce(subst(yq,x,xq,y,yq),I);
   ideal J;
   int ql=q mod l;
   if(ql==0){ERROR("q is not prime");}
   int t;
   poly F1,F2,G1,G2,P1,P2,Q1,Q2,H1,H2,L1,L2;

   if(pr>=5)
   {
      "===================================================================";
      "q=",q;
      "l=",l;
      "q mod l=",ql;
      "the Groebner basis for C[l]:";I;
      "x^q mod I = ",xq;
      "x^(q^2) mod I = ",xq2;
      "y^q mod I = ",yq;
      "y^(q^2) mod I = ",yq2;
      pause();
   }
   //==== l=2 =============================================================
   if(l==2)
   {
      xq=powerX(q,2,std(x3+a*x+b));
      J=std(ideal(xq-x,x3+a*x+b));
      if(deg(J[1])==0){t=1;}
      if(pr>=5)
      {
         "===================================================================";
         "the case l=2";
         "the gcd(x^q-x,x^3+ax+b)=",J[1];
         pause();
      }
      return(t);
   }
   //=== (F1/G1,F2/G2)=[ql](x,y) ==========================================
   if(ql==1)
   {
      F1=x;G1=1;F2=y;G2=1;
   }
   else
   {
      G1=reduce(generateG(a,b,ql)^2,I);
      F1=reduce(x*G1-generateG(a,b,ql-1)*generateG(a,b,ql+1),I);
      G2=reduce(4*y*generateG(a,b,ql)^3,I);
      F2=reduce(generateG(a,b,ql+2)*generateG(a,b,ql-1)^2
               -generateG(a,b,ql-2)*generateG(a,b,ql+1)^2,I);

   }
   if(pr>=5)
   {
      "===================================================================";
      "the point ql*(x,y)=(F1/G1,F2/G2)";
      "F1=",F1;
      "G1=",G1;
      "F2=",F2;
      "G2=",G2;
      pause();
   }
   //==== the case t=0 :  the equations for (x,y)^(q^2)=-[ql](x,y) ===
   J[1]=xq2*G1-F1;
   J[2]=yq2*G2+F2;
   if(pr>=5)
   {
      "===================================================================";
      "the case t=0 mod l";
      "the equations for (x,y)^(q^2)=-[ql](x,y) :";
      J;
      "the test, if they vanish for all points in C[l]:";
      reduce(J,I);
      pause();
   }
   //=== test if all points of C[l] satisfy  (x,y)^(q^2)=-[ql](x,y)
   //=== if so: t mod l =0 is returned
   if(size(reduce(J,I))==0){return(0);}

   //==== test for (x,y)^(q^2)=[ql](x,y) for some point

   J=xq2*G1-F1,yq2*G2-F2;
   J=std(J+I);
   if(pr>=5)
   {
      "===================================================================";
      "test if (x,y)^(q^2)=[ql](x,y) for one point";
      "if so, the Frobenius has an eigenvalue 2ql/t: (x,y)^q=(2ql/t)*(x,y)";
      "it follows that t^2=4q mod l";
      "if w is one square root of q mod l";
      "t =2w mod l or -2w mod l ";
      "-------------------------------------------------------------------";
      "the equations for (x,y)^(q^2)=[ql](x,y) :";
      xq2*G1-F1,yq2*G2-F2;
      "the test if one point satisfies them";
      J;
      pause();
   }
   if(deg(J[1])>0)
   {
      int w=int(squareRoot(q,l));
      //=== +/-2w mod l zurueckgeben, wenn (x,y)^q=+/-[w](x,y)
      //==== the case t>0 :  (Q1/P1,Q2/P2)=[w](x,y) ==============
      if(w==1)
      {
         Q1=x;P1=1;Q2=y;P2=1;
      }
      else
      {
         P1=reduce(generateG(a,b,w)^2,I);
         Q1=reduce(x*G1-generateG(a,b,w-1)*generateG(a,b,w+1),I);
         P2=reduce(4*y*generateG(a,b,w)^3,I);
         Q2=reduce(generateG(a,b,w+2)*generateG(a,b,w-1)^2
               -generateG(a,b,w-2)*generateG(a,b,w+1)^2,I);
      }
      J=xq*P1-Q1,yq*P2-Q2;
      J=std(I+J);
      if(pr>=5)
      {
      "===================================================================";
      "the Frobenius has an eigenvalue, one of the roots of  w^2=q mod l:";
      "one root is:";w;
      "test, if it is the eigenvalue (if not it must be -w):";
      "the equations for (x,y)^q=w*(x,y)";I;xq*P1-Q1,yq*P2-Q2;
      "the Groebner basis";
       J;
         pause();
      }
      if(deg(J[1])>0){return(2*w mod l);}
      return(-2*w mod l);
   }

   //==== the case t>0 :  (Q1/P1,Q2/P2)=(x,y)^(q^2)+[ql](x,y) =====
   P1=reduce(G1*G2^2*(F1-xq2*G1)^2,I);
   Q1=reduce((F2-yq2*G2)^2*G1^3-F1*G2^2*(F1-xq2*G1)^2-xq2*P1,I);
   P2=reduce(P1*G2*(F1-xq2*G1),I);
   Q2=reduce((xq2*P1-Q1)*(F2-yq2*G2)*G1-yq2*P2,I);

   if(pr>=5)
   {
      "we are in the general case:";
      "(x,y)^(q^2)!=ql*(x,y) and (x,y)^(q^2)!=-ql*(x,y) ";
      "the point (Q1/P1,Q2/P2)=(x,y)^(q^2)+[ql](x,y)";
      "Q1=",Q1;
      "P1=",P1;
      "Q2=",Q2;
      "P2=",P2;
      pause();
   }
   while(t<(l-1) div 2)
   {
      t++;
      //====  (H1/L1,H2/L2)=[t](x,y)^q ===============================
      if(t==1)
      {
         H1=xq;L1=1;
         H2=yq;L2=1;
      }
      else
      {
         H1=x*generateG(a,b,t)^2-generateG(a,b,t-1)*generateG(a,b,t+1);
         H1=subst(H1,x,xq,y,yq);
         H1=reduce(H1,I);
         L1=generateG(a,b,t)^2;
         L1=subst(L1,x,xq,y,yq);
         L1=reduce(L1,I);
         H2=generateG(a,b,t+2)*generateG(a,b,t-1)^2
           -generateG(a,b,t-2)*generateG(a,b,t+1)^2;
         H2=subst(H2,x,xq,y,yq);
         H2=reduce(H2,I);
         L2=4*y*generateG(a,b,t)^3;
         L2=subst(L2,x,xq,y,yq);
         L2=reduce(L2,I);
      }
      J=Q1*L1-P1*H1,Q2*L2-P2*H2;
      if(pr>=5)
      {
      "we test now the different t, 0<t<=(l-1)/2:";
      "the point (H1/L1,H2/L2)=[t](x,y)^q :";
      "H1=",H1;
      "L1=",L1;
      "H2=",H2;
      "L2=",L2;
      "the equations for (x,y)^(q^2)+[ql](x,y)=[t](x,y)^q :";J;
      "the test";reduce(J,I);
      "the test for l-t (the x-cordinate is the same):";
       Q1*L1-P1*H1,Q2*L2+P2*H2;
       reduce(ideal(Q1*L1-P1*H1,Q2*L2+P2*H2),I);
         pause();
      }
      if(size(reduce(J,I))==0){return(t);}
      J=Q1*L1-P1*H1,Q2*L2+P2*H2;
      if(size(reduce(J,I))==0){return(l-t);}
   }
   ERROR("something is wrong in testElliptic");
}
example
{ "EXAMPLE:"; echo = 2;
   testElliptic(1267985441,338474977,64740730,3);
}

//============================================================================
//================== Factorization and Primality Test ========================
//============================================================================

//============= Lenstra's ECM Factorization ==================================

proc factorLenstraECM(bigint N, list S, int B, list #)
"USAGE:  factorLenstraECM(N,S,B); optional: factorLenstraECM(N,S,B,d);
         d+1 the number of loops in the algorithm (default d=0)
RETURN: a factor of N or the message no factor found
NOTE: - computes a factor of N using Lenstra's ECM factorization@*
      - the idea is that the fact that N is not prime is dedected using
        the operations on the elliptic curve
      - is similarly to Pollard's p-1-factorization
EXAMPLE:example factorLenstraECM; shows an example
"
{
   list L,P;
   bigint g,M,w;
   int i,j,k,d;
   int l=size(S);
   if(size(#)>0)
   {
      d=#[1];
   }

   while(i<=d)
   {
      L=ellipticRandomCurve(N);
      if(L[3]>1){return(L[3]);} //the discriminant was not invertible
      P=list(0,L[2],1);
      j=0;
      M=1;
      while(j<l)
      {
         j++;
         w=S[j];
         if(w>B) break;
         while(w*S[j]<B)
         {
           w=w*S[j];
         }
         M=M*w;
         P=ellipticMult(N,L[1],L[2]^2,P,w);
         if(size(P)==4){return(P[4]);}  //some inverse did not exsist
         if(P[3]==0){break;}            //the case M*P=0
      }
      i++;
   }
   return("no factor found");
}
example
{ "EXAMPLE:"; echo = 2;
   list L=primList(1000);
   factorLenstraECM(181*32003,L,10,5);
   bigint h=10;
   h=h^30+25;
   factorLenstraECM(h,L,4,3);
}

//================= ECPP (Goldwasser-Kilian) a primaly-test =============

proc ECPP(bigint N)
"USAGE:  ECPP(N);
RETURN: message:N is not prime or {L,P,m,q} as certificate for N being prime@*
         L a list (y^2=x^3+L[1]*x+L[2] defines an elliptic curve C)@*
         P a list ((P[1]:P[2]:P[3]) is a point of C)@*
         m,q integers
ASSUME: gcd(N,6)=1
NOTE:   The basis of the algorithm is the following theorem:
         Given C, an elliptic curve over Z/N, P a point of C(Z/N),
         m an integer, q a prime with the following properties:
         - q|m
         - q>(4-th root(N) +1)^2
         - m*P=0=(0:1:0)
         - (m/q)*P=(x:y:z) and z a unit in Z/N
         Then N is prime.
EXAMPLE:example ECPP; shows an example
"
{
   list L,S,P;
   bigint m,q;
   int i;

   bigint n=intRoot(intRoot(N));
   n=(n+1)^2;                         //lower bound for q
   while(1)
   {
      L=ellipticRandomCurve(N);       //a random elliptic curve C
      m=ShanksMestre(N,L[1],L[2],3);  //number of points of the curve C
      S=PollardRho(m,10000,1);        //factorization of m
      for(i=1;i<=size(S);i++)         //search for q between the primes
      {
         q=S[i];
         if(n<q){break;}
      }
      if(n<q){break;}
   }
   bigint u=m/q;
   while(1)
   {
      P=ellipticRandomPoint(N,L[1],L[2]);  //a random point on C
      "P=",P;
      if(ellipticMult(N,L[1],L[2],P,m)[3]!=0){"N is not prime";return(-5);}
      if(ellipticMult(N,L[1],L[2],P,u)[3]!=0)
      {
         L=delete(L,3);
         return(list(L,P,m,q));
      }
   }
}
example
{ "EXAMPLE:"; echo = 2;
   bigint N=1267985441;
   ECPP(N);
}

static proc wordToNumber(string s)
{
   int i;
   intvec v;
   bigint n;
   bigint t=27;
   for(i=size(s);i>0;i--)
   {
      if(s[i]=="a"){v[i]=0;}
      if(s[i]=="b"){v[i]=1;}
      if(s[i]=="c"){v[i]=2;}
      if(s[i]=="d"){v[i]=3;}
      if(s[i]=="e"){v[i]=4;}
      if(s[i]=="f"){v[i]=5;}
      if(s[i]=="g"){v[i]=6;}
      if(s[i]=="h"){v[i]=7;}
      if(s[i]=="i"){v[i]=8;}
      if(s[i]=="j"){v[i]=9;}
      if(s[i]=="k"){v[i]=10;}
      if(s[i]=="l"){v[i]=11;}
      if(s[i]=="m"){v[i]=12;}
      if(s[i]=="n"){v[i]=13;}
      if(s[i]=="o"){v[i]=14;}
      if(s[i]=="p"){v[i]=15;}
      if(s[i]=="q"){v[i]=16;}
      if(s[i]=="r"){v[i]=17;}
      if(s[i]=="s"){v[i]=18;}
      if(s[i]=="t"){v[i]=19;}
      if(s[i]=="u"){v[i]=20;}
      if(s[i]=="v"){v[i]=21;}
      if(s[i]=="w"){v[i]=22;}
      if(s[i]=="x"){v[i]=23;}
      if(s[i]=="y"){v[i]=24;}
      if(s[i]=="z"){v[i]=25;}
      if(s[i]==" "){v[i]=26;}
   }
   for(i=1;i<=size(s);i++)
   {
      n=n+v[i]*t^(i-1);
   }
   return(n);
}

static proc numberToWord(bigint n)
{
   int i,j;
   string v;
   list s;
   bigint t=27;
   bigint mm;
   bigint nn=n;
   while(nn>t)
   {
      j++;
      mm=nn mod t;
      s[j]=mm;
      nn=(nn-mm) div t;
   }
   j++;
   s[j]=nn;
   for(i=1;i<=j;i++)
   {
      if(s[i]==0){v=v+"a";}
      if(s[i]==1){v=v+"b";}
      if(s[i]==2){v=v+"c";}
      if(s[i]==3){v=v+"d";}
      if(s[i]==4){v=v+"e";}
      if(s[i]==5){v=v+"f";}
      if(s[i]==6){v=v+"g";}
      if(s[i]==7){v=v+"h";}
      if(s[i]==8){v=v+"i";}
      if(s[i]==9){v=v+"j";}
      if(s[i]==10){v=v+"k";}
      if(s[i]==11){v=v+"l";}
      if(s[i]==12){v=v+"m";}
      if(s[i]==13){v=v+"n";}
      if(s[i]==14){v=v+"o";}
      if(s[i]==15){v=v+"p";}
      if(s[i]==16){v=v+"q";}
      if(s[i]==17){v=v+"r";}
      if(s[i]==18){v=v+"s";}
      if(s[i]==19){v=v+"t";}
      if(s[i]==20){v=v+"u";}
      if(s[i]==21){v=v+"v";}
      if(s[i]==22){v=v+"w";}
      if(s[i]==23){v=v+"x";}
      if(s[i]==24){v=v+"y";}
      if(s[i]==25){v=v+"z";}
      if(s[i]==26){v=v+" ";}
   }
   return(v);
}

proc code(string s)
"USAGE:  code(s); s a string
ASSUME:  s contains only small letters and space
COMPUTE: a bigint, RSA-coding of the string s
RETURN:  return RSA-coding of the string s as string
EXAMPLE: code;  shows an example
"
{
   ring r=0,x,dp;
   bigint
p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
   bigint
q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
   bigint n=p*q;
   bigint phi=(p-1)*(q-1);
   bigint e=1234567891;
   //bigint d=extgcd(e,phi)[2];
   bigint m=wordToNumber(s);
   bigint c=powerN(m,e,n);
   string cc=string(c);
   return(cc);
}
example
{"EXAMPLE:";  echo = 2;
  string s="i go to school";
  code(s);
}

proc decodeString(string g)
"USAGE:  decodeString(s); s a string
ASSUME:  s is a string of a bigint, the output of code
COMPUTE: a string, RSA-decoding of the string s
RETURN:  return RSA-decoding of the string s as string
EXAMPLE: decodeString;  shows an example
"
{
   bigint
p=398075086424064937397125500550386491199064362342526708406385189575946388957261768583317;
   bigint
q=472772146107435302536223071973048224632914695302097116459852171130520711256363590397527;
   bigint n=p*q;
   bigint phi=(p-1)*(q-1);
   bigint e=1234567891;
   bigint d=extgcd(e,phi)[2];
   execute("bigint c="+g+";");
   bigint f=powerN(c,d,n);
   string s=numberToWord(f);
   return(s);
}
example
{"EXAMPLE:";  echo = 2;
  string
s="78638618599886548153321853785991541374544958648147340831959482696082179852616053583234149080198937632782579537867262780982185252913122030800897193851413140758915381848932565";
  string t=decodeString(s);
  t;
}
/*----------------------------------------------------------------------------
 * set stuff
 * -------------------------------------------------------------------------*/
static proc set_multiply_list_content(list h)
"USAGE:   set_multiply_list_content(h)
RETURN:  An integer c als product of all elements in h
EXAMPLE: example set_multiply_list_content; shows an example;
"
{
  int c = 1;
  for (int i=1;i<=size(h);i++)
  {
    c = c*h[i];
  }
  return(c);
}
example
{
  "EXAMPLE:"; echo = 2;
  list h=2,4,5;
  set_multiply_list_content(h);
}

static proc set_delete_certain_element(list a, int e)
"USAGE:   set_delete_certain_element(a,e)
RETURN:  A list a without element e. If e was not in the list before, a will not be changed
EXAMPLE: example set_delete_certain_element; shows an example.
"
{
  list output_list = a;
  for (int i=1;i<=size(a);i++)
  {
    if (a[i]==e)
    {
      output_list = delete(output_list,i);
    }
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list h=2,4,5;
  set_delete_certain_element(h,4);
  set_delete_certain_element(h,10);
}

static proc set_bubblesort_int(list output_list)
"USAGE:   set_bubblesort_int(output_list)
RETURN:  An ascending sorted list with integer values.
EXAMPLE: set_bubblesort_int; shows an example.
"
{
  output_list = bubblesort(output_list);
  //Cast every value into an integer
  for (int i=1; i<=size(output_list);i++)
  {
    output_list[i] = int(output_list[i]);
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list output_list=10,4,5,24,9;
  set_bubblesort_int(output_list);
}

static proc set_is_set(list a)
"USAGE:   set_is_set(a)
RETURN:  1 if the list is a set, 0 the list contains any duplicated elements
EXAMPLE: set_is_set; shows an example.
"
{
  int i,v;
  for (v=1; v<=size(a); v++)
  {
    for (i=v+1; i<=size(a); i++)
    {
      if (a[i]==a[v])
      {
        return(0);
      }
    }
  }
  return(1);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set = 1,5,10,2;
  list noset = 1,5,10,2,5;
  set_is_set(set);
  set_is_set(noset);
}

static proc set_contains(list a, int e)
"USAGE:   set_contains(a,e)
RETURN:  1 if the list contains e, 0 otherwise
EXAMPLE: set_contains; shows an example.
"
{
  for (int v=1; v<=size(a); v++)
  {
    if (a[v]==e)
    {
      return(1);
    }
  }
  return(0);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set = 1,5,10,2;
  set_contains(set,5);
  set_contains(set,40);
}

static proc set_delete_duplicates(list a)
"USAGE:   set_delete_duplicates(a)
RETURN:  a list a without any duplicated elements
EXAMPLE: set_delete_duplicates; shows an example.
"
{
  int i;
  list output_list = a[1];
  for (i=1; i<=size(a); i++)
  {
    if (set_contains(output_list,a[i])==0)
    {
      output_list = insert(output_list,a[i]);
    }
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set = 1,5,10,2,10,5;
  set_delete_duplicates(set);
}

static proc set_equals(list a,list b)
"USAGE:   set_equals(a, b)
RETURN:  1 if the lists are equal from a set-structure standpoint, 0 otherwise
EXAMPLE: set_equals; shows an example.
"
{
  //Checks if the lists have the same length
  if (size(a)!=size(b))
  {
    return(0);
  }

  //Sorts the lists
  a = set_bubblesort_int(a);
  b = set_bubblesort_int(b);

  //Checks every single element of both lists
  for (int i=1; i<=size(a); i++)
  {
    if (a[i]!=b[i])
    {
      return(0);
    }
  }
  return(1);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set1 = 1,5,10,2;
  list set2 = 10,2,5,1;
  list set3 = 1,5,9,2;
  set_equals(set1,set2);
  set_equals(set1,set3);
}

static proc set_insert(list a, int e)
"USAGE:   set_insert(a,e)
RETURN:  list a containing element e
EXAMPLE: set_insert; shows an example.
"
{
  if(set_contains(a,e))
  {
    return(a);
  }
  else
  {
    a=insert(a,e);
    return(a);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  list set = 1,5,10,2;
  set_insert(set,5);
  set_insert(set,22);
}

static proc set_union(list a, list b)
"USAGE:   set_union(a, b)
RETURN:  list a as union of a and b
EXAMPLE: set_union; shows an example.
"
{
  for (int i=1; i<=size(b); i++)
  {
    if (set_contains(a,b[i])==0)
    {
      a = insert(a,b[i]);
    }
  }
  return(a);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set1 = 1,5,10,2;
  list set2 = 5,10,93,58,29;
  set_union(set1,set2);
}

static proc set_section(list a, list b)
"USAGE:   set_section(a, b)
RETURN:  list output_list as intersection of a and b
EXAMPLE: set_section; shows an example.
"
{
  list output_list;
  for (int i=1; i<=size(a); i++)
  {
    if (set_contains(b,a[i])==1)
    {
      output_list = insert(output_list,a[i]);
    }
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set1 = 1,5,10,2;
  list set2 = 5,10,93,58,29;
  set_section(set1,set2);
}

static proc set_list_delete_duplicates(list a)
"USAGE:   set_list_delete_duplicates(a)
RETURN:  list output_list with no duplicated lists
EXAMPLE: set_list_delete_duplicates; shows an example.
"
{
  int v;
  int i;
  int counter;
  int out_size;
  list output_list=insert(output_list,a[1]);
  //Create a new list and try to insert every list element from a into that list. If a list is already inserted into the
  //new list, do nothing.
  for (i=2; i<=size(a); i++)
  {
    out_size = size(output_list);
    counter = 0;

    for (v=1; v<=out_size; v++)
    {
      if (set_equals(output_list[v],a[i])==1)
      {
        counter++;
      }
    }
    if (counter==0)
    {
      output_list = insert(output_list,a[i]);
    }
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set1 = 1,5,10,2;
  list set2 = 1,10,2,5;
  list set3 = 1,2,3;
  list superset = set1,set2,set3;
  set_list_delete_duplicates(superset);
}

static proc set_construct_increasing_set(int maxelement)
"USAGE:   set_construct_increasing_set(maxelement)
RETURN:  list output_list with increasing elements from 1 to maxelement
EXAMPLE: set_construct_increasing_set; shows an example.
"
{
  list output_list;
  for (int i=1; i<=maxelement; i++)
  {
    output_list = insert(output_list, i);
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  set_construct_increasing_set(10);
}

static proc set_addtoall(list a, int element)
"USAGE:   set_addtoall(a,alement)
RETURN:  Transformed list with e_i+element for every element e_i in list a
EXAMPLE: set_addtoall; shows an example.
"
{
  list output_list;
  int c;
  for (int i=1; i<=size(a); i++)
  {
    c = a[i]+element;
    output_list = insert(output_list,c);
  }
  return(set_turn(output_list));
}
example
{
  "EXAMPLE:"; echo = 2;
  list set = 1,5,10,2;
  set_addtoall(set,5);
}

static proc set_turn(list a)
"USAGE:   set_turn(a)
RETURN:  Turned list a
EXAMPLE: set_turn; shows an example.
"
{
  list output_list;
  for (int i=1; i<=size(a); i++)
  {
    output_list[size(a)+1-i] = a[i];
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set = 1,5,10;
  set_turn(set);
}

static proc set_subset_set(list a)
"USAGE:   set_subset_set(a)
RETURN:  Set of subsets
EXAMPLE: set_subset_set; shows an example.
"
{
  int v;
  int choice;
  list output_list;
  list start = a[1];
  list insertion_list;
  int output_length;
  output_list = insert(output_list,start);

  for (int i=2; i<=size(a); i++)
  {
    choice = 0;
    start = a[i];
    output_list = insert(output_list,start);
    output_length = size(output_list);
    for (v=2; v<=output_length; v++)
    {
      insertion_list = set_insert(output_list[v],a[i]);
      output_list = insert(output_list, insertion_list,size(output_list));
    }
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list set = 1,5,10;
  set_subset_set(set);
}
/* -----------------------------------------------------------------
 * knapsack_utilities: Utility functions needed for knapsack
 * -----------------------------------------------------------------*/
proc calculate_ordering(bigint num1, bigint primitive, bigint mod1)
"USAGE:   calculate_ordering(num1, primitive, mod1)
RETURN:  x so that primitive^x == num1 mod mod1
EXAMPLE: example calculate_ordering; shows an example;
"
{
  for (int i=1;i<=int((mod1-2));i++)
  {
    if ((primitive^i%mod1)==num1)
    {
      return(i);
    }
  }
  return(0);
}
example
{
  "EXAMPLE:"; echo = 2;
  bigint mod1 = 33;
  bigint primitive = 14;
  bigint num1 = 5;
  calculate_ordering(num1,primitive,mod1);
}

proc is_primitive_root(bigint primitive, bigint mod1)
"USAGE:   is_primitive_root(primitive, mod1)
RETURN:  1 if primitive is a primitive root modulo mod1, 0 otherwise
EXAMPLE: example is_primitive_root; shows an example;
"
{
  list output_list;
  for (int i=1;i<=int((mod1-1));i++)
  {
    output_list = set_insert(output_list,int((primitive^i%mod1)));
  }
  if (bigint(size(output_list))==bigint(mod1-1))
  {
    return(1);
  }
  else
  {
    return(0);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  is_primitive_root(3,7);
  is_primitive_root(2,7);
}

proc find_first_primitive_root(bigint mod1)
"USAGE:   find_first_primitive_root(mod1)
RETURN:  First primitive root modulo mod1, 0 if no root can be found.
EXAMPLE: example find_first_primitive_root; shows an example;
"
{
  for (int i=0;i<=int(mod1-1);i++)
  {
    if (is_primitive_root(bigint(i),mod1)==1)
    {
      return(i);
    }
  }
  return(0);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x,lp;
  find_first_primitive_root(7);
  find_first_primitive_root(557);
}

proc binary_add(list binary_list)
"USAGE:   binary_add(binary_list)
RETURN:  binary encoded list, increased by 1
EXAMPLE: example binary_add; shows an example;
"
{
  int residual=1;
  int position = size(binary_list);
  while((residual==1)&&(position!=0))
  {
    if(binary_list[position]==0)
    {
      binary_list[position]=1;
      residual=0;
    }
    else
    {
      binary_list[position]=0;
      position--;
    }
  }
  if (position==0)
  {
    binary_list = insert(binary_list,1);
  }
  return(binary_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x,lp;
  list binary_list = 1,0,1,1,1;
  binary_add(binary_list);
}

proc inverse_modulus(int num, int mod1)
"USAGE:   inverse_modulus(num, mod1)
RETURN:  inverse element of num modulo mod1
EXAMPLE: example inverse_modulus; shows an example;
"
{
  if (num>=mod1)
  {
    return(0);
  }
  else
  {
    for (int i=1;i<mod1;i++)
    {
      if ((i*num%mod1)==1)
      {
        return(i);
      }
    }
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x,lp;
  int mod1 = 13;
  int num = 5;
  inverse_modulus(num,mod1);
}

proc is_prime(int n)
"USAGE:   is_prime(n)
RETURN:  1 if n is prime, 0 otherwise
EXAMPLE: example is_prime; shows an example;
"
{
  int prime1 = 1;
  for (int i=n-1;i>1;i--)
  {
    if(n%i==0)
    {
        prime1 = 0;
    }
  }
  return(prime1);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x,lp;
  is_prime(10);
  is_prime(7);
}

proc find_biggest_index(list a)
"USAGE:   find_biggest_index( a)
RETURN:  Returns the index of the biggest element of a
EXAMPLE: example find_biggest_index; shows an example;
"
{
  list sortedlist = bubblesort(a);
  return(find_index(a,sortedlist[1]));
}

proc find_index(list a, bigint e)
"USAGE:   find_index(a, e)
RETURN:  Returns the list index of element e in list a. Returns 0 if e is not in a
EXAMPLE: example find_index; shows an example;
"
{
  for(int i=1;i<=size(a);i++)
  {
    if (bigint(a[i])==e)
    {
      return(i);
    }
  }
  return(0);
}
example
{
  "EXAMPLE:"; echo = 2;
  list a = 1,5,20,6,37;
  find_index(a,20);
  find_index(a,6);
  find_index(a,100);
}
/* ------------------------------------------------------------------
 * Knapsack Algorithmus such as solving several knapsack problems,
 * kryptographic algorithms and algorithms for creatings suitable knapsacks
 * ---------------------------------------------------------------- */
proc subset_sum01(list knapsack, int solution)
"USAGE:   subset_sum01(knapsack,solution)
RETURN:  binary list of the positions of the elements included in the subset sum or 0 if no solution exists
NOTE: This will return the first solution of the ssk-problem, given be the smallest binary encoding. It wont return several solutions if they exist
EXAMPLE: example subset_sum01; shows an example;
"
{
  int i;
  int v;
  int comparable;
  //Check if the knapsack is a set
  if (set_is_set(knapsack)==1)
  {
    //Create a binary list full of zeroes
    list binary_list;
    for (i=1;i<=size(knapsack);i++)
    {
      binary_list = insert(binary_list,0);
    }
    binary_list = binary_add(binary_list);

    for(i=1;i<=2^(size(knapsack));i++)
    {
      comparable = 0;
      //Create the Subset-Sum for the actual binary coding of binary_list
      for (v=1;v<=size(knapsack);v++)
      {
        comparable = comparable+knapsack[v]*binary_list[v];
      }
      //Check if the sum equals the solution
      if (comparable==solution)
      {
        return(binary_list);
      }
      else
      {
        binary_list = binary_add(binary_list);
      }
    }
    return(0);
  }
  else
  {
    return(0);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  list h=1,4,7,32;
  subset_sum01(h,20);
  subset_sum01(h,11);
  subset_sum01(h,33);
}

proc subset_sum02(list knapsack, int sol)
"USAGE:   subset_sum02(knapsack,sol)
RETURN:  binary list of the positions of the elements included in the subset sum or 0 if no solution exists
EXAMPLE: example subset_sum02; shows an example;
"
{
  list summands;
  int calcu;
  int i;
  //Create a sorted copy of the knapsack, calling it worksack
  list worksack = set_bubblesort_int(knapsack);
  int counter = 1;
  while((counter<=size(worksack))&&(sol>0))
  {
    //Try to substract an element of the knapsack from the capacity. Create a list with all the summands used
    calcu = sol-worksack[counter];
    if (calcu>=0)
    {
      sol = sol-worksack[counter];
      summands = insert(summands,int(worksack[counter]));
    }
    counter++;
  }
  if(sol>0)
  {
    return(0);
  }

  //Get the index of the summands of the original knapsack and change the bits in the binary list wich will be the solution
  list binary_list;
  for (i=1;i<=size(knapsack);i++)
  {
    binary_list = insert(binary_list,0);
  }

  for (i=1; i<=size(knapsack);i++)
  {
    if (set_contains(summands,knapsack[i])==1)
    {
      binary_list[i]=1;
    }
  }
  return(binary_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list h=1,4,7,32;
  subset_sum02(h,20);
  subset_sum02(h,11);
  subset_sum02(h,33);
}

proc unbounded_knapsack(list knapsack, list profit, int capacity)
"USAGE:   unbounded_knapsack(knapsack,profit,capacity)
RETURN:  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.
EXAMPLE: unbounded_knapsack; shows an example;
"
{
  int i;
  int v;
  list output_list;
  for (i=1;i<=capacity+1;i++)
  {
    output_list = insert(output_list,0);
  }
  for (i=1;i<=capacity+1;i++)
  {
    for(v=1;v<=size(knapsack);v++)
    {
      if (knapsack[v]<i)
      {
        if(output_list[i]<(output_list[i-knapsack[v]]+profit[v]))
        {
          output_list[i] = output_list[i-knapsack[v]]+profit[v];
        }
      }
    }
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list h=1,4,7,32;
  list knapsack = 5,2;
  list profit = 10,3;
  int capacity = 5;
  unbounded_knapsack(knapsack,profit,capacity);
}

proc multidimensional_knapsack(matrix m, list capacities, list profits)
"USAGE:   multidimensional_knapsack(m,capacities,profits)
RETURN:  binary list of the positions of the elements included in the optimal selection
EXAMPLE: example multidimensional_knapsack; shows an example;
"
{
  int index;
  list output_list;
  list nolist;
  list y_list;
  list minmax_list;
  int i;
  int v;
  int checkint;
  //create the output_list full of zeroes with the length of all given selections
  for (i=1;i<=size(profits);i++)
  {
    output_list = insert(output_list,0);
  }

  //Create the List E with all indices of the output_list that haven't been used yet.
  list E = set_turn(set_construct_increasing_set(size(profits)));

  //Repeat till every index in E is used.
  while(size(E)>0)
  {
    y_list = nolist;
    for (i=1;i<=size(E);i++)
    {
      //Create all possible elements of y_i (y_list). minimax_list will be replaced of an empty list in each iteration
      minmax_list = nolist;
      for (v=1; v<=size(capacities);v++)
      {
        if (set_contains(E,v)==1)
        {
        minmax_list = insert(minmax_list, bigint(capacities[v])/bigint(m[v,E[i]]));
        }
      }
      //Sort the elements so that it is easy to pick the smallest one
      minmax_list = bubblesort(minmax_list);

      //insert Element y_i into y_list, wich is the smallest of (b_i/w_ij) like in the PECH algorithm description.
      y_list = insert(y_list,minmax_list[size(minmax_list)],size(y_list));


    }

    //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.
    checkint=0;
    for(i=1;i<=size(y_list);i++)
    {
      if (y_list[i]>=1)
      {
        checkint=1;
      }
    }
    if (checkint==0)
    {
      return(output_list);
    }


    //Find the index of the selection and update the binary output_list
    minmax_list = nolist;
    for (i=1;i<=size(E);i++)
    {
      minmax_list = insert(minmax_list, profits[E[i]]*y_list[i],size(minmax_list));
    }
    index = find_biggest_index(minmax_list);

    output_list[E[index]]=1;

    //Update the capacities by substracting the weights of the selection
    for (i=1;i<=size(capacities);i++)
    {
      capacities[i] = capacities[i]- m[i,E[index]];
    }
    E = set_delete_certain_element(E,index);

  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x,lp;
  matrix m[3][3] = 1,4,10,7,8,3,1,9,7;
  list c = 12,17,10;
  list p = 3,2,5;
  multidimensional_knapsack(m,c,p);
}

proc naccache_stern_generation(int key, int primenum)
"USAGE:   naccache_stern_generation(key, primenum)
RETURN:  a hard knapsack list
EXAMPLE: example naccache_stern_generation; shows an example;
"
{
  //Check if primenum is a prime and the gcd-Condition holds
  if ((is_prime(primenum)==0)||(gcd(key,(primenum-1))!=1))
  {
    return(0);
  }
  else
  {
    int i;
    int p;
    list primelist;
    int primecounter=2;
    //Generate the knapsack containing the smallest prime numbers so that primenum exceeds the product of all of them
    while(set_multiply_list_content(primelist)<primenum)
    {
      if (is_prime(primecounter)==1)
      {
      primelist = insert(primelist,primecounter);
      }
      primecounter++;
    }
    primelist = delete(primelist,1);


    //Generate the hard knapsack of the length of the prime numbers containing zeroes.
    list hardknapsack;
    for (i=1;i<=size(primelist);i++)
    {
      hardknapsack = insert(hardknapsack,0);
    }

    //Create the elements of the hard knapsack
    primecounter = 1;
    bigint calcu;
    i=0;
    while (i<size(primelist))
    {
      //Create some v_i^key%primenum and store it in calcu
      calcu = bigint(primecounter)^key%bigint(primenum);

      //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
      if(set_contains(primelist,int(calcu))==1)
      {
        p=find_index(primelist,int(calcu));
        hardknapsack[p] = primecounter;
        i++;
      }
      primecounter++;
    }
  return(hardknapsack);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  naccache_stern_generation(5,292);
  naccache_stern_generation(5,293);
}

proc naccache_stern_encryption(list knapsack, list message, int primenum)
"USAGE:   naccache_stern_encryption(knapsack, message, primenum)
RETURN:  an encrypted message as integer
EXAMPLE: example naccache_stern_encryption; shows an example;
"
{
  bigint solution = 1;
  if (size(knapsack)==size(message))
  {
    for(int i=1;i<=size(knapsack);i++)
    {
      solution = solution*((bigint(knapsack[i])^message[i])%bigint(primenum));
    }
    return(solution);
  }
  else
  {
    return(0);
  }

}
example
{
  "EXAMPLE:"; echo = 2;
  //Please note that the values for primenum and hardknapsack have been obtained from the example of naccache_stern_generation!
  list hardknapsack = 85,164,117,44;
  int primenum = 293;
  list message = 1,0,1,0;
  naccache_stern_encryption(hardknapsack,message,primenum);
}

proc naccache_stern_decryption(list knapsack, int key, int primenum, int message)
"USAGE:   naccache_stern_decryption(knapsack, key, primenum, message)
RETURN:  decrypted binary list
EXAMPLE: example naccache_stern_decryption; shows an example;
"
{
  //create a binary list with the length of the knapsack containing zeros
  int k = int(bigint(message)^key%bigint(primenum));
  int i;
  list binary_list;
  for (i=1;i<=size(knapsack);i++)
  {
    binary_list = insert(binary_list,0);
  }

  //create primelist like in int naccache_stern_generation process
  list primelist;
  int primecounter=2;
  while(size(primelist)<size(knapsack))
  {
    if (is_prime(primecounter)==1)
    {
      primelist = insert(primelist,primecounter);
    }
      primecounter++;
  }

  //find divisors of k and update the binarylist in a way that the positions of the divisors in primelist are marked
  for (i=1;i<=size(primelist);i++)
  {
    if(k%primelist[i]==0)
    {
      binary_list[i]=1;
    }
  }
  return(binary_list);

}
example
{
  "EXAMPLE:"; echo = 2;
  //Please note that the values have been obtained from the example of naccache_stern_generation and naccache_stern_encryption!
  int primenum = 293;
  int message = 9945;
  int key = 5;
  list hardknapsack = 85,164,117,44;
  naccache_stern_decryption(hardknapsack,key,primenum,message);
}

proc m_merkle_hellman_transformation(list knapsack, int primitive, int mod1)
"USAGE:   m_merkle_hellman_transformation(knapsack, primitive, mod1)
RETURN:  list containing a hard knapsack
EXAMPLE: example m_merkle_hellman_transformation; shows an example;
"
{
  bigint new_element;
  list output_list;
  //calculate the primitiv root of every element in knapsack and insert it into a new knapsack
  for (int i=size(knapsack);i>=1;i--)
  {
    new_element = calculate_ordering(knapsack[i],primitive,mod1);
    output_list = insert(output_list,int(new_element));
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  //Please note that the values for primenum and hardknapsack have been obtained from the example of naccache_stern_generation and naccache_stern_encryption!
  list knapsack = 2,3,5,7;
  int mod1 = 211;
  int primitive = 2;
  m_merkle_hellman_transformation(knapsack,primitive,mod1);
}

proc m_merkle_hellman_encryption(list knapsack, list message)
"USAGE:    m_merkle_hellman_encryption(knapsack, message)
RETURN:  an encrypted message as integer
NOTE: 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.
EXAMPLE: example m_merkle_hellman_encryption; shows an example;
"
{
  return(merkle_hellman_encryption(knapsack,message));
}
example
{
  "EXAMPLE:"; echo = 2;
  //Please note that the values for primenum and hardknapsack have been obtained from the example of m_merkle_hellman_transformation!
  list knapsack = 1,43,132,139;
  list message = 1,0,0,1;
  m_merkle_hellman_encryption(knapsack,message);
}

proc m_merkle_hellman_decryption(list knapsack, bigint primitive, bigint mod1, int message)
"USAGE:  m_merkle_hellman_decryption(knapsack, primitive, mod1, message)
RETURN:  decrypted binary list
EXAMPLE: example merkle_hellman_decryption; shows an example;
"
{
  //Convert message
  int factorizing = int((primitive^message)%mod1);
  int i;

  //Create binary list of length of the knapsack, containing zeroes.
  list binary_list;
  for (i=1;i<=size(knapsack);i++)
  {
    binary_list = insert(binary_list,0);
  }

  //factorize the converted message, mark the factor positions in knapsack as bits in binary_list
  for (i=1;i<=size(knapsack);i++)
  {
    if(factorizing%knapsack[i]==0)
    {
      binary_list[i]=1;
    }
  }
  return(binary_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  //Please note that the values  have been obtained from the example of m_merkle_hellman_encryption and m_merkle_hellman_transformation!
  list knapsack = 2,3,5,7;
  int message = 140;
  bigint primitive = 2;
  bigint mod1 = 211;
  m_merkle_hellman_decryption(knapsack,primitive,mod1,message);
}

proc merkle_hellman_transformation(list knapsack, int key, int mod1)
"USAGE:   merkle_hellman_transformation(knapsack, key, mod1)
RETURN:  hard knapsack
EXAMPLE: example merkle_hellman_transformation; shows an example;
"
{
  list output_list;
  int new_element;
  //transform every element in the knapsack with normal strong modular multiplication
  for (int i=size(knapsack);i>=1;i--)
  {
    new_element=knapsack[i]*key%mod1;
    output_list = insert(output_list,new_element);
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  list knapsack = 1,3,5,12;
  int key = 3;
  int mod1 = 23;
  merkle_hellman_transformation(knapsack,key,mod1);
}

proc merkle_hellman_encryption(list knapsack, list message)
"USAGE:   merkle_hellman_encryption(knapsack, message)
RETURN:  encrypted integer
EXAMPLE: example merkle_hellman_encryption; shows an example;
"
{
  int solution = 0;
  if (size(knapsack)!=size(message)||(set_is_set(knapsack)==0))
  {
    return(0);
  }
  else
  {
    for (int i=1;i<=size(knapsack);i++)
    {
      solution = solution+knapsack[i]*message[i];
    }
    return(solution);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  //Please note that the values  have been obtained from the example of merkle_hellman_transformation!
  list hardknapsack =3,9,15,13;
  list message = 0,1,0,1;
  merkle_hellman_encryption(hardknapsack,message);
}

proc merkle_hellman_decryption(list knapsack, int key, int mod1, int message)
"USAGE:   merkle_hellman_decryption(knapsack, key, mod1, message)
RETURN:  decrypted binary list
EXAMPLE: example merkle_hellman_decryption; shows an example;
"
{
  int new_element;
  int t = inverse_modulus(key,mod1);
  int transformed_message;
  list binary_list;
  if ((set_is_set(knapsack)==1)&&(key<mod1))
  {
    //reconstruct easy knapsack be multiplying with the inverse modulus t
    list easy_knapsack;
    for (int i=size(knapsack);i>=1;i--)
    {
      new_element=knapsack[i]*t%mod1;
      easy_knapsack = insert(easy_knapsack,new_element);
    }

    //solve the easy knapsack problem with subset_sum01 or subset_sum02
    transformed_message = (message*t)%mod1;
    transformed_message;
    binary_list = subset_sum01(easy_knapsack,transformed_message);
    return(binary_list)
  }
  else
  {
    return(0)
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  //Please note that the values  have been obtained from the example of merkle_hellman_decryption and merkle_hellman_transformation!
  list hardknapsack =3,9,15,13;
  int key = 3;
  int message = 22;
  int mod1 = 23;
  merkle_hellman_decryption(hardknapsack, key, mod1, message);
}

proc super_increasing_knapsack(int ksize)
"USAGE:   super_increasing_knapsack(ksize)
RETURN:  super-increasing knapsack list
EXAMPLE: super_increasing_knapsack; shows an example;
"
{
  list output_list = insert(output_list,1);
  int next_element;

  for (int i=2; i<=ksize; i++)
  {
    next_element = calculate_max_sum(output_list)+1;
    output_list = insert(output_list,next_element);
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  super_increasing_knapsack(10);
}

proc h_increasing_knapsack(int ksize, int h)
"USAGE:   h_increasing_knapsack(ksize, h)
RETURN:  h-increasing knapsack list
EXAMPLE: h_increasing_knapsack; shows an example;
"
{
  int v;
  if (ksize<=h+1)
  {
    return(set_turn(super_increasing_knapsack(ksize)))
  }
  else
  {
    list out = set_turn(super_increasing_knapsack(h+1));
    int next_element;
    for (int i=h+2; i<=ksize; i++)
    {
      next_element = 0;
      for (v=i-h; v<=i-1; v++)
      {
        next_element = next_element+out[v];
      }
      next_element++;
      out = insert(out,next_element,size(out));
    }
    return(out);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  h_increasing_knapsack(10,5);
}

proc injective_knapsack(int ksize, int kmaxelement)
"USAGE:   injective_knapsack(ksize, kmaxelement)
RETURN:  list of injective knapsacks with maximal element kmaxelement and size ksize
EXAMPLE: injective_knapsack; shows an example;
"
{
  //Create a List of size ksize with the greatest possible elements keeping the set structure
  list list_of_lists;
  list A = insert(A,kmaxelement);
  int i;
  for (i=2;i<=ksize;i++)
  {
    A = insert(A,kmaxelement-(i-1));
  }
  A = set_turn(A);
  list_of_lists = insert(list_of_lists,A);

  //Create all possible sets containing the possible elements of A
  int residual;
  int position;
  while(A[1]==kmaxelement)
  {
    residual=3;
    position = ksize;
    while((residual!=0))
    {
      if(A[position]==1)
      {
        A[position]=kmaxelement-position+1;
        residual=1;
        position--;
      }
      else
      {
        A[position]=A[position]-1;
        residual=0;
      }
    }
    //Insert the list into the overall list if its a set
    if (set_is_set(A)==1)
    {
      list_of_lists = insert(list_of_lists,A);
    }
  }
  //delete the first element since it is smaller than kmaxelement
  list_of_lists = delete(list_of_lists,1);
  //delete duplicates
  list_of_lists = set_list_delete_duplicates(list_of_lists);

  //Check if the remaining knapsacks are injective
  list output_list;
  for(i=1;i<=size(list_of_lists);i++)
  {
    if (set_is_injective(list_of_lists[i])==1)
    {
      output_list=insert(output_list,list_of_lists[i]);
    }
  }
  return(output_list);
}
example
{
  "EXAMPLE:"; echo = 2;
  injective_knapsack(3,9);
}

proc calculate_max_sum(list a)
"USAGE:   calculate_max_sum(a)
RETURN:  sum of all elements in a
EXAMPLE: calculate_max_sum; shows an example;
"
{
  int sum = a[1];
  for (int i=2; i<=size(a);i++)
  {
    sum = sum+a[i];
  }
  return(sum);
}
example
{
  "EXAMPLE:"; echo = 2;
  list a = 1,5,3,2,12;
  calculate_max_sum(a);
}

proc set_is_injective(list a)
"USAGE:   set_is_injective(a)
RETURN:  1 if a is injective, 0 otherwise
EXAMPLE: set_is_injective; shows an example;
"
{
  //Create all subsets of the set a
  list subsum = set_subset_set(a);
  list checklist=calculate_max_sum(subsum[1]);
  int calculator;
  for (int i=2; i<=size(subsum);i++)
  {
    //calculate the maximal subset_sum for every subset. Check if there are duplicated subset_sums. If so, a is not injective
    calculator = calculate_max_sum(subsum[i]);
    if (set_contains(checklist, calculator))
    {
      return(0);
    }
    else
    {
      checklist = insert(checklist,calculator);
    }
  }
  return(1);
}
example
{
  "EXAMPLE:"; echo = 2;
  list inj = 1,5,7,41;
  list non_inj = 1,2,3,4;
  set_is_injective(inj);
  set_is_injective(non_inj);
}

proc is_h_injective(list a, int h)
"USAGE:   is_h_injective(a, h)
RETURN:  1 if a is h-injective, 0 otherwise
EXAMPLE: is_h_injective; shows an example;
"
{
  //Create all sets of subsets
  list subsetlist = set_subset_set(a);
  list h_subsetlist;
  //delete every list with elements more than h+1 since they are not needed to check h-injectivity
  for (int i=1; i<=size(subsetlist); i++)
  {
    if(size(subsetlist[i])<=h)
    {
      h_subsetlist = insert(h_subsetlist,subsetlist[i]);
    }
  }

  //Check if the remaining max_sums do not occure more than once
  list checklist=calculate_max_sum(h_subsetlist[1]);
  int calculator;
  for (i=2; i<=size(h_subsetlist);i++)
  {
    calculator = calculate_max_sum(h_subsetlist[i]);
    if (set_contains(checklist, calculator)==1)
    {
      return(0);
    }
    else
    {
      checklist = insert(checklist,calculator);
    }
  }
  return(1);

}
example
{
  "EXAMPLE:"; echo = 2;
  list h_inj = 1,2,4,10,17;
  is_h_injective(h_inj,3);
  //1+2+4+10=17
  is_h_injective(h_inj,4);
}

proc is_fix_injective(list a)
"USAGE:   is_fix_injective(a)
RETURN:  1 if a is fix-injective, 0 otherwise
EXAMPLE: is_fix_injective; shows an example;
"
{
  //Generation of the list-list-list
  list subsetlist = set_subset_set(a);
  list alreadycreatedlist;
  list listoflists;
  list emptylist1;
  list worklist;
  int i;
  int v;
  list checklist;
  int calculator;

  int set_destination;

  //create list of lists wich contain the lists of a certain length as elements
  for (i = 1; i<= size(subsetlist); i++)
  {
    //Determine the size of the acutal list to choose where to insert it in the listoflists
    set_destination = size(subsetlist[i]);
    if (set_contains(alreadycreatedlist,set_destination)==1)
    {
      //There is already an element with the same set size, so just insert it
      listoflists[set_destination] = insert(listoflists[set_destination],subsetlist[i]);
    }
    else
    {
      //There is not yet an element with the same set size, so create a new one
      listoflists[set_destination] = insert(emptylist1,subsetlist[i]);
      alreadycreatedlist = set_insert(alreadycreatedlist,set_destination );
    }
  }

  //Check for injectivity of each seperate list. Works as in injectivity or h-injectivity
  for (v=1; v<=size(listoflists); v++)
  {
    worklist = listoflists[v];

    checklist=calculate_max_sum(worklist[1]);
    for (i=2; i<=size(worklist); i++)
    {
      calculator = calculate_max_sum(worklist[i]);
      if (set_contains(checklist, calculator)==1)
      {
        return(0);
      }
      else
      {
        checklist = insert(checklist,calculator);
      }
    }
  }
  return(1);
}
example
{
  "EXAMPLE:"; echo = 2;
  //this is fix-injective because 17=10+2+4+1 with different numbers of addens.
  list fix_inj = 1,2,4,10,17;
  //this is not fix-injective because 4+1=2+3.
  list not_fix_inj = 1,2,3,4;
  is_fix_injective(fix_inj);
  is_fix_injective(not_fix_inj);
}

proc three_elements(list out, int iterations)
"USAGE:   three_elements(out, iterations)
RETURN:  Injective_knapsack created with the three elements method
EXAMPLE: three_elements; shows an example;
"
{
  int a;
  int b;
  int c;
  int subsum;
  int adden = 1;
  int condition = 0;
  out = set_turn(out);
  if (set_is_injective(out)==0)
  {
    return(0);
  }
  else
  {
    for (int i=1; i<=iterations; i++)
    {
      while(condition==0)
      {
        subsum = calculate_max_sum(out);
        a = 2*subsum+adden;
        b = a+subsum+adden;
        c = b+subsum+1;
        if ((a+b)>(c+subsum))
        {
          condition=1;
        }
        else
        {
          adden++;
        }
      }
      adden =1;
      condition=0;
      out=set_insert(out, a);
      out=set_insert(out, b);
      out=set_insert(out, c);
    }
    return(out);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  //this is fix-injective because 17=10+2+4+1 with different numbers of addens.
  list super_increasing = 1,2,4,10,20;
  list a = three_elements(super_increasing,2);
  a;
  set_is_injective(a);
}
/*
//===============================================================
//=======  Example for DSA  =====================================
//===============================================================
Suppose a file test is given.It contains "Oscar".

//Hash-function MD5 under Linux

md5sum test           8edfe37dae96cfd2466d77d3884d4196

//================================================================

ring R=0,x,dp;

number q=2^19+21;          //524309
number o=2*3*23*number(7883)*number(16170811);

number p=o*q+1;            //9223372036869000547
number b=2;
number g=power(2,o,p);     //8308467587808723131

number a=111111;
number A=power(g,a,p);     //8566038811843553785

number h =decimal("8edfe37dae96cfd2466d77d3884d4196");

                           //189912871665444375716340628395668619670
h= h mod q;                //259847

number k=123456;

number ki=exgcd(k,q)[1];    //50804
                            //inverse von k mod q

number r= power(g,k,p) mod q;    //76646

number s=ki*(h+a*r) mod q;       //2065

//========== signatur is (r,s)=(76646,2065) =====================
//==================== verification  ============================

number si=exgcd(s,q)[1];  //inverse von s mod q
number e1=si*h mod q;
number e2=si*r mod q;
number rr=((power(g,e1,p)*power(A,e2,p)) mod p) mod q;    //76646

//===============================================================
//=======  Example for knapsack  ================================
//===============================================================
ring R=(5^5,t),x,dp;
R;
//   # ground field : 3125
//   primitive element : t
//   minpoly        : 1*t^5+4*t^1+2*t^0
//   number of vars : 1
//        block   1 : ordering dp
//                  : names    x
//        block   2 : ordering C

proc findEx(number n, number g)
{
   int i;
   for(i=0;i<=size(basering)-1;i++)
   {
      if(g^i==n){return(i);}
   }
}

number g=t^3;                       //choice of the primitive root

findEx(t+1,g);
//2091
findEx(t+2,g);
//2291
findEx(t+3,g);
//1043

intvec b=1,2091,2291,1043;          // k=4
int z=199;
intvec v=1043+z,1+z,2091+z,2291+z;  //permutation   pi=(0123)
v;
1242,200,2290,2490

//(1101)=(e_3,e_2,e_1,e_0)
//encoding 2490+2290+1242=6022 und 1+1+0+1=3

//(6022,3) decoding: c-z*c'=6022-199*3=5425

ring S=5,x,dp;
poly F=x5+4x+2;
poly G=reduce((x^3)^5425,std(F));
G;
//x3+x2+x+1

factorize(G);
//[1]:
//   _[1]=1
//   _[2]=x+1
//   _[3]=x-2
//   _[4]=x+2
//[2]:
//   1,1,1,1

//factors x+1,x+2,x+3, i.e.  (1110)=(e_pi(3),e_pi(2),e_pi(1),e_pi(0))

//pi(0)=1,pi(1)=2,pi(2)=3,pi(3)=0 gives:  (1101)

*/