source: git/Singular/LIB/AtkinsTest.lib @ 53e03a6

///////////////////////////////////////////////////////////////////////////////
version="$Id: AtkinsTest.lib,v 1.3 2006-12-11 18:51:12 Singular Exp $";
category="Teaching";
info="
LIBRARY:  AtkinsTest.lib     Procedures for teaching cryptography
AUTHOR:                      Stefan Steidel, Stefan.Steidel@gmx.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:
  new(L,D)                     checks if number D already exists in list L
  bubblesort(L)                sorts elements (out of Z) of the list L in decreasing order
  disc(N,k)                    generates a sequence of negative discriminants D with |D|<4N, sort in decreasing order
  Cornacchia(d,p)              computes solution (x,y) for the Diophantine equation x^2+d*y^2=p with p prime and 0<d<p
  CornacchiaModified(D,p)      computes solution (x,y) for the Diophantine equation x^2+|D|*y^2=4p with p prime
  pFactor1(n,B,P)              Pollard's p-factorization
  maximum(L)                   computes the maximal number contained in list L
  cmod(x,y)                    computes x mod y while working in the complex numbers, e.g. ring C=(complex,30,i),x,dp;
  sqr(w,k)                     computes the square root of w
  e(z,k)                       computes e^z, i.e. the exponential function of z to the order k
  jot(t,k)                     computes the j-invariant of the complex number t
  round(r)                     rounds r to the nearest number out of Z
  HilbertClassPolynomial(D,k)  computes the monic polynomial of degree h(D) in Z[X] of which jot((D+sqr(D))/2) is a root
  RootsModp(p,P)               computes roots of the polynomial P modulo p with p prime and p>=3
  w(D)                         computes the number of roots of unity in the quadratic order of discriminant D
  Atkin(N,K,B)                 tries to prove that N is prime
";

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

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

proc new(list L, number D)
"USAGE: new(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;
    new(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, sort 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)&&(new(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),x(1),a,b,l,r,c,i;
           int j;

           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]
                x(1)=-x(0) mod p;
                while(1)
                   {
                     while((p/2>=x(0))||(p<=x(0)))
                        {
                          x(0)=x(0)+p;
                          if(p<=x(0))
                             {
                               x(0)=-x(0)+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))
                        {
                          if(j==1)
                             {
                               return(-1);
                               // ERROR("The Diophantine equation has no solution!");
                             }

                          else
                             {
                               j=j+1;
                               x(0)=x(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),x(1),a,b,l,r,c,i;
                int j;

                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]
                     x(1)=-x(0) mod p;
                     while(1)
                        {
                          while((0>x(0))||(p<=x(0)))
                             {
                               x(0)=x(0)+p;
                               if(p<x(0))
                                  {
                                    x(0)=-x(0)+p;
                                  }
                             }

                          if((x(0) mod 2)!=D)
                             {
                               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))
                             {
                               if(j==1)
                                  {
                                    return(-1);
                                    // ERROR("The Diophantine equation has no solution!");
                                  }

                               else
                                  {
                                    j=j+1;
                                    x(0)=x(1);
                                  }
                             }

                          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);
}



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 e(number z, int k)
"USAGE: e(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 e; 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;
    e(z,1000);
}



proc jot(number t, int k)
"USAGE: jot(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 "e"
EXAMPLE:example jot; 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=e(2*i*pi*tr,10*k);
         }

      qi1=e(-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=n+1;
         }

      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;
    jot(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 HilbertClassPolynomial(number D, int k)
"USAGE: HilbertClassPolynomial(D,k);
RETURN: the monic polynomial of degree h(D) in Z[X] of which jot((D+sqr(D))/2) is a root
ASSUME: D is a negative discriminant
NOTE:   k is input for the procedure "jot",
        5*k is input for the procedure "sqr",
        10*k describes the number of decimals being calculated in the complex numbers
EXAMPLE:example HilbertClassPolynomial; 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=jot(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;
    HilbertClassPolynomial(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 w(number D)
"USAGE: w(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;
    w(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 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);}
      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=HilbertClassPolynomial(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)
                   {
                     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)++)                         // (2)[Is N(i) small??]
                        {
                          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];}
                        }

                     if(D==-3)
                        {
                          if(absValue(L1[1]+3*L1[2])^2<=4*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<=4*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=pFactor1(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)];
                               if((q>(intRoot(intRoot(N(i)))+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=HilbertClassPolynomial(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)
                   {
                     k=k+1;
                     if(k>=w(D))
                        {
                          if(printlevel>=1) {"Da k=w(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();}
                          return(-1);
                        }
                   }
              }
         }
}
example
{ "EXAMPLE:"; echo = 2;
    ring R = 0,x,dp;
    printlevel=1;
    Atkin(7691,100,5);
    Atkin(8543,100,4);
    Atkin(100019,100,5);
    Atkin(10000079,100,2);
}