///////////////////////////////////////////////////////////////////////////////
version="$Id$";
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 means: D does not already exist in L
  int i;
  for(i=1;i<=size(L);i++)
  {
    if(D==L[i])
    {
      a=-1;                 // a=-1 means: D does already exist in L
      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)!=-1)&&(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+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)==-1)||(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))) // D is <0
        {
          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]
        number root_c=intRoot(c);
        if((((4*p-b^2) mod absValue(D))!=0)||(c!=root_c^2))
        {
          return(-1);
          // ERROR("The Diophantine equation has no solution!");
        }
        else
        {
          list L=b,root_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=gcd(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,k div 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;}
  else
  {
    if(tr==0){qr1=1;}
    else
    {
      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) div 2)+q1^(n*(3*n+1) div 2));
    m2=m2+(-1)^n*(q2^(n*(3*n-1) div 2)+q2^(n*(3*n+1) div 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 squareRoot(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) div 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 few 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 "HilbertClassPoly",@*
        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(gcd(N,6)!=1)
  {
    if(printlevel>=1) {"gcd(N,6) = "+string(gcd(N,6));pause();"";}
    return(-1);
  }
  else
  {
    int i;                         // (1)[Initialize]
    int n(i);
    number N(i)=N;
    if(printlevel>=1) {"Set i = 0, n = 0 and 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) {"List H of possibly suitable discriminants will be calculated.";}
    H=disc(N,K div 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))+" is divisible by "+string(L[j(0)])+".";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) {"Algorithm is not applicable, since there are not enough suitable discriminants.";
                             "Increase the parameter of accuracy K and start the algorithm again.";pause();"";}
          return(0);
        }
        D=H[n(i)];
        if(printlevel>=1) {"Next discriminant D will be chosen. 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) {"The solution (x,y) of the equation x^2+|D|y^2 = 4N("+string(i)+") is";L1;pause();"";}
            step=4;
          }
          else
          {
            if(L1[1]==-1)
            {
              if(printlevel>=1) {"The equation x^2+|D|y^2 = 4N("+string(i)+") has no solution.";pause();"";}
              continue;
            }
            if(L1[1]==0)
            {
              if(printLevel>=1) {"Algorithm is not applicable for N("+string(i)+") = "+string(N(i))+",";
                                 "since there are not enough suitable discriminants.";pause();"";}
              step=14;
            }
          }
        }
      }

      if(step==4)                  // (4)[Factor m]
      {
        if(printlevel>=1) {"List L2 of possible m = |E(Z/N("+string(i)+")Z)| will be calculated.";}
        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];}
        }
// At this point "<=4*N(i)" has been replaced by "<=16*N(i)".
        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) {"Due to the theorem of Hasse there were no possible m = |E(Z/N("+string(i)+")Z)|";
                             "found for D = "+string(D)+".";}
          step=3;
          continue;
        }
        else
        {
          if(printlevel>=1) {"L2 = ";L2;pause();"";}
        }

        if(printlevel>=1) {"List S of factors of all possible m will be calculated.";}
        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) replaces 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) {"Suitable pair (m,q) has been found such that q|m,";
                               "q > (4-th root(N("+string(i)+"))+1)^2 and q passes the Miller-Rabin-Test.";
                               "m = "+string(m)+",";"q = "+string(q);pause();"";}
            break;
          }
          else
          {
            step=3;
          }
        }
        if(step==3)
        {
          if(printlevel>=1) {"No suitable pair (m,q) has been found such that q|m,";
                             "q > (4-th root(N("+string(i)+"))+1)^2 and q passes the Miller-Rabin-Test.";
                              pause();"";}
          continue;
        }
      }

      if(step==6)                  // (6)[Compute elliptic curve]
      {
        if(D==-4)
        {
          a=-1;
          b=0;
          if(printlevel>=1) {"Since D = -4, set a = -1 and b = 0.";pause();"";}
        }
        if(D==-3)
        {
          a=0;
          b=-1;
          if(printlevel>=1) {"Since D = -3, set a = 0 and b = -1.";pause();"";}
        }
        if(D<-4)
        {
          if(printlevel>=1) {"The minimal polynomial T of j((D+sqr(D))/2) in Z[X] will be calculated for D="+string(D)+".";}
          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) {"Set T = T mod N("+string(i)+").";"T = "+string(T);pause();"";}

          R=rootsModp(int(N(i)),T);
          if(deg(T)>size(R))
          {
            ERROR("The polynomial T does not completely split into linear factors modulo N("+string(i)+")."
                  "Increase the parameter of accuracy K and start the algorithm again.");
          }
          if(printlevel>=1) {if(deg(T)>1) {"The "+string(deg(T))+" zeroes of T modulo N("+string(i)+") are";
                                           R;pause();"";}
                             if(deg(T)==1){"The zero of T modulo N("+string(i)+") is";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) {"Choose the zero j = "+string(j)+" and set"; "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) {"A random point P on the elliptic curve corresponding";
                           "to the equation y^2 = x^3+ax+b for";"N("+string(i)+") = "+string(N(i))+",";
                           "   a = "+string(a)+",";"   b = "+string(b);"will be chosen.";}
        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) {"The points P2 = (m/q)*P and P1 = q*P2 on the curve will be calculated.";}
        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) {"Since P1 != (0:1:0), it holds m != |E(Z/N("+string(i)+")Z)| for the coefficients a = "+string(a)+" and b = "+string(b)+".";
                             "Therefore choose new coefficients a and b.";pause();"";}
          step=10;
        }
      }

      if(step==10)                 // (10)[Change coefficients]
      {
        k=k+1;
        if(k>=wUnit(D))
        {
          if(printlevel>=1) {"Since k = wUnit(D) = "+string(k)+", it holds that N("+string(i)+") = "+string(N(i))+" is not prime.";pause();"";}
          step=14;
        }
        else
        {
          if(D<-4) {a=a*g^2 mod N(i); b=b*g^3 mod N(i);
                    if(printlevel>=1) {"Since D < -4, set a = a*g^2 mod N("+string(i)+") and 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) {"Since D = -4, set 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) {"Since D = -3, set 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) {"A new random point P on the elliptic curve will be chosen,";
                           "since also 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) {"The points P2 = (m/q)*P and P1 = q*P2 on the curve will be calculated.";}
          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) {"Since P1 != (0:1:0), it holds m != |E(Z/N("+string(i)+")Z)| for the coefficients a = "+string(a)+" and b = "+string(b)+".";
                               "Therefore choose new coefficients a and 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)+". Recursion:";"";
                             "N("+string(i)+") = "+string(N(i))+" suffices the conditions of the underlying theorem,";
                             "since P1 = (0:1:0) and P2[3] in (Z/N("+string(i)+")Z)*.";"";
                             "Now check if also the found factor q="+string(q)+" suffices these assumptions.";
                             "Therefore set i = i+1, N("+string(i+1)+") = q = "+string(q)+" and restart the algorithm.";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))+" suffices the conditions of the underlying theorem,";
                             "since P1 = (0:1:0) and P2[3] in (Z/N("+string(i)+")Z)*.";
                             "In particular N = "+string(N)+" is prime.";pause();"";}
          return(1);
        }
      }

      if(step==14)                 // (14)[Backtrack]
      {
        if(i>0)
        {
          if(printlevel>=1) {"Set i = i-1 and restart the algorithm for N("+string(i-1)+") = "+string(N(i-1))+" with";
                             "a new discriminant.";pause();"";}
          i=i-1;
          k=0;
          step=3;
        }
        else
        {
          if(printlevel>=1) {"N(0) = N = "+string(N)+" and therefore N is not prime.";pause();"";}
          return(-1);
        }
      }
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = 0,x,dp;
    Atkin(7691,100,5);
    Atkin(3473,10,2);
    printlevel=1;
    Atkin(10000079,100,2);
}