- Timestamp:
- Jan 8, 2007, 2:14:02 PM (17 years ago)
- Branches:
- (u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b9f50b373314e74e83c7c060a651dd2913e1f033')
- Children:
- 2472afeb87b451c05f3b1487b341a41c01255d7a
- Parents:
- ea3c28a21d6e4e4287e0d97fa6b6dea3f5abf2d9
- File:
-
- 1 edited
-
Singular/LIB/AtkinsTest.lib (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/AtkinsTest.lib
rea3c28a rd1b0065 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: AtkinsTest.lib,v 1. 6 2006-12-14 17:30:46 Singular Exp $";2 version="$Id: AtkinsTest.lib,v 1.7 2007-01-08 13:14:02 pfister Exp $"; 3 3 category="Teaching"; 4 4 info=" … … 14 14 15 15 PROCEDURES: 16 new(L,D) checks if number D already exists in list L17 bubblesort(L) sorts elements (out of Z) of the list L in decreasing order18 disc(N,k) generates a sequence of negative discriminants D with |D|<4N, sort in decreasing order19 Cornacchia(d,p) computes solution (x,y) for the Diophantine equation x^2+d*y^2=p with p prime and 0<d<p20 CornacchiaModified(D,p) computes solution (x,y) for the Diophantine equation x^2+|D|*y^2=4p with p prime21 pFactor1(n,B,P) Pollard's p-factorization22 maximum(L) computes the maximal number contained in list L23 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 w25 e(z,k) computes e^z, i.e. the exponential function of z to the order k26 jot(t,k) computes the j-invariant of the complex number t27 round(r) rounds r to the nearest number out of Z28 HilbertClassPolynomial(D,k) computes the monic polynomial of degree h(D) in Z[X] of which jot((D+sqr(D))/2) is a root29 RootsModp(p,P) computes roots of the polynomial P modulo p with p prime and p>=330 w(D) computes the number of roots of unity in the quadratic order of discriminant D31 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 32 32 "; 33 33 … … 809 809 B describes the number of recursions being calculated 810 810 - 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. 813 814 (1) There exists a prime divisor q of m such that q>(4-th root(N)+1)^2. 814 815 (2) m*P=O(E)=(0:1:0). … … 818 819 " 819 820 { 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) 864 900 { 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; 1365 903 } 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 1368 1229 example 1369 1230 { "EXAMPLE:"; echo = 2;
Note: See TracChangeset
for help on using the changeset viewer.
