Changeset d1b0065 in git


Ignore:
Timestamp:
Jan 8, 2007, 2:14:02 PM (17 years ago)
Author:
Gerhard Pfister <pfister@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b21a664aa22dc6e196223af8a74ad4885e83547c')
Children:
2472afeb87b451c05f3b1487b341a41c01255d7a
Parents:
ea3c28a21d6e4e4287e0d97fa6b6dea3f5abf2d9
Message:
Hauptprozedur ausgetauscht


git-svn-id: file:///usr/local/Singular/svn/trunk@9640 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/AtkinsTest.lib

    rea3c28a rd1b0065  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: AtkinsTest.lib,v 1.6 2006-12-14 17:30:46 Singular Exp $";
     2version="$Id: AtkinsTest.lib,v 1.7 2007-01-08 13:14:02 pfister Exp $";
    33category="Teaching";
    44info="
     
    1414
    1515PROCEDURES:
    16   new(L,D)        checks if number D already exists in list L
    17   bubblesort(L)   sorts elements (out of Z) of the list L in decreasing order
    18   disc(N,k)       generates a sequence of negative discriminants D with |D|<4N, sort in decreasing order
    19   Cornacchia(d,p) computes solution (x,y) for the Diophantine equation x^2+d*y^2=p with p prime and 0<d<p
    20   CornacchiaModified(D,p)      computes solution (x,y) for the Diophantine equation x^2+|D|*y^2=4p with p prime
    21   pFactor1(n,B,P) Pollard's p-factorization
    22   maximum(L)      computes the maximal number contained in list L
    23   cmod(x,y)       computes x mod y while working in the complex numbers, e.g. ring C=(complex,30,i),x,dp;
    24   sqr(w,k)        computes the square root of w
    25   e(z,k)          computes e^z, i.e. the exponential function of z to the order k
    26   jot(t,k)        computes the j-invariant of the complex number t
    27   round(r)        rounds r to the nearest number out of Z
    28   HilbertClassPolynomial(D,k)  computes the monic polynomial of degree h(D) in Z[X] of which jot((D+sqr(D))/2) is a root
    29   RootsModp(p,P)  computes roots of the polynomial P modulo p with p prime and p>=3
    30   w(D)            computes the number of roots of unity in the quadratic order of discriminant D
    31   Atkin(N,K,B)    tries to prove that N is prime
     16  new(L,D)                     checks if number D already exists in list L
     17  bubblesort(L)                sorts elements of the list L
     18  disc(N,k)                    generates a list of negative discriminants
     19  Cornacchia(d,p)              computes solution (x,y) for x^2+d*y^2=p
     20  CornacchiaModified(D,p)      computes solution (x,y) for x^2+|D|*y^2=4p
     21  maximum(L)                   computes the maximal number contained in L
     22  cmod(x,y)                    computes x mod y
     23  sqr(w,k)                     computes the square root of w
     24  e(z,k)                       computes exp(z)
     25  jot(t,k)                     computes the j-invariant of t
     26  round(r)                     rounds r to the nearest number out of Z
     27  HilbertClassPolynomial(D,k)  computes the Hilbert Class Polynomial
     28  RootsModp(p,P)               computes roots of the polynomial P modulo p
     29  w(D)                         computes the number of units in Q(sqr(D))
     30  Atkin(N,K,B)                 tries to prove that N is prime
     31
    3232";
    3333
     
    809809          B describes the number of recursions being calculated
    810810        - The basis of the the algorithm is the following theorem:
    811           Let N be an integer coprime to 6 and different from 1 and E be an ellipic curve modulo N.
    812           Assume that we know an integer m and a point P of E(Z/NZ) satisfying the following conditions.
     811          Let N be an integer coprime to 6 and different from 1 and E be an
     812          ellipic curve modulo N. Assume that we know an integer m and a
     813          point P of E(Z/NZ) satisfying the following conditions.
    813814           (1) There exists a prime divisor q of m such that q>(4-th root(N)+1)^2.
    814815           (2) m*P=O(E)=(0:1:0).
     
    818819"
    819820{
    820       if(N==1)           {return(-1);}
    821       if((N==2)||(N==3)) {return(1);}
    822       if(gcdN(N,6)!=1)
    823       {
    824         if(printlevel>=1)
    825         {
    826           "ggT(N,6)="+string(gcdN(N,6));
    827           pause();
    828         }
    829         return(-1);
    830       }
    831       else
    832       {
    833          int i;                                                             // (1)[Initialize]
    834          int n(i);
    835          number N(i)=N;
    836          if(printlevel>=1)
    837          {
    838            "Setze i=0, n=0 und N(i)=N(0)="+string(N(i))+".";
    839            pause();
    840          }
    841 
    842          // declarations:
    843          int j(0),j(1),j(2),j(3),j(4),k;                                    // running indices
    844          list L;                                                            // all primes smaller than 1000
    845          list H;                                                            // sequence of negative discriminants
    846          number D;                                                          // discriminant out of H
    847          list L1,L2,S,S1,S2,R;                                              // lists of relevant elements
    848          list P,P1,P2;                                                      // elliptic points on E(Z/N(i)Z)
    849          number m,q;                                                        // m=|E(Z/N(i)Z)| and q|m
    850          number a,b,j,c;                                                    // characterize E(Z/N(i)Z)
    851          number g,u;                                                        // g out of Z/N(i)Z, u=Jacobi(g,N(i))
    852          poly T;                                                            // T=HilbertClassPolynomial(D,K)
    853          matrix M;                                                          // M contains the coefficients of T
    854 
    855          if(printlevel>=1)
    856          {
    857            "Liste H der moeglichen geeigneten Diskriminanten wird berechnet.";
    858          }
    859          H=disc(N,K/2);
    860          if(printlevel>=1) {"H="+string(H);pause();}
    861 
    862          int step=2;
    863          while(1)
     821  if(N==1)           {return(-1);} // (0)[Test if assumptions well-defined]
     822  if((N==2)||(N==3)) {return(1);}
     823  if(gcdN(N,6)!=1)
     824  {
     825    if(printlevel>=1) {"ggT(N,6)="+string(gcdN(N,6));pause();}
     826    return(-1);
     827  }
     828  else
     829  {
     830    int i;                         // (1)[Initialize]
     831    int n(i);
     832    number N(i)=N;
     833    if(printlevel>=1) {"Setze i=0, n=0 und N(i)=N(0)="+string(N(i))+".";pause();}
     834
     835    // declarations:
     836    int j(0),j(1),j(2),j(3),j(4),k;    // running indices
     837    list L;                            // all primes smaller than 1000
     838    list H;                            // sequence of negative discriminants
     839    number D;                          // discriminant out of H
     840    list L1,L2,S,S1,S2,R;              // lists of relevant elements
     841    list P,P1,P2;                      // elliptic points on E(Z/N(i)Z)
     842    number m,q;                        // m=|E(Z/N(i)Z)| and q|m
     843    number a,b,j,c;                    // characterize E(Z/N(i)Z)
     844    number g,u;                        // g out of Z/N(i)Z, u=Jacobi(g,N(i))
     845    poly T;                            // T=HilbertClassPolynomial(D,K)
     846    matrix M;                          // M contains the coefficients of T
     847
     848    if(printlevel>=1) {"Liste H der moeglichen geeigneten Diskriminanten wird berechnet.";}
     849    H=disc(N,K/2);
     850    if(printlevel>=1) {"H="+string(H);pause();}
     851
     852    int step=2;
     853    while(1)
     854    {
     855      if(step==2)                  // (2)[Is N(i) small??]
     856      {
     857        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;
     858        for(j(0)=1;j(0)<=size(L);j(0)++)                         
     859        {
     860          if(((N(i) mod L[j(0)])==0)&&(N(i)!=L[j(0)]))
     861          {
     862            if(printlevel>=1) {"N("+string(i)+")="+string(N(i))+" ist durch "+string(L[j(0)])+" teilbar.";pause();}
     863            step=14;
     864            break;
     865          }
     866        }
     867        if(step==2)
     868        {
     869          step=3;
     870        }
     871      }
     872
     873      if(step==3)                  // (3)[Choose next discriminant]
     874      {
     875        n(i)=n(i)+1;
     876        if(n(i)==size(H)+1)
     877        {
     878          if(printlevel>=1) {"Algorithmus nicht anwendbar, da zu wenige geeignete Diskriminanten existieren.";
     879                             "Erhoehe den Genauigkeitsparameter K und starte den Algorithmus erneut.";pause();}
     880          return(0);
     881        }
     882        D=H[n(i)];
     883        if(printlevel>=1) {"Naechste Diskriminante D wird gewaehlt. D="+string(D)+".";pause();}
     884        if(Jacobi(D,N(i))!=1)
     885        {
     886          if(printlevel>=1) {"Jacobi(D,N("+string(i)+"))="+string(Jacobi(D,N(i)));pause();}
     887          continue;
     888        }
     889        else
     890        {
     891          L1=CornacchiaModified(D,N(i));
     892          if(size(L1)>1)
     893          {
     894            if(printlevel>=1) {"Die Loesung (x,y) der Gleichung x^2+|D|y^2=4N("+string(i)+") lautet";L1;pause();}
     895            step=4;
     896          }
     897          else
     898          {
     899            if(L1[1]==-1)
    864900            {
    865               if(step==2)
    866               {
    867                 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;
    868                 for(j(0)=1;j(0)<=size(L);j(0)++)                         // (2)[Is N(i) small??]
    869                 {
    870                   if(((N(i) mod L[j(0)])==0)&&(N(i)!=L[j(0)]))
    871                   {
    872                     if(printlevel>=1)
    873                     {
    874                       "N("+string(i)+")="+string(N(i))+" ist durch "+string(L[j(0)])+" teilbar.";pause();
    875                     }
    876                     step=14;
    877                     break;
    878                   }
    879                 }
    880 
    881                 if(step==2)
    882                 {
    883                   step=3;
    884                 }
    885               }
    886 
    887               if(step==3)                                                   // (3)[Choose next discriminant]
    888                  {
    889                    n(i)=n(i)+1;
    890                    if(n(i)==size(H)+1)
    891                    {
    892                      if(printlevel>=1)
    893                      {
    894                        "Algorithmus nicht anwendbar, da zu wenige geeignete Diskriminanten existieren.";
    895                        "Erhoehe den Genauigkeitsparameter K und starte den Algorithmus erneut.";pause();
    896                      }
    897                      return(0);
    898                    }
    899 
    900                    D=H[n(i)];
    901                    if(printlevel>=1) {"Naechste Diskriminante D wird gewaehlt. D="+string(D)+".";pause();}
    902 
    903                    if(Jacobi(D,N(i))!=1)
    904                       {
    905                         if(printlevel>=1) {"Jacobi(D,N("+string(i)+"))="+string(Jacobi(D,N(i)));pause();}
    906                         continue;
    907                       }
    908 
    909                    else
    910                    {
    911                      L1=CornacchiaModified(D,N(i));
    912                      if(size(L1)>1)
    913                      {
    914                        if(printlevel>=1)
    915                        {
    916                          "Die Loesung (x,y) der Gleichung x^2+|D|y^2=4N("+string(i)+") lautet";
    917                          L1;
    918                          pause();
    919                        }
    920                        step=4;
    921                      }
    922                      else
    923                      {
    924                        if(L1[1]==-1)
    925                        {
    926                          if(printlevel>=1)
    927                          {
    928                            "Die Gleichung x^2+|D|y^2=4N("+string(i)+") hat keine Loesung.";
    929                            pause();
    930                          }
    931                          continue;
    932                        }
    933 
    934                        if(L1[1]==0)
    935                        {
    936                          if(printLevel>=1)
    937                          {
    938                            "Algorithmus fuer N("+string(i)+")="+string(N(i))+" nicht anwendbar, da zu wenige geeignete Diskriminanten existieren.";
    939                            pause();
    940                          }
    941                          step=14;
    942                        }
    943                      }
    944                    }
    945                  }
    946 
    947               if(step==4)                                                   // (4)[Factor m]
    948                  {
    949                    if(printlevel>=1) {"Die Liste L2 der moeglichen m=|E(Z/N("+string(i)+")Z)| wird berechnet.";}
    950                    if(absValue(L1[1])^2<=4*N(i)) {L2=N(i)+1+L1[1],N(i)+1-L1[1];}
    951                    if(D==-4)
    952                    {
    953                      if(absValue(2*L1[2])^2<=4*N(i))
    954                      {
    955                        L2[size(L2)+1]=N(i)+1+2*L1[2];
    956                        L2[size(L2)+1]=N(i)+1-2*L1[2];
    957                      }
    958                    }
    959 
    960                    if(D==-3)
    961                    {
    962                      if(absValue(L1[1]+3*L1[2])^2<=4*N(i))
    963                      {
    964                        L2[size(L2)+1]=N(i)+1+(L1[1]+3*L1[2])/2;
    965                        L2[size(L2)+1]=N(i)+1-(L1[1]+3*L1[2])/2;
    966                      }
    967                      if(absValue(L1[1]-3*L1[2])^2<=4*N(i))
    968                      {
    969                        L2[size(L2)+1]=N(i)+1+(L1[1]-3*L1[2])/2;
    970                        L2[size(L2)+1]=N(i)+1-(L1[1]-3*L1[2])/2;
    971                      }
    972                    }
    973 
    974                    if(size(L2)==0)
    975                    {
    976                      if(printlevel>=1)
    977                      {
    978                        "Nach dem Satz von Hasse wurden keine moeglichen m=|E(Z/N("+string(i)+")Z)|";
    979                        "fuer D="+string(D)+" gefunden.";
    980                      }
    981                      step=3;
    982                      continue;
    983                    }
    984                    else
    985                    {
    986                      if(printlevel>=1)
    987                      {
    988                        "L2=";L2;
    989                        pause();
    990                      }
    991                    }
    992 
    993                    if(printlevel>=1) {"Die Liste S der Faktoren aller moeglichen m wird berechnet.";}
    994                    S=list();
    995                    for(j(1)=1;j(1)<=size(L2);j(1)++)
    996                    {
    997                      m=L2[j(1)];
    998                      if(m!=0)
    999                      {
    1000                        S1=PollardRho(m,10000,1,L);
    1001                        S2=pFactor1(m,100,L);
    1002                        S[size(S)+1]=list(m,S1+S2);
    1003                      }
    1004                    }
    1005                    if(printlevel>=1)
    1006                    {
    1007                      "S=";S;
    1008                      pause();
    1009                    }
    1010                    step=5;
    1011                  }
    1012 
    1013               if(step==5)                                                   // (5)[Does a suitable m exist??]
    1014               {
    1015                 for(j(2)=1;j(2)<=size(S);j(2)++)
    1016                 {
    1017                   m=L2[j(2)];
    1018                   for(j(3)=1;j(3)<=size(S[j(2)][2]);j(3)++)
    1019                   {
    1020                     q=S[j(2)][2][j(3)];
    1021                     if((q>(intRoot(intRoot(N(i)))+1)^2) && (MillerRabin(q,5)==1))
    1022                     {
    1023                       step=6;
    1024                       break;
    1025                     }
    1026                   }
    1027 
    1028                   if(step==6)
    1029                   {
    1030                     if(printlevel>=1)
    1031                     {
    1032                       "Geeignetes Paar (m,q) gefunden, so dass q|m,";
    1033                       "q>(4-th root(N("+string(i)+"))+1)^2 und q den Miller-Rabin-Test passiert.";
    1034                       "m="+string(m)+",";"q="+string(q);
    1035                       pause();
    1036                     }
    1037                     break;
    1038                   }
    1039                   else
    1040                   {
    1041                     step=3;
    1042                   }
    1043                 }
    1044 
    1045                 if(step==3)
    1046                 {
    1047                   if(printlevel>=1)
    1048                   {
    1049                     "Kein geeignetes Paar (m,q), so dass q|m,";
    1050                     "q>(4-th root(N("+string(i)+"))+1)^2 und q den Miller-Rabin-Test passiert, gefunden.";
    1051                     pause();
    1052                   }
    1053                   continue;
    1054                 }
    1055               }
    1056 
    1057               if(step==6)                                                   // (6)[Compute elliptic curve]
    1058                  {
    1059                    if(D==-4)
    1060                    {
    1061                      a=-1;
    1062                      b=0;
    1063                      if(printlevel>=1)
    1064                      {
    1065                        "Da D=-4, setze a=-1 und b=0.";
    1066                        pause();
    1067                      }
    1068                    }
    1069 
    1070                    if(D==-3)
    1071                    {
    1072                      a=0;
    1073                      b=-1;
    1074                      if(printlevel>=1)
    1075                      {
    1076                        "Da D=-3, setze a=0 und b=-1.";
    1077                        pause();
    1078                      }
    1079                    }
    1080 
    1081                    if(D<-4)
    1082                    {
    1083                      if(printlevel>=1)
    1084                      {
    1085                        "Das Minimalpolynom T von j((D+sqr(D))/2) aus Z[X] fuer D="+string(D)+" wird berechnet.";
    1086                      }
    1087                      T=HilbertClassPolynomial(D,K);
    1088                      if(printlevel>=1)
    1089                      {
    1090                        "T="+string(T);
    1091                        pause();
    1092                      }
    1093 
    1094                      M=coeffs(T,var(1));
    1095                      T=0;
    1096 
    1097                      for(j(4)=1;j(4)<=nrows(M);j(4)++)
    1098                      {
    1099                        M[j(4),1]=leadcoef(M[j(4),1]) mod N(i);
    1100                        T=T+M[j(4),1]*var(1)^(j(4)-1);
    1101                      }
    1102                      if(printlevel>=1)
    1103                      {
    1104                        "Setze T=T mod N("+string(i)+").";"T="+string(T);
    1105                        pause();
    1106                      }
    1107 
    1108                      R=RootsModp(int(N(i)),T);
    1109                      if(deg(T)>size(R))
    1110                      {
    1111                        ERROR("Das Polynom T zerfaellt modulo N("+string(i)+") nicht vollstaendig in Linearfaktoren."+
    1112                              "Erhoehe den Genauigkeitsparameter K und starte den Algorithmus erneut.");
    1113                      }
    1114                      if(printlevel>=1)
    1115                      {
    1116                        if(deg(T)>1)
    1117                        {
    1118                          "Die "+string(deg(T))+" Nullstellen von T modulo N("+string(i)+") sind";
    1119                          R;
    1120                          pause();
    1121                        }
    1122                        if(deg(T)==1)
    1123                        {
    1124                          "Die Nullstelle von T modulo N("+string(i)+") ist";
    1125                          R;
    1126                          pause();
    1127                        }
    1128                      }
    1129 
    1130                      j=R[1];
    1131                      c=j*exgcdN(j-1728,N(i))[1];
    1132                      a=-3*c mod N(i);
    1133                      b=2*c mod N(i);
    1134                      if(printlevel>=1)
    1135                      {
    1136                        "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)+").";
    1137                        "a="+string(a)+",";"b="+string(b);
    1138                        pause();
    1139                      }
    1140                    }
    1141 
    1142                    step=7;
    1143                  }
    1144 
    1145               if(step==7)                                                   // (7)[Find g]
    1146               {
    1147                 if(D==-3)
    1148                 {
    1149                   while(1)
    1150                   {
    1151                     g=random(1,2147483647) mod N(i);
    1152                     u=Jacobi(g,N(i));
    1153                     if((u==-1)&&(powerN(g,(N(i)-1)/3,N(i))!=1))
    1154                     {
    1155                       if(printlevel>=1)
    1156                       {
    1157                         "g="+string(g);
    1158                         pause();
    1159                       }
    1160                       break;
    1161                     }
    1162                   }
    1163                 }
    1164                 else
    1165                 {
    1166                   while(1)
    1167                   {
    1168                     g=random(1,2147483647) mod N(i);
    1169                     u=Jacobi(g,N(i));
    1170                     if(u==-1)
    1171                     {
    1172                       if(printlevel>=1)
    1173                       {
    1174                         "g="+string(g);
    1175                         pause();
    1176                       }
    1177                       break;
    1178                     }
    1179                   }
    1180                 }
    1181 
    1182                 step=8;
    1183               }
    1184 
    1185               if(step==8)                                                   // (8)[Find P]
    1186               {
    1187                 if(printlevel>=1)
    1188                 {
    1189                   "Ein zufaelliger Punkt P auf der Elliptischen Kurve";
    1190                   "mit der Gleichung y^2=x^3+ax+b fuer";"N("+string(i)+")="+string(N(i))+",";"   a="+string(a)+",";"   b="+string(b);"wird gewaehlt.";
    1191                 }
    1192                 P=ellipticRandomPoint(N(i),a,b);
    1193                 if(printlevel>=1) {"P=("+string(P)+")";pause();}
    1194 
    1195                 if(size(P)==1)
    1196                 {
    1197                   step=14;
    1198                 }
    1199                 else
    1200                 {
    1201                   step=9;
    1202                 }
    1203               }
    1204 
    1205               if(step==9)                                                   // (9)[Find right curve]
    1206                  {
    1207                    if(printlevel>=1) {"Die Punkte P2=(m/q)*P und P1=q*P2 auf der Kurve werden berechnet.";}
    1208                    P2=ellipticMult(N(i),a,b,P,m/q);
    1209                    P1=ellipticMult(N(i),a,b,P2,q);
    1210                    if(printlevel>=1) {"P1=("+string(P1)+"),";"P2=("+string(P2)+")";pause();}
    1211 
    1212                    if((P1[1]==0)&&(P1[2]==1)&&(P1[3]==0))
    1213                    {
    1214                      step=12;
    1215                    }
    1216                    else
    1217                    {
    1218                      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)|.";
    1219                                            "Waehle daher neue Koeffizienten a und b.";pause();}
    1220                      step=10;
    1221                    }
    1222                  }
    1223 
    1224               if(step==10)
    1225                  {
    1226                    k=k+1;
    1227                    if(k>=w(D))
    1228                       {
    1229                         if(printlevel>=1) {"Da k=w(D)="+string(k)+", ist N("+string(i)+")="+string(N(i))+" nicht prim.";pause();}
    1230                         step=14;
    1231                       }
    1232 
    1233                    else
    1234                       {
    1235                         if(D<-4) {a=a*g^2 mod N(i); b=b*g^3 mod N(i);
    1236                                   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();}}
    1237                         if(D==-4){a=a*g mod N(i);
    1238                                   if(printlevel>=1) {"Da D=-4, setze a=a*g mod N("+string(i)+").";"a="+string(a)+",";"b="+string(b)+",";"k="+string(k);pause();}}
    1239                         if(D==-3){b=b*g mod N(i);
    1240                                   if(printlevel>=1) {"Da D=-3, setze b=b*g mod N("+string(i)+").";"a="+string(a)+",";"b="+string(b)+",";"k="+string(k);pause();}}
    1241                         step=8;
    1242                         continue;
    1243                       }
    1244                  }
    1245 
    1246               if(step==11)                                                  // (11)[Find a new P]
    1247               {
    1248                 if(printlevel>=1)
    1249                 {
    1250                   "Ein neuer zufaelliger Punkt P auf der Elliptischen Kurve wird gewaehlt,";
    1251                   "da auch P2=(0:1:0).";
    1252                 }
    1253                 P=ellipticRandomPoint(N(i),a,b);
    1254                 if(printlevel>=1)
    1255                 {
    1256                   "P=("+string(P)+")";
    1257                   pause();
    1258                 }
    1259 
    1260                 if(size(P)==1)
    1261                 {
    1262                   step=14;
    1263                 }
    1264                 else
    1265                 {
    1266                   if(printlevel>=1)
    1267                   {
    1268                     "Die Punkte P2=(m/q)*P und P1=q*P2 auf der Kurve werden berechnet.";
    1269                   }
    1270                   P2=ellipticMult(N(i),a,b,P,m/q);
    1271                   P1=ellipticMult(N(i),a,b,P2,q);
    1272                   if(printlevel>=1)
    1273                   {
    1274                     "P1=("+string(P1)+"),";"P2=("+string(P2)+")";
    1275                     pause();
    1276                   }
    1277 
    1278                   if((P1[1]!=0)||(P1[2]!=1)||(P1[3]!=0))
    1279                   {
    1280                     if(printlevel>=1)
    1281                     {
    1282                       "Da P1!=(0:1:0), ist, fuer die Koeffizienten a="+string(a)+" und b="+string(b)+", m!=|E(Z/N("+string(i)+")Z)|.";
    1283                       "Waehle daher neue Koeffizienten a und b.";
    1284                       pause();
    1285                     }
    1286                     step=10;
    1287                     continue;
    1288                   }
    1289                   else
    1290                   {
    1291                     step=12;
    1292                   }
    1293                 }
    1294               }
    1295 
    1296               if(step==12)                                                  // (12)[Check P]
    1297               {
    1298                 if((P2[1]==0)&&(P2[2]==1)&&(P2[3]==0))
    1299                 {
    1300                   step=11;
    1301                   continue;
    1302                 }
    1303                 else
    1304                 {
    1305                   step=13;
    1306                 }
    1307               }
    1308 
    1309               if(step==13)                                                  // (13)[Recurse]
    1310               {
    1311                 if(i<B)
    1312                 {
    1313                   if(printlevel>=1)
    1314                   {
    1315                     string(i+1)+". Rekursion:";"";
    1316                     "N("+string(i)+")="+string(N(i))+" erfuellt die Bedingungen des zugrunde liegenden Satzes,";
    1317                     "da P1=(0:1:0) und P2[3] aus (Z/N("+string(i)+")Z)*.";"";
    1318                     "Untersuche nun, ob auch der gefundene Faktor q="+string(q)+" diese Bedingungen erfuellt.";
    1319                     "Setze dazu i=i+1, N("+string(i+1)+")=q="+string(q)+" und beginne den Algorithmus von vorne.";
    1320                     pause();
    1321                   }
    1322                   i=i+1;
    1323                   int n(i);
    1324                   number N(i)=q;
    1325                   k=0;
    1326                   step=2;
    1327                   continue;
    1328                 }
    1329                 else
    1330                 {
    1331                   if(printlevel>=1)
    1332                   {
    1333                     "N(B)=N("+string(i)+")="+string(N(i))+" erfuellt die Bedingungen des zugrunde liegenden Satzes,";
    1334                     "da P1=(0:1:0) und P2[3] aus (Z/N("+string(i)+")Z)*.";
    1335                     "Insbesondere ist N="+string(N)+" prim.";
    1336                     pause();
    1337                   }
    1338                   return(1);
    1339                 }
    1340               }
    1341 
    1342               if(step==14)                                                  // (14)[Backtrack]
    1343               {
    1344                 if(i>0)
    1345                 {
    1346                   if(printlevel>=1)
    1347                   {
    1348                     "Setze i=i-1 und starte den Algorithmus fuer N("+string(i-1)+")="+string(N(i-1))+" mit neuer Diskriminanten von vorne.";
    1349                     pause();
    1350                   }
    1351                   i=i-1;
    1352                   k=0;
    1353                   step=3;
    1354                 }
    1355                 else
    1356                 {
    1357                   if(printlevel>=1)
    1358                   {
    1359                     "N(0)=N="+string(N)+" und daher ist N nicht prim.";
    1360                     pause();
    1361                   }
    1362                   return(-1);
    1363                 }
    1364               }
     901              if(printlevel>=1) {"Die Gleichung x^2+|D|y^2=4N("+string(i)+") hat keine Loesung.";pause();}
     902              continue;
    1365903            }
    1366          }
    1367 }
     904            if(L1[1]==0)
     905            {
     906              if(printLevel>=1) {"Algorithmus fuer N("+string(i)+")="+string(N(i))+" nicht anwendbar,";
     907                                 "da zu wenige geeignete Diskriminanten existieren.";pause();}
     908              step=14;
     909            }
     910          }
     911        }
     912      }
     913
     914      if(step==4)                  // (4)[Factor m]
     915      {
     916        if(printlevel>=1) {"Die Liste L2 der moeglichen m=|E(Z/N("+string(i)+")Z)| wird berechnet.";}
     917        if(absValue(L1[1])^2<=4*N(i)) {L2=N(i)+1+L1[1],N(i)+1-L1[1];}
     918        if(D==-4)
     919        {
     920          if(absValue(2*L1[2])^2<=4*N(i)) {L2[size(L2)+1]=N(i)+1+2*L1[2];
     921                                           L2[size(L2)+1]=N(i)+1-2*L1[2];}
     922        }
     923// An dieser Stelle wurde "<=4*N(i)" durch "<=16*N(i)" ersetzt.
     924        if(D==-3)
     925        {
     926          if(absValue(L1[1]+3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]+3*L1[2])/2;
     927                                                  L2[size(L2)+1]=N(i)+1-(L1[1]+3*L1[2])/2;}
     928          if(absValue(L1[1]-3*L1[2])^2<=16*N(i)) {L2[size(L2)+1]=N(i)+1+(L1[1]-3*L1[2])/2;
     929                                                  L2[size(L2)+1]=N(i)+1-(L1[1]-3*L1[2])/2;}
     930        }
     931///////////////////////////////////////////////////////////////
     932        if(size(L2)==0)
     933        {
     934          if(printlevel>=1) {"Nach dem Satz von Hasse wurden keine moeglichen m=|E(Z/N("+string(i)+")Z)|";
     935                             "fuer D="+string(D)+" gefunden.";}
     936          step=3;
     937          continue;
     938        }
     939        else
     940        {
     941          if(printlevel>=1) {"L2=";L2;pause();}
     942        }
     943
     944        if(printlevel>=1) {"Die Liste S der Faktoren aller moeglichen m wird berechnet.";}
     945        S=list();
     946        for(j(1)=1;j(1)<=size(L2);j(1)++)
     947        {
     948          m=L2[j(1)];
     949          if(m!=0)
     950          {
     951            S1=PollardRho(m,10000,1,L);
     952            S2=pFactor(m,100,L);
     953            S[size(S)+1]=list(m,S1+S2);
     954          }
     955        }
     956        if(printlevel>=1) {"S=";S;pause();}
     957        step=5;
     958      }
     959
     960      if(step==5)                  // (5)[Does a suitable m exist??]
     961      {
     962        for(j(2)=1;j(2)<=size(S);j(2)++)
     963        {
     964          m=L2[j(2)];
     965          for(j(3)=1;j(3)<=size(S[j(2)][2]);j(3)++)
     966          {
     967            q=S[j(2)][2][j(3)];
     968// sqr(sqr(N(i),50),50) ersetzt intRoot(intRoot(N(i)))
     969            if((q>(sqr(sqr(N(i),50),50)+1)^2) && (MillerRabin(q,5)==1))
     970            {
     971              step=6;
     972              break;
     973            }
     974//////////////////////////////////////////////////////
     975          }
     976          if(step==6)
     977          {
     978            if(printlevel>=1) {"Geeignetes Paar (m,q) gefunden, so dass q|m,";
     979                               "q>(4-th root(N("+string(i)+"))+1)^2 und q den Miller-Rabin-Test passiert.";
     980                               "m="+string(m)+",";"q="+string(q);pause();}
     981            break;
     982          }
     983          else
     984          {
     985            step=3;
     986          }
     987        }
     988        if(step==3)
     989        {
     990          if(printlevel>=1) {"Kein geeignetes Paar (m,q), so dass q|m,";
     991                             "q>(4-th root(N("+string(i)+"))+1)^2 und q den Miller-Rabin-Test passiert, gefunden.";
     992                              pause();}
     993          continue;
     994        }
     995      }
     996
     997      if(step==6)                  // (6)[Compute elliptic curve]
     998      {
     999        if(D==-4)
     1000        {
     1001          a=-1;
     1002          b=0;
     1003          if(printlevel>=1) {"Da D=-4, setze a=-1 und b=0.";pause();}
     1004        }
     1005        if(D==-3)
     1006        {
     1007          a=0;
     1008          b=-1;
     1009          if(printlevel>=1) {"Da D=-3, setze a=0 und b=-1.";pause();}
     1010        }
     1011        if(D<-4)
     1012        {
     1013          if(printlevel>=1) {"Das Minimalpolynom T von j((D+sqr(D))/2) aus Z[X] fuer D="+string(D)+" wird berechnet.";}
     1014          T=HilbertClassPolynomial(D,K);
     1015          if(printlevel>=1) {"T="+string(T);pause();}
     1016
     1017          M=coeffs(T,var(1));
     1018          T=0;
     1019
     1020          for(j(4)=1;j(4)<=nrows(M);j(4)++)
     1021          {
     1022            M[j(4),1]=leadcoef(M[j(4),1]) mod N(i);
     1023            T=T+M[j(4),1]*var(1)^(j(4)-1);
     1024          }
     1025          if(printlevel>=1) {"Setze T=T mod N("+string(i)+").";"T="+string(T);pause();}
     1026
     1027          R=RootsModp(int(N(i)),T);
     1028          if(deg(T)>size(R))
     1029          {
     1030            ERROR("Das Polynom T zerfaellt modulo N("+string(i)+") nicht vollstaendig in Linearfaktoren."
     1031                  "Erhoehe den Genauigkeitsparameter K und starte den Algorithmus erneut.");
     1032          }                                                                       
     1033          if(printlevel>=1) {if(deg(T)>1) {"Die "+string(deg(T))+" Nullstellen von T modulo N("+string(i)+") sind";
     1034                                           R;pause();}
     1035                             if(deg(T)==1){"Die Nullstelle von T modulo N("+string(i)+") ist";R;pause();}}
     1036
     1037          j=R[1];
     1038          c=j*exgcdN(j-1728,N(i))[1];
     1039          a=-3*c mod N(i);
     1040          b=2*c mod N(i);
     1041          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)+").";
     1042                             "a="+string(a)+",";"b="+string(b);pause();}
     1043        }
     1044        step=7;
     1045      }
     1046
     1047      if(step==7)                  // (7)[Find g]
     1048      {
     1049        if(D==-3)
     1050        {
     1051          while(1)
     1052          {
     1053            g=random(1,2147483647) mod N(i);
     1054            u=Jacobi(g,N(i));
     1055            if((u==-1)&&(powerN(g,(N(i)-1)/3,N(i))!=1))
     1056            {
     1057              if(printlevel>=1) {"g="+string(g);pause();}
     1058              break;
     1059            }
     1060          }
     1061        }
     1062        else
     1063        {
     1064          while(1)
     1065          {
     1066            g=random(1,2147483647) mod N(i);
     1067            u=Jacobi(g,N(i));
     1068            if(u==-1)
     1069            {
     1070              if(printlevel>=1) {"g="+string(g);pause();}
     1071              break;
     1072            }
     1073          }
     1074        }
     1075        step=8;
     1076      }
     1077
     1078      if(step==8)                  // (8)[Find P]
     1079      {
     1080        if(printlevel>=1) {"Ein zufaelliger Punkt P auf der Elliptischen Kurve";
     1081                           "mit der Gleichung y^2=x^3+ax+b fuer";"N("+string(i)+")="+string(N(i))+",";
     1082                           "   a="+string(a)+",";"   b="+string(b);"wird gewaehlt.";}
     1083        P=ellipticRandomPoint(N(i),a,b);
     1084        if(printlevel>=1) {"P=("+string(P)+")";pause();}
     1085
     1086        if(size(P)==1)
     1087        {
     1088          step=14;
     1089        }
     1090        else
     1091        {
     1092          step=9;
     1093        }
     1094      }
     1095
     1096      if(step==9)                  // (9)[Find right curve]
     1097      {
     1098        if(printlevel>=1) {"Die Punkte P2=(m/q)*P und P1=q*P2 auf der Kurve werden berechnet.";}
     1099        P2=ellipticMult(N(i),a,b,P,m/q);
     1100        P1=ellipticMult(N(i),a,b,P2,q);
     1101        if(printlevel>=1) {"P1=("+string(P1)+"),";"P2=("+string(P2)+")";pause();}
     1102
     1103        if((P1[1]==0)&&(P1[2]==1)&&(P1[3]==0))
     1104        {
     1105          step=12;
     1106        }
     1107        else
     1108        {
     1109          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)|.";
     1110                             "Waehle daher neue Koeffizienten a und b.";pause();}
     1111          step=10;
     1112        }
     1113      }
     1114
     1115      if(step==10)                 // (10)[Change coefficients]     
     1116      {
     1117        k=k+1;
     1118        if(k>=w(D))
     1119        {
     1120          if(printlevel>=1) {"Da k=w(D)="+string(k)+", ist N("+string(i)+")="+string(N(i))+" nicht prim.";pause();}
     1121          step=14;
     1122        }
     1123        else
     1124        {
     1125          if(D<-4) {a=a*g^2 mod N(i); b=b*g^3 mod N(i);
     1126                    if(printlevel>=1) {"Da D<-4, setze a=a*g^2 mod N("+string(i)+") und b=b*g^3 mod N("+string(i)+").";
     1127                                       "a="+string(a)+",";"b="+string(b)+",";"k="+string(k);pause();}}
     1128          if(D==-4){a=a*g mod N(i);
     1129                    if(printlevel>=1) {"Da D=-4, setze a=a*g mod N("+string(i)+").";"a="+string(a)+",";
     1130                                       "b="+string(b)+",";"k="+string(k);pause();}}
     1131          if(D==-3){b=b*g mod N(i);
     1132                    if(printlevel>=1) {"Da D=-3, setze b=b*g mod N("+string(i)+").";"a="+string(a)+",";
     1133                                       "b="+string(b)+",";"k="+string(k);pause();}}
     1134          step=8;
     1135          continue;
     1136        }
     1137      }
     1138
     1139      if(step==11)                 // (11)[Find a new P]
     1140      {
     1141        if(printlevel>=1) {"Ein neuer zufaelliger Punkt P auf der Elliptischen Kurve wird gewaehlt,";
     1142                           "da auch P2=(0:1:0).";}
     1143        P=ellipticRandomPoint(N(i),a,b);
     1144        if(printlevel>=1) {"P=("+string(P)+")";pause();}
     1145
     1146        if(size(P)==1)
     1147        {
     1148          step=14;
     1149        }
     1150        else
     1151        {
     1152          if(printlevel>=1) {"Die Punkte P2=(m/q)*P und P1=q*P2 auf der Kurve werden berechnet.";}
     1153          P2=ellipticMult(N(i),a,b,P,m/q);
     1154          P1=ellipticMult(N(i),a,b,P2,q);
     1155          if(printlevel>=1) {"P1=("+string(P1)+"),";"P2=("+string(P2)+")";pause();}
     1156
     1157          if((P1[1]!=0)||(P1[2]!=1)||(P1[3]!=0))
     1158          {
     1159            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)|.";
     1160                               "Waehle daher neue Koeffizienten a und b.";pause();}
     1161            step=10;
     1162            continue;
     1163          }
     1164          else
     1165          {
     1166            step=12;
     1167          }
     1168        }
     1169      }
     1170
     1171      if(step==12)                 // (12)[Check P]
     1172      {
     1173        if((P2[1]==0)&&(P2[2]==1)&&(P2[3]==0))
     1174        {
     1175          step=11;
     1176          continue;
     1177        }
     1178        else
     1179        {
     1180          step=13;
     1181        }
     1182      }
     1183
     1184      if(step==13)                 // (13)[Recurse]
     1185      {
     1186        if(i<B)
     1187        {
     1188          if(printlevel>=1) {string(i+1)+". Rekursion:";"";
     1189                             "N("+string(i)+")="+string(N(i))+" erfuellt die Bedingungen des zugrunde liegenden Satzes,";
     1190                             "da P1=(0:1:0) und P2[3] aus (Z/N("+string(i)+")Z)*.";"";
     1191                             "Untersuche nun, ob auch der gefundene Faktor q="+string(q)+" diese Bedingungen erfuellt.";
     1192                             "Setze dazu i=i+1, N("+string(i+1)+")=q="+string(q)+" und beginne den Algorithmus von vorne.";pause();}
     1193          i=i+1;
     1194          int n(i);
     1195          number N(i)=q;
     1196          k=0;
     1197          step=2;
     1198          continue;
     1199        }
     1200        else
     1201        {
     1202          if(printlevel>=1) {"N(B)=N("+string(i)+")="+string(N(i))+" erfuellt die Bedingungen des zugrunde liegenden Satzes,";
     1203                             "da P1=(0:1:0) und P2[3] aus (Z/N("+string(i)+")Z)*.";
     1204                             "Insbesondere ist N="+string(N)+" prim.";pause();}
     1205          return(1);
     1206        }
     1207      }
     1208
     1209      if(step==14)                 // (14)[Backtrack]
     1210      {
     1211        if(i>0)
     1212        {
     1213          if(printlevel>=1) {"Setze i=i-1 und starte den Algorithmus fuer N("+string(i-1)+")="+string(N(i-1))+" mit";
     1214                             "neuer Diskriminanten von vorne.";pause();}
     1215          i=i-1;
     1216          k=0;
     1217          step=3;
     1218        }
     1219        else
     1220        {
     1221          if(printlevel>=1) {"N(0)=N="+string(N)+" und daher ist N nicht prim.";pause();}
     1222          return(-1);
     1223        }
     1224      }
     1225    }
     1226  }
     1227}
     1228
    13681229example
    13691230{ "EXAMPLE:"; echo = 2;
Note: See TracChangeset for help on using the changeset viewer.