///////////////////////////////////////////////////////////////////////////////
version="$Id: atkins.lib,v 1.6 2009-04-06 12:39:02 seelisch Exp $";
category="Teaching";
info="
LIBRARY:  atkins.lib     Procedures for teaching cryptography
AUTHOR:                  Stefan Steidel, steidel@mathematik.uni-kl.de

NOTE: The library contains auxiliary procedures to compute the elliptic
       curve primality test of Atkin and the Atkin's Test itself.
       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:
  newTest(L,D)              checks if number D already exists in list L
  bubblesort(L)             sorts elements of the list L
  disc(N,k)                 generates a list of negative discriminants
  Cornacchia(d,p)           computes solution (x,y) for x^2+d*y^2=p
  CornacchiaModified(D,p)   computes solution (x,y) for x^2+|D|*y^2=4p
  maximum(L)                computes the maximal number contained in L
  sqr(w,k)                  computes the square root of w w.r.t. accuracy k
  expo(z,k)                 computes exp(z)
  jOft(t,k)                 computes the j-invariant of t
  round(r)                  rounds r to the nearest number out of Z
  HilbertClassPoly(D,k)     computes the Hilbert Class Polynomial
  rootsModp(p,P)            computes roots of the polynomial P modulo p
  wUnit(D)                  computes the number of units in Q(sqr(D))
  Atkin(N,K,B)              tries to prove that N is prime

";

LIB "crypto.lib";
LIB "general.lib";
LIB "ntsolve.lib";
LIB "inout.lib";

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

proc newTest(list L, number D)
"USAGE: newTest(L,D);
RETURN:  1, if D does not already exist in L,
        -1, if D does already exist in L
EXAMPLE:example new; shows an example
"
{
  number a=1;                // a=1 bedeutet: D noch nicht in L vorhanden
  int i;
  for(i=1;i<=size(L);i++)
  {
    if(D==L[i])
    {
      a=-1;            // a=-1 bedeutet: D bereits in L vorhanden
      break;
    }
  }
 return(a);
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = 0,x,dp;
    list L=8976,-223456,556,-778,3,-55603,45,766677;
    number D=-55603;
    newTest(L,D);
}

proc bubblesort(list L)
"USAGE: bubblesort(L);
RETURN: list L, sort in decreasing order
EXAMPLE:example bubblesort; shows an example
"
{
  number b;
  int n,i,j;
  while(j==0)
  {
    i=i+1;
    j=1;
    for(n=1;n<=size(L)-i;n++)
    {
      if(L[n]<L[n+1])
      {
        b=L[n];
        L[n]=L[n+1];
        L[n+1]=b;
        j=0;
      }
    }
  }
  return(L);
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = 0,x,dp;
    list L=-567,-233,446,12,-34,8907;
    bubblesort(L);
}

proc disc(number N, int k)
"USAGE: disc(N,k);
RETURN: list L of negative discriminants D, sorted in decreasing order
ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4N
NOTE:   D=b^2-4*a, where 0<=b<=k and intPart((b^2)/4)+1<=a<=k for each b
EXAMPLE:example disc; shows an example
"
{
  list L=-3,-4,-7;
  number D;
  number B;
  int a,b;
  for(b=0;b<=k;b++)
  {
    B=b^2;
    for(a=int(intPart(B/4))+1;a<=k;a++)
    {
      D=-4*a+B;
      if((D<0)&&((D mod 4)!=2)&&((D mod 4)!=3)&&(absValue(D)<4*N)&&(newTest(L,D)==1))
      {
        L[size(L)+1]=D;
      }
    }
  }
  L=bubblesort(L);
  return(L);
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = 0,x,dp;
    disc(2003,50);
}



proc Cornacchia(number d, number p)
"USAGE: Cornacchia(d,p);
RETURN: x,y such that x^2+d*y^2=p with p prime,
        -1, if the Diophantine equation has no solution,
         0, if the parameters are wrong selected
ASSUME: 0<d<p
EXAMPLE:example Cornacchia; shows an example
"
{
  if((d<0)||(p<d))            // (0)[Test if assumptions well-defined]
  {
    return(0);
    // ERROR("Parameters wrong selected! It has to be 0<d<p!");
  }
  else
  {
    number k,x(0),a,b,l,r,c,i;

    k=Jacobi(-d,p);          // (1)[Test if residue]
    if(k==-1)
    {
      return(-1);
      // ERROR("The Diophantine equation has no solution!");
    }
    else
    {
       x(0)=squareRoot(-d,p);   // (2)[Compute square root]
       if((p/2>=x(0))||(p<=x(0)))
       {
           x(0)=-x(0) mod p;
       }
       a=p;
       b=x(0);
       l=intRoot(p);
       while(b>l)       // (3)[Euclidean algorithm]
       {
          r=a mod b;
          a=b;
          b=r;
       }
       c=(p-b^2)/d;     // (4)[Test solution]
       i=intRoot(c);
       if((((p-b^2) mod d)!=0)||(c!=i^2))
       {
          return(-1);
       }
       else
       {
          list L=b,i;
          return(L);
       }
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = 0,x,dp;
    Cornacchia(55,9551);
}

proc CornacchiaModified(number D, number p)
"USAGE: CornacchiaModified(D,p);
RETURN: x,y such that x^2+|D|*y^2=p with p prime,
        -1, if the Diophantine equation has no solution,
         0, if the parameters are wrong selected
ASSUME: D<0, D kongruent 0 or 1 modulo 4 and |D|<4p
EXAMPLE:example CornacchiaModified; shows an example
"
{
  if((D>=0)||((D mod 4)==2)||((D mod 4)==3)||(absValue(D)>=4*p))
  {// (0)[Test if assumptions well-defined]
    return(0);
    // ERROR("Parameters wrong selected!");
  }
  else
  {
    if(p==2)                          // (1)[Case p=2]
    {
      if((D+8)==intRoot(D+8)^2)
      {
        return(intRoot(D+8),1);
      }
      else
      {
        return(-1);
        // ERROR("The Diophantine equation has no solution!");
      }
    }
    else
    {
      number k,x(0),a,b,l,r,c;
      k=Jacobi(D,p);                // (2)[Test if residue]
      if(k==-1)
      {
        return(-1);
        // ERROR("The Diophantine equation has no solution!");
      }
      else
      {
        x(0)=squareRoot(D,p);    // (3)[Compute square root]
        if((x(0) mod 2)!=(D mod 2))
        {
          x(0)=p-x(0);
        }
        a=2*p;
        b=x(0);
        l=intRoot(4*p);
        while(b>l)     // (4)[Euclidean algorithm]
        {
          r=a mod b;
          a=b;
          b=r;
        }
        c=(4*p-b^2)/absValue(D);// (5)[Test solution]
        if((((4*p-b^2) mod absValue(D))!=0)||(c!=intRoot(c)^2))
        {
          return(-1);
          // ERROR("The Diophantine equation has no solution!");
        }
        else
        {
          list L=b,intRoot(c);
          return(L);
        }
      }
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = 0,x,dp;
    CornacchiaModified(-107,1319);
}

proc pFactor1(number n,int B, list P)
"USAGE: pFactor1(n,B,P); n to be factorized, B a bound , P a list of primes
RETURN: a list of factors of n or the message: 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 p-1|k
        then p devides gcd(a^k-1,n) for some random a
EXAMPLE:example pFactor1; shows an example
"
{
  int i;
  number k=1;
  number 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;
  }
  number a=random(2,2147483629);
  number d=gcdN(powerN(a,k,n)-1,n);
  if((d>1)&&(d<n)){return(d);}
  return(n);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,z,dp;
   list L=primList(1000);
   pFactor1(1241143,13,L);
   number h=10;
   h=h^30+25;
   pFactor1(h,20,L);
}

proc maximum(list L)
"USAGE: maximum(list L);
RETURN: the maximal number contained in list L
EXAMPLE:example maximum; shows an example
"
{
  number max=L[1];
  int i;
  for(i=2;i<=size(L);i++)
  {
    if(L[i]>max)
    {
      max=L[i];
    }
  }
  return(max);
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = 0,x,dp;
    list L=465,867,1233,4567,776544,233445,2334,556;
    maximum(L);
}

static proc cmod(number x, number y)
"USAGE: cmod(x,y);
RETURN: x mod y
ASSUME: x,y out of Z and x,y<=2147483647
NOTE:   this algorithm is a helping procedure to be able to calculate
        x mod y with x,y out of Z while working in the complex field
EXAMPLE:example cmod; shows an example
"
{
  int rest=int(x-y*int(x/y));
  if(rest<0)
  {
    rest=rest+int(y);
  }
  return(rest);
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = (complex,30,i),x,dp;
    number x=-1004456;
    number y=1233;
    cmod(x,y);
}

proc sqr(number w, int k)
"USAGE: sqr(w,k);
RETURN: the square root of w
ASSUME: w>=0
NOTE:   k describes the number of decimals being calculated in the real numbers,
        k, intPart(k/5) are inputs for the procedure "nt_solve"
EXAMPLE:example sqr; shows an example
"
{
  poly f=var(1)^2-w;
  def S=basering;
  ring R=(real,k),var(1),dp;
  poly f=imap(S,f);
  ideal I=nt_solve(f,1.1,list(k,int(intPart(k/5))));
  number c=leadcoef(I[1]);
  setring S;
  number c=imap(R,c);
  return(c);
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = (real,60),x,dp;
    number ww=288469650108669535726081;
    sqr(ww,60);
}



proc expo(number z, int k)
"USAGE: expo(z,k);
RETURN: e^z to the order k
NOTE:   k describes the number of summands being calculated in the exponential power series
EXAMPLE:example expo; shows an example
"
{
  number q=1;
  number e=1;
  int n;
  for(n=1;n<=k;n++)
  {
    q=q*z/n;
    e=e+q;
  }
  return(e);
}

example
{ "EXAMPLE:"; echo = 2;
    ring r = (real,30),x,dp;
    number z=40.35;
    expo(z,1000);
}



proc jOft(number t, int k)
"USAGE: jOft(t,k);
RETURN: the j-invariant of t
ASSUME: t is a complex number with positive imaginary part
NOTE:   k describes the number of summands being calculated in the power series,
        10*k is input for the procedure @code{expo}
EXAMPLE:example jOft; shows an example
"
{
  number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j;

  number pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989;

  tr=repart(t);
  ti=impart(t);
  if(tr==-1/2){qr1=-1;}
  if(tr==0){qr1=1;}
  if((tr!=-1/2)&&(tr!=0))
  {
    tr=tr-round(tr);
    qr1=expo(2*i*pi*tr,10*k);
  }

  qi1=expo(-pi*ti,10*k);
  q1=qr1*qi1^2;
  q2=q1^2;

  int n=1;
  while(n<=k)
  {
    m1=m1+(-1)^n*(q1^(n*(3*n-1)/2)+q1^(n*(3*n+1)/2));
    m2=m2+(-1)^n*(q2^(n*(3*n-1)/2)+q2^(n*(3*n+1)/2));
    n++;
  }

  f=q1*((1+m2)/(1+m1))^24;

  j=(256*f+1)^3/f;
  return(j);
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = (complex,30,i),x,dp;
    number t=(-7+i*sqr(7,250))/2;
    jOft(t,50);
}

proc round(number r)
"USAGE: round(r);
RETURN: the nearest number to r out of Z
ASSUME: r should be a rational or a real number
EXAMPLE:example round; shows an example
"
{
  number a=absValue(r);
  number v=r/a;

  number d=10;
  int e;
  while(1)
  {
    e=e+1;
    if(a-d^e<0)
    {
      e=e-1;
      break;
    }
  }

  number b=a;
  int k;
  for(k=0;k<=e;k++)
  {
    while(1)
    {
      b=b-d^(e-k);
      if(b<0)
      {
        b=b+d^(e-k);
        break;
      }
    }
  }

  if(b<1/2)
  {
    return(v*(a-b));
  }
  else
  {
    return(v*(a+1-b));
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = (real,50),x,dp;
    number r=7357683445788723456321.6788643224;
    round(r);
}



proc HilbertClassPoly(number D, int k)
"USAGE: HilbertClassPoly(D,k);
RETURN: the monic polynomial of degree h(D) in Z[X] of which jOft((D+sqr(D))/2) is a root
ASSUME: D is a negative discriminant
NOTE:   k is input for the procedure "jOft",
        5*k is input for the procedure "sqr",
        10*k describes the number of decimals being calculated in the complex numbers
EXAMPLE:example HilbertClassPoly; shows an example
"
{
  if(D>=0)                // (0)[Test if assumptions well-defined]
  {
    ERROR("Parameter wrong selected!");
  }
  else
  {
    def S=basering;
    ring R=0,x,dp;

    string s1,s2,s3;
    number a1,b1,t1,g1;
    number D=imap(S,D);
    number B=intRoot(absValue(D)/3);

    ring C=(complex,10*k,i),x,dp;
    number D=imap(S,D);

    poly P=1;          // (1)[Initialize]
    number b=cmod(D,2);
    number B=imap(R,B);

    number t,a,g,tau,j;
    list L;

    int step=2;
    while(1)
    {
      if(step==2)              // (2)[Initialize a]
      {
        t=(b^2-D)/4;
        L=b,1;
        a=maximum(L);
        step=3;
      }

      if(step==3)               // (3)[Test]
      {
        if((cmod(t,a)!=0))
        {
          step=4;
        }
        else
        {
          s1=string(a);
          s2=string(b);
          s3=string(t);

          setring R;
          execute("a1="+s1+";");
          execute("b1="+s2+";");
          execute("t1="+s3+";");
          g1=gcd(gcd(a1,b1),t1/a1);
          setring C;
          g=imap(R,g1);

          if(g!=1)
          {
            step=4;
          }
          else
          {
            tau=(-b+i*sqr(absValue(D),5*k))/(2*a);
            j=jOft(tau,k);
            if((a==b)||(a^2==t)||(b==0))
            {
              P=P*(var(1)-repart(j));
              step=4;
            }
            else
            {
              P=P*(var(1)^2-2*repart(j)*var(1)+repart(j)^2+impart(j)^2);
              step=4;
            }
          }
        }
      }

      if(step==4)                  // (4)[Loop on a]
      {
        a=a+1;
        if(a^2<=t)
        {
          step=3;
          continue;
        }
        else
        {
          step=5;
        }
      }

      if(step==5)                 // (5)[Loop on b]
      {
        b=b+2;
        if(b<=B)
        {
          step=2;
        }
        else
        {
          break;
        }
      }
    }

    matrix M=coeffs(P,var(1));

    list liste;
    int n;
    for(n=1;n<=nrows(M);n++)
    {
      liste[n]=round(repart(number(M[n,1])));
    }

    poly Q;
    int m;
    for(m=1;m<=size(liste);m++)
    {
      Q=Q+liste[m]*var(1)^(m-1);
    }

    string s=string(Q);
    setring S;
    execute("poly Q="+s+";");
    return(Q);
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = 0,x,dp;
    number D=-23;
    HilbertClassPoly(D,50);
}

proc rootsModp(int p, poly P)
"USAGE: rootsModp(p,P);
RETURN: list of roots of the polynomial P modulo p with p prime
ASSUME: p>=3
NOTE:   this algorithm will be called recursively, and it is understood
        that all the operations are done in Z/pZ (excepting sqareRoot(d,p))
EXAMPLE:example rootsModp; shows an example
"
{
  if(p<3)                         // (0)[Test if assumptions well-defined]
  {
    ERROR("Parameter wrong selected, since p<3!");
  }
  else
  {
    def S=basering;
    ring R=p,var(1),dp;

    poly P=imap(S,P);
    number d;
    int a;
    list L;

    poly A=gcd(var(1)^p-var(1),P);         // (1)[Isolate roots in Z/pZ]
    if(subst(A,var(1),0)==0)
    {
      L[1]=0;
      A=A/var(1);
    }

    if(deg(A)==0)                          // (2)[Small degree?]
    {
      return(L);
    }

    if(deg(A)==1)
    {
      matrix M=coeffs(A,var(1));
      L[size(L)+1]=-leadcoef(M[1,1])/leadcoef(M[2,1]);
      setring S;
      list L=imap(R,L);
      return(L);
    }

    if(deg(A)==2)
    {
      matrix M=coeffs(A,var(1));
      d=leadcoef(M[2,1])^2-4*leadcoef(M[1,1])*leadcoef(M[3,1]);

      ring T=0,var(1),dp;
      number d=imap(R,d);
      number e=squareRoot(d,p);
      setring R;
      number e=imap(T,e);

      L[size(L)+1]=(-leadcoef(M[2,1])+e)/(2*leadcoef(M[3,1]));
      L[size(L)+1]=(-leadcoef(M[2,1])-e)/(2*leadcoef(M[3,1]));
      setring S;
      list L=imap(R,L);
      return(L);
    }

    poly B=1;                       // (3)[Random splitting]
    poly C;
    while((deg(B)==0)||(deg(B)==deg(A)))
    {
      a=random(0,p-1);
      B=gcd((var(1)+a)^((p-1)/2)-1,A);
      C=A/B;
    }

    setring S;                      // (4)[Recurse]
    poly B=imap(R,B);
    poly C=imap(R,C);
    list l=L+rootsModp(p,B)+rootsModp(p,C);
    return(l);
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = 0,x,dp;
    poly f=x4+2x3-5x2+x;
    rootsModp(7,f);
    poly g=x5+112x4+655x3+551x2+1129x+831;
    rootsModp(1223,g);
}

proc wUnit(number D)
"USAGE: wUnit(D);
RETURN: the number of roots of unity in the quadratic order of discriminant D
ASSUME: D<0 a discriminant kongruent to 0 or 1 modulo 4
EXAMPLE:example w; shows an example
"
{
  if((D>=0)||((D mod 4)==2)||((D mod 4)==3))
  {
    ERROR("Parameter wrong selected!");
  }
  else
  {
    if(D<-4) {return(2);}
    if(D==-4){return(4);}
    if(D==-3){return(6);}
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring r = 0,x,dp;
    number D=-3;
    wUnit(D);
}



proc Atkin(number N, int K, int B)
"USAGE: Atkin(N,K,B);
RETURN:  1, if N is prime,
        -1, if N is not prime,
         0, if the algorithm is not applicable, since there are too little discriminants
ASSUME: N is coprime to 6 and different from 1
NOTE:   K/2 is input for the procedure "disc",@*
        K is input for the procedure "HilbertClassPolynomial",@*
        B describes the number of recursions being calculated.@*
        The basis of the algorithm is the following theorem:
          Let N be an integer coprime to 6 and different from 1 and E be an
          ellipic curve modulo N. Assume that we know an integer m and a
          point P of E(Z/NZ) satisfying the following conditions.@*
           (1) There exists a prime divisor q of m such that q>(4-th root(N)+1)^2.@*
           (2) m*P=O(E)=(0:1:0).@*
           (3) (m/q)*P=(x:y:t) with t element of (Z/NZ)*.@*
          Then N is prime.
EXAMPLE:example Atkin; shows an example
"
{
  if(N==1)           {return(-1);} // (0)[Test if assumptions well-defined]
  if((N==2)||(N==3)) {return(1);}
  if(gcdN(N,6)!=1)
  {
    if(printlevel>=1) {"ggT(N,6)="+string(gcdN(N,6));pause();}
    return(-1);
  }
  else
  {
    int i;                         // (1)[Initialize]
    int n(i);
    number N(i)=N;
    if(printlevel>=1) {"Setze i=0, n=0 und N(i)=N(0)="+string(N(i))+".";pause();}

    // declarations:
    int j(0),j(1),j(2),j(3),j(4),k;    // running indices
    list L;                            // all primes smaller than 1000
    list H;                            // sequence of negative discriminants
    number D;                          // discriminant out of H
    list L1,L2,S,S1,S2,R;              // lists of relevant elements
    list P,P1,P2;                      // elliptic points on E(Z/N(i)Z)
    number m,q;                        // m=|E(Z/N(i)Z)| and q|m
    number a,b,j,c;                    // characterize E(Z/N(i)Z)
    number g,u;                        // g out of Z/N(i)Z, u=Jacobi(g,N(i))
    poly T;                            // T=HilbertClassPoly(D,K)
    matrix M;                          // M contains the coefficients of T

    if(printlevel>=1) {"Liste H der moeglichen geeigneten Diskriminanten wird berechnet.";}
    H=disc(N,K/2);
    if(printlevel>=1) {"H="+string(H);pause();}

    int step=2;
    while(1)
    {
      if(step==2)                  // (2)[Is N(i) small??]
      {
        L=5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997;
        for(j(0)=1;j(0)<=size(L);j(0)++)
        {
          if(N(i)==L[j(0)]){return(1);}
          if(((N(i) mod L[j(0)])==0)&&(N(i)!=L[j(0)]))
          {
            if(printlevel>=1) {"N("+string(i)+")="+string(N(i))+" ist durch "+string(L[j(0)])+" teilbar.";pause();}
            step=14;
            break;
          }
        }
        if(step==2)
        {
          step=3;
        }
      }

      if(step==3)                  // (3)[Choose next discriminant]
      {
        n(i)=n(i)+1;
        if(n(i)==size(H)+1)
        {
          if(printlevel>=1) {"Algorithmus nicht anwendbar, da zu wenige geeignete Diskriminanten existieren.";
                             "Erhoehe den Genauigkeitsparameter K und starte den Algorithmus erneut.";pause();}
          return(0);
        }
        D=H[n(i)];
        if(printlevel>=1) {"Naechste Diskriminante D wird gewaehlt. D="+string(D)+".";pause();}
        if(Jacobi(D,N(i))!=1)
        {
          if(printlevel>=1) {"Jacobi(D,N("+string(i)+"))="+string(Jacobi(D,N(i)));pause();}
          continue;
        }
        else
        {
          L1=CornacchiaModified(D,N(i));
          if(size(L1)>1)
          {
            if(printlevel>=1) {"Die Loesung (x,y) der Gleichung x^2+|D|y^2=4N("+string(i)+") lautet";L1;pause();}
            step=4;
          }
          else
          {
            if(L1[1]==-1)
            {
              if(printlevel>=1) {"Die Gleichung x^2+|D|y^2=4N("+string(i)+") hat keine Loesung.";pause();}
              continue;
            }
            if(L1[1]==0)
            {
              if(printLevel>=1) {"Algorithmus fuer N("+string(i)+")="+string(N(i))+" nicht anwendbar,";
                                 "da zu wenige geeignete Diskriminanten existieren.";pause();}
              step=14;
            }
          }
        }
      }

      if(step==4)                  // (4)[Factor m]
      {
        if(printlevel>=1) {"Die Liste L2 der moeglichen m=|E(Z/N("+string(i)+")Z)| wird berechnet.";}
        if(absValue(L1[1])^2<=4*N(i)) {L2=N(i)+1+L1[1],N(i)+1-L1[1];}
        if(D==-4)
        {
          if(absValue(2*L1[2])^2<=4*N(i)) {L2[size(L2)+1]=N(i)+1+2*L1[2];
                                           L2[size(L2)+1]=N(i)+1-2*L1[2];}
        }
// An dieser Stelle wurde "<=4*N(i)" durch "<=16*N(i)" ersetzt.
        if(D==-3)
        {
          if(absValue(L1[1]+3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]+3*L1[2])/2;
                                                  L2[size(L2)+1]=N(i)+1-(L1[1]+3*L1[2])/2;}
          if(absValue(L1[1]-3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]-3*L1[2])/2;
                                                  L2[size(L2)+1]=N(i)+1-(L1[1]-3*L1[2])/2;}
        }
///////////////////////////////////////////////////////////////
        if(size(L2)==0)
        {
          if(printlevel>=1) {"Nach dem Satz von Hasse wurden keine moeglichen m=|E(Z/N("+string(i)+")Z)|";
                             "fuer D="+string(D)+" gefunden.";}
          step=3;
          continue;
        }
        else
        {
          if(printlevel>=1) {"L2=";L2;pause();}
        }

        if(printlevel>=1) {"Die Liste S der Faktoren aller moeglichen m wird berechnet.";}
        S=list();
        for(j(1)=1;j(1)<=size(L2);j(1)++)
        {
          m=L2[j(1)];
          if(m!=0)
          {
            S1=PollardRho(m,10000,1,L);
            S2=pFactor(m,100,L);
            S[size(S)+1]=list(m,S1+S2);
          }
        }
        if(printlevel>=1) {"S=";S;pause();}
        step=5;
      }

      if(step==5)                  // (5)[Does a suitable m exist??]
      {
        for(j(2)=1;j(2)<=size(S);j(2)++)
        {
          m=L2[j(2)];
          for(j(3)=1;j(3)<=size(S[j(2)][2]);j(3)++)
          {
            q=S[j(2)][2][j(3)];
// sqr(sqr(N(i),50),50) ersetzt intRoot(intRoot(N(i)))
            if((q>(sqr(sqr(N(i),50),50)+1)^2) && (MillerRabin(q,5)==1))
            {
              step=6;
              break;
            }
//////////////////////////////////////////////////////
          }
          if(step==6)
          {
            if(printlevel>=1) {"Geeignetes Paar (m,q) gefunden, so dass q|m,";
                               "q>(4-th root(N("+string(i)+"))+1)^2 und q den Miller-Rabin-Test passiert.";
                               "m="+string(m)+",";"q="+string(q);pause();}
            break;
          }
          else
          {
            step=3;
          }
        }
        if(step==3)
        {
          if(printlevel>=1) {"Kein geeignetes Paar (m,q), so dass q|m,";
                             "q>(4-th root(N("+string(i)+"))+1)^2 und q den Miller-Rabin-Test passiert, gefunden.";
                              pause();}
          continue;
        }
      }

      if(step==6)                  // (6)[Compute elliptic curve]
      {
        if(D==-4)
        {
          a=-1;
          b=0;
          if(printlevel>=1) {"Da D=-4, setze a=-1 und b=0.";pause();}
        }
        if(D==-3)
        {
          a=0;
          b=-1;
          if(printlevel>=1) {"Da D=-3, setze a=0 und b=-1.";pause();}
        }
        if(D<-4)
        {
          if(printlevel>=1) {"Das Minimalpolynom T von j((D+sqr(D))/2) aus Z[X] fuer D="+string(D)+" wird berechnet.";}
          T=HilbertClassPoly(D,K);
          if(printlevel>=1) {"T="+string(T);pause();}

          M=coeffs(T,var(1));
          T=0;

          for(j(4)=1;j(4)<=nrows(M);j(4)++)
          {
            M[j(4),1]=leadcoef(M[j(4),1]) mod N(i);
            T=T+M[j(4),1]*var(1)^(j(4)-1);
          }
          if(printlevel>=1) {"Setze T=T mod N("+string(i)+").";"T="+string(T);pause();}

          R=rootsModp(int(N(i)),T);
          if(deg(T)>size(R))
          {
            ERROR("Das Polynom T zerfaellt modulo N("+string(i)+") nicht vollstaendig in Linearfaktoren."
                  "Erhoehe den Genauigkeitsparameter K und starte den Algorithmus erneut.");
          }
          if(printlevel>=1) {if(deg(T)>1) {"Die "+string(deg(T))+" Nullstellen von T modulo N("+string(i)+") sind";
                                           R;pause();}
                             if(deg(T)==1){"Die Nullstelle von T modulo N("+string(i)+") ist";R;pause();}}

          j=R[1];
          c=j*exgcdN(j-1728,N(i))[1];
          a=-3*c mod N(i);
          b=2*c mod N(i);
          if(printlevel>=1) {"Waehle die Nullstelle j="+string(j)+" aus und setze";"c=j/(j-1728) mod N("+string(i)+"), a=-3c mod N("+string(i)+"), b=2c mod N("+string(i)+").";
                             "a="+string(a)+",";"b="+string(b);pause();}
        }
        step=7;
      }

      if(step==7)                  // (7)[Find g]
      {
        if(D==-3)
        {
          while(1)
          {
            g=random(1,2147483647) mod N(i);
            u=Jacobi(g,N(i));
            if((u==-1)&&(powerN(g,(N(i)-1)/3,N(i))!=1))
            {
              if(printlevel>=1) {"g="+string(g);pause();}
              break;
            }
          }
        }
        else
        {
          while(1)
          {
            g=random(1,2147483647) mod N(i);
            u=Jacobi(g,N(i));
            if(u==-1)
            {
              if(printlevel>=1) {"g="+string(g);pause();}
              break;
            }
          }
        }
        step=8;
      }

      if(step==8)                  // (8)[Find P]
      {
        if(printlevel>=1) {"Ein zufaelliger Punkt P auf der Elliptischen Kurve";
                           "mit der Gleichung y^2=x^3+ax+b fuer";"N("+string(i)+")="+string(N(i))+",";
                           "   a="+string(a)+",";"   b="+string(b);"wird gewaehlt.";}
        P=ellipticRandomPoint(N(i),a,b);
        if(printlevel>=1) {"P=("+string(P)+")";pause();}

        if(size(P)==1)
        {
          step=14;
        }
        else
        {
          step=9;
        }
      }

      if(step==9)                  // (9)[Find right curve]
      {
        if(printlevel>=1) {"Die Punkte P2=(m/q)*P und P1=q*P2 auf der Kurve werden berechnet.";}
        P2=ellipticMult(N(i),a,b,P,m/q);
        P1=ellipticMult(N(i),a,b,P2,q);
        if(printlevel>=1) {"P1=("+string(P1)+"),";"P2=("+string(P2)+")";pause();}

        if((P1[1]==0)&&(P1[2]==1)&&(P1[3]==0))
        {
          step=12;
        }
        else
        {
          if(printlevel>=1) {"Da P1!=(0:1:0), ist fuer die Koeffizienten a="+string(a)+" und b="+string(b)+" m!=|E(Z/N("+string(i)+")Z)|.";
                             "Waehle daher neue Koeffizienten a und b.";pause();}
          step=10;
        }
      }

      if(step==10)                 // (10)[Change coefficients]
      {
        k=k+1;
        if(k>=wUnit(D))
        {
          if(printlevel>=1) {"Da k=wUnit(D)="+string(k)+", ist N("+string(i)+")="+string(N(i))+" nicht prim.";pause();}
          step=14;
        }
        else
        {
          if(D<-4) {a=a*g^2 mod N(i); b=b*g^3 mod N(i);
                    if(printlevel>=1) {"Da D<-4, setze a=a*g^2 mod N("+string(i)+") und b=b*g^3 mod N("+string(i)+").";
                                       "a="+string(a)+",";"b="+string(b)+",";"k="+string(k);pause();}}
          if(D==-4){a=a*g mod N(i);
                    if(printlevel>=1) {"Da D=-4, setze a=a*g mod N("+string(i)+").";"a="+string(a)+",";
                                       "b="+string(b)+",";"k="+string(k);pause();}}
          if(D==-3){b=b*g mod N(i);
                    if(printlevel>=1) {"Da D=-3, setze b=b*g mod N("+string(i)+").";"a="+string(a)+",";
                                       "b="+string(b)+",";"k="+string(k);pause();}}
          step=8;
          continue;
        }
      }

      if(step==11)                 // (11)[Find a new P]
      {
        if(printlevel>=1) {"Ein neuer zufaelliger Punkt P auf der Elliptischen Kurve wird gewaehlt,";
                           "da auch P2=(0:1:0).";}
        P=ellipticRandomPoint(N(i),a,b);
        if(printlevel>=1) {"P=("+string(P)+")";pause();}

        if(size(P)==1)
        {
          step=14;
        }
        else
        {
          if(printlevel>=1) {"Die Punkte P2=(m/q)*P und P1=q*P2 auf der Kurve werden berechnet.";}
          P2=ellipticMult(N(i),a,b,P,m/q);
          P1=ellipticMult(N(i),a,b,P2,q);
          if(printlevel>=1) {"P1=("+string(P1)+"),";"P2=("+string(P2)+")";pause();}

          if((P1[1]!=0)||(P1[2]!=1)||(P1[3]!=0))
          {
            if(printlevel>=1) {"Da P1!=(0:1:0), ist, fuer die Koeffizienten a="+string(a)+" und b="+string(b)+", m!=|E(Z/N("+string(i)+")Z)|.";
                               "Waehle daher neue Koeffizienten a und b.";pause();}
            step=10;
            continue;
          }
          else
          {
            step=12;
          }
        }
      }

      if(step==12)                 // (12)[Check P]
      {
        if((P2[1]==0)&&(P2[2]==1)&&(P2[3]==0))
        {
          step=11;
          continue;
        }
        else
        {
          step=13;
        }
      }

      if(step==13)                 // (13)[Recurse]
      {
        if(i<B)
        {
          if(printlevel>=1) {string(i+1)+". Rekursion:";"";
                             "N("+string(i)+")="+string(N(i))+" erfuellt die Bedingungen des zugrunde liegenden Satzes,";
                             "da P1=(0:1:0) und P2[3] aus (Z/N("+string(i)+")Z)*.";"";
                             "Untersuche nun, ob auch der gefundene Faktor q="+string(q)+" diese Bedingungen erfuellt.";
                             "Setze dazu i=i+1, N("+string(i+1)+")=q="+string(q)+" und beginne den Algorithmus von vorne.";pause();}
          i=i+1;
          int n(i);
          number N(i)=q;
          k=0;
          step=2;
          continue;
        }
        else
        {
          if(printlevel>=1) {"N(B)=N("+string(i)+")="+string(N(i))+" erfuellt die Bedingungen des zugrunde liegenden Satzes,";
                             "da P1=(0:1:0) und P2[3] aus (Z/N("+string(i)+")Z)*.";
                             "Insbesondere ist N="+string(N)+" prim.";pause();}
          return(1);
        }
      }

      if(step==14)                 // (14)[Backtrack]
      {
        if(i>0)
        {
          if(printlevel>=1) {"Setze i=i-1 und starte den Algorithmus fuer N("+string(i-1)+")="+string(N(i-1))+" mit";
                             "neuer Diskriminanten von vorne.";pause();}
          i=i-1;
          k=0;
          step=3;
        }
        else
        {
          if(printlevel>=1) {"N(0)=N="+string(N)+" und daher ist N nicht prim.";pause(n);}
          return(-1);
        }
      }
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = 0,x,dp;
    printlevel=1;
    Atkin(7691,100,5);
    Atkin(10000079,100,2);
}