1 /*
2 * program to generate MNT Elliptic curves
3 *
4 * Compile as (for example)
5 * cl /O2 /GX mnt.cpp big.cpp miracl.lib
6 *
7 * See Miyaji, Nakabayashi & Takano "New explicit conditions of elliptic curve
8 * traces for FR-reduction", IEICE Trans. Fundementals., Vol. E84A, No. 5,May
9 * 2001
10 *
11 * Invoke as mnt N where k= 3,4 or 6
12 *
13 * Security level is preportional to the difficulty of the associated
14 * discrete log problem in bits. e.g. if p is 160 bits and k=6, then
15 * discrete log problem is 160*6=960 bits.
16 *
17 * To find actual elliptic curve use CM utility
18 * For example from output:-
19 *
20 * *** Found one - 165 bits
21 * D= 998
22 * p= 28795013727143986241793993458326894706197949151527
23 * NP= 28795013727143986241793989663922216291941998151802
24 * = 78*369166842655692131305051149537464311435153822459
25 * (a 159-bit prime)
26 *
27 * Execute
28 * cm 28795013727143986241793993458326894706197949151527 -D998 -K78
29 *
30 * New! Now also finds "extended" MNT curves. By allowing a co-factor > 1
31 * many more curves will be found.
32 *
33 * However this can cause a problem for our Pell code. A "." indicates that
34 * only a partial search was possible.
35 *
36 * Mike Scott 4/07/2002
37 *
38 */
39
40 #include <cstring>
41 #include <iostream>
42 #include <fstream>
43 #include <iomanip>
44 #include "big.h"
45
46 using namespace std;
47
48 // Solve the Pell equation
49
pell(int max,Big D,Big F,Big & X,Big & Y,Big * SX,Big * SY,BOOL & complete)50 int pell(int max,Big D,Big F,Big& X,Big& Y,Big *SX,Big *SY,BOOL& complete)
51 {
52 int i,j,ns;
53 BOOL SMALLD;
54 Big A,P,Q,SD,G,B,OG,OB,NG,NB,T;
55
56 complete=FALSE;
57 SMALLD=FALSE;
58 if (D<F*F) SMALLD=TRUE;
59
60 SD=sqrt(D);
61
62 if (SD*SD==D) return 0;
63 P=0; Q=1;
64 OG=-P;
65 OB=1;
66 G=Q;
67 B=0;
68 ns=0;
69 X=Y=0;
70
71 for (i=1;;i++)
72 {
73 A=(P+SD)/Q;
74 P=A*Q-P;
75 Q=(D-P*P)/Q;
76 NG=A*G+OG; // G is getting bigger all the time....
77 NB=A*B+OB;
78 OG=G; OB=B;
79 G=NG; B=NB;
80 if (!SMALLD && bits(G)>max) return ns;
81 // abort - these are only solutions
82 T=G*G-D*B*B;
83 if (T == F/4)
84 {
85 SX[ns]=2*G; SY[ns]=2*B;
86 ns++;
87 }
88
89 if (T == F)
90 {
91 SX[ns]=G; SY[ns]=B;
92 ns++;
93 }
94 if (i>0 && Q==1) break;
95 }
96
97 if (i%2==1)
98 for (j=0;j<i;j++)
99 {
100 A=(P+SD)/Q;
101 P=A*Q-P;
102 Q=(D-P*P)/Q;
103 NG=A*G+OG;
104 NB=A*B+OB;
105 OG=G; OB=B;
106 G=NG; B=NB;
107 if (!SMALLD && bits(G)>max) return ns;
108
109 T=G*G-D*B*B;
110 if (T == F/4)
111 {
112 SX[ns]=2*G; SY[ns]=2*B;
113 ns++;
114 }
115 if (T == F)
116 {
117 SX[ns]=G; SY[ns]=B;
118 ns++;
119 }
120 }
121
122 complete=TRUE; // we got to the end....
123 X=G; Y=B; // minimal solution of x^2-dy^2=1
124 // can be used to find more solutions....
125
126 if (SMALLD)
127 { // too small - do it the hard way
128 Big ylim1,ylim2,R;
129 ns=0;
130 if (F<0)
131 {
132 ylim1=sqrt(-F/D);
133 ylim2=sqrt(-F*(X+1)/(2*D));
134 }
135 else
136 {
137 ylim1=0;
138 ylim2=sqrt(F*(X-1)/(2*D));
139 }
140
141 // This might take too long...
142 // Should really implement LMM method here
143
144 if (ylim2-ylim1>300)
145 {
146 cout << "." << flush;
147 ylim2=ylim1+300;
148 }
149 for (B=ylim1;B<ylim2;B+=1)
150 {
151 R=F+D*B*B;
152 if (R<0) continue;
153 G=sqrt(R);
154 if (G*G==R)
155 {
156 SX[ns]=G; SY[ns]=B;
157 ns++;
158 }
159 }
160 }
161
162 return ns;
163 }
164
multiply(Big & td,Big & a,Big & b,Big & c,Big & d)165 void multiply(Big& td,Big& a,Big& b,Big& c,Big& d)
166 { // (c+td.d) = (a+td.b)(c+td.d)
167 Big t;
168 t=a*c+b*d*td;
169 d=c*b+a*d;
170 c=t;
171 }
172
squfree(int d)173 BOOL squfree(int d)
174 { // check if d is square-free
175 int i,s;
176 miracl *mip=get_mip();
177
178 for (i=0;;i++)
179 {
180 s=mip->PRIMES[i];
181 if ((d%(s*s))==0) return FALSE;
182 if ((s*s)>d) break;
183 }
184
185 return TRUE;
186 }
187
results(BOOL fout,ofstream & ofile,int d,Big p,Big nrp,Big ord,BOOL prime_ord)188 void results(BOOL fout,ofstream& ofile,int d,Big p,Big nrp,Big ord,BOOL prime_ord)
189 {
190 Big cf=nrp/ord;
191
192 cout << "*** Found one - " << bits(p) << " bits" << endl;
193 cout << " D= " << d << endl;
194 cout << " p= " << p << " p%24= " << p%24 << endl;
195 cout << "NP= " << nrp << endl;
196 if (cf>1) cout << " = " << cf << "*" << ord << endl;
197 if (prime_ord)
198 {
199 if (cf==1)
200 cout << " (a " << bits(ord) << "-bit prime) !\n" << endl;
201 else
202 cout << " (a " << bits(ord) << "-bit prime)\n" << endl;
203
204 }
205 else cout << endl;
206
207 if (fout)
208 {
209 ofile << "*** Found one - " << bits(p) << " bits" << endl;
210 ofile << " D= " << d << endl;
211 ofile << " p= " << p << " p%24= " << p%24 << endl;
212 ofile << "NP= " << nrp << endl;
213 if (cf>1) ofile << " = " << cf << "*" << ord << endl;
214 if (prime_ord)
215 {
216 if (cf==1)
217 ofile << " (a " << bits(ord) << "-bit prime) !" << endl;
218 else
219 ofile << " (a " << bits(ord) << "-bit prime)" << endl;
220 ofile << " cm " << p << " -D" << d << " -K" << cf << endl << endl;
221 }
222 else ofile << endl;
223 }
224 }
225
main(int argc,char ** argv)226 int main(int argc,char **argv)
227 {
228 ofstream ofile;
229 int MIN_SECURITY,MAX_SECURITY,MIN_D,MAX_D,MIN_P;
230 int ip,j,d,m,mm,solutions,mnt,start,N,precision,max;
231 BOOL fout,complete;
232
233 // Defaults
234 MIN_SECURITY=768;
235 MAX_SECURITY=1536;
236 MIN_D=1;
237 MAX_D=10000;
238 MIN_P=128;
239 fout=FALSE;
240 precision=100;
241 m=1;
242
243 argc--; argv++;
244 if (argc<1)
245 {
246 cout << "Wrong number of parameters" << endl;
247 cout << "Utility to find MNT elliptic curves" << endl;
248 cout << "The output from this program must be processed by the" << endl;
249 cout << "CM (Complex Multiplication) utility to find the actual" << endl;
250 cout << "curve parameters" << endl;
251 cout << "See comments at start of source file mnt.cpp for more details" << endl;
252 cout << "mnt k, where k=3,4 or 6 is the security multiplier" << endl;
253 cout << "Flags allowed, e.g. mnt 6 -c2 -d100 -D1000 -m1000 -M2000" << endl;
254 cout << "-c: cofactor (number of points on curve will be multiple of cofactor)" << endl;
255 cout << "-d: min value of D, -D: max value of D" << endl;
256 cout << "-m: min security level, -M: max security level (RSA equivalent bits)" << endl;
257 cout << "-p: min prime size/field size in bits" << endl;
258 cout << "Increase co-factor for more curves" << endl;
259 cout << "Defaults -c1 -d0 -D10000 -m768 -M1536 -p128" << endl;
260 cout << "To output to a file, use flag -o <filename>" << endl;
261 exit(0);
262 }
263
264 N=atoi(argv[0]);
265 if (N!=3 && N!=4 && N!=6)
266 {
267 cout << "Wrong parameter - should be 3,4 or 6" << endl;
268 exit(0);
269 }
270
271 ip=1;
272 while (ip<argc)
273 { // look for flags
274 if (strncmp(argv[ip],"-d",2)==0)
275 {
276 MIN_D=atoi(argv[ip]+2);
277 ip++;
278 continue;
279 }
280 if (strncmp(argv[ip],"-D",2)==0)
281 {
282 MAX_D=atoi(argv[ip]+2);
283 ip++;
284 continue;
285
286 }
287 if (strncmp(argv[ip],"-m",2)==0)
288 {
289 MIN_SECURITY=atoi(argv[ip]+2);
290 ip++;
291 continue;
292 }
293 if (strncmp(argv[ip],"-M",2)==0)
294 {
295 MAX_SECURITY=atoi(argv[ip]+2);
296 ip++;
297 continue;
298 }
299 if (strncmp(argv[ip],"-p",2)==0)
300 {
301 MIN_P=atoi(argv[ip]+2);
302 ip++;
303 continue;
304 }
305 if (strncmp(argv[ip],"-c",2)==0)
306 {
307 m=atoi(argv[ip]+2);
308 ip++;
309 continue;
310 }
311 if (strcmp(argv[ip],"-o")==0)
312 {
313 ip++;
314 if (ip<argc)
315 {
316 fout=TRUE;
317 ofile.open(argv[ip++]);
318 continue;
319 }
320 else
321 {
322 cout << "Error in command line" << endl;
323 return 0;
324 }
325
326 }
327 cout << "Error in command line" << endl;
328 return 0;
329 }
330
331 if (m<1 || m>8)
332 {
333 cout << "Co-factor must be > 0 and <= 8" << endl;
334 exit(0);
335 }
336
337 miracl *mip=mirsys(precision,0);
338 Big F,T,M,D,W,K,C,mmax,td,X,Y,x,y,w,xn,yn,t0,u0,p,k,nrp,ord,R;
339 Big SX[20],SY[20];
340
341 start=1;
342 if (1<MIN_D) start=MIN_D;
343 mmax=m;
344 mnt=0;
345 for (M=1;M<=mmax;M+=1)
346 {
347 for (R=1;;R+=1)
348 {
349 if (N==4)
350 { // R=1,2,5,... 1 or 2 mod 4
351 if (R%4!=1 && R%4!=2) continue;
352 if (R%3==0) continue;
353 if (R>=4*M) break; // Hasse limit
354 if (gcd(R,M)!=1) continue;
355 T=R*R-4*M*M;
356 }
357 else
358 { // R=1,3,7,9,13,... 1 or 3 mod 6
359 if (R%6!=1 && R%6!=3) continue;
360 if (R%5==0 || R%2==0) continue;
361 if (R>=4*M) break; // Hasse limit
362 if (gcd(R,M)!=1) continue;
363 if (N==3) T=(R+M)*(R+M)-4*M*M;
364 else T=(R-M)*(R-M)-4*M*M;
365 }
366
367 // Check that quadratic polynomial for p does not factor
368 // i.e. b^2-4ac is not a square
369 // The prime p will be of the form M*Phi_N(x)/R + x
370
371 if (T>=0)
372 {
373 W=sqrt(T);
374 if (W*W==T) continue; // p cannot be prime - has factors
375 }
376
377 // There are probably more conditions like these .....
378
379 if (N==6)
380 {
381 if (M==4 && R==3) continue; // p cannot be prime - multiple of 3
382 if (R>3 && (R/3)%3==0) continue;
383 }
384
385 cout << "\nM= " << M << " R= " << R << endl << endl;
386 if (fout) ofile << "\nM= " << M << " R= " << R << endl << endl;
387
388 K=(4*M-R);
389
390 if (N==3) C=2*M+R;
391 if (N==4) C=R;
392 if (N==6) C=-(2*M-R);
393
394 F=C*C-K*K;
395
396 // CM equation is R.D.V^2 = 4*M*Phi_N(x)-R*(x-1)^2
397 // Substitute x=(X-C)/K (to get rid of x term)
398 // Pell equation is then X^2-R*K*D.V^2 = F
399
400 // Our Pell code runs into problems if F is too big
401
402 // This hack makes F a bit smaller - better for Pell
403 // If M is even x must be odd...
404
405 mm=1;
406 if (M%2==0 && F%4==0)
407 {
408 mm=2; F/=4;
409 }
410
411 max=10+MAX_SECURITY/(2*N);
412 gprime(100000); // > 2^16
413
414
415 for (d=start;d<=MAX_D;d++)
416 {
417 // D must be square-free
418 if (!squfree(d)) continue;
419 td=R*K*d;
420 solutions=pell(max,td,F,t0,u0,SX,SY,complete);
421
422 if (!solutions) continue;
423 for (j=0;j<solutions;j++)
424 {
425 X=SX[j];
426 Y=SY[j];
427
428 forever
429 {
430 if (bits(mm*X) > max) break; // primes are too big
431
432 if ((mm*X-C)%K == 0)
433 {
434 x=(mm*X-C)/K;
435 if (N==3) p=x*x+x+1;
436 if (N==4) p=x*x+1;
437 if (N==6) p=x*x-x+1;
438 if (R==1 || p%R==0)
439 {
440 p=(M*p)/R+x;
441 if (N*bits(p) > MAX_SECURITY) break;
442 if (bits(p)>=MIN_P && N*bits(p)>= MIN_SECURITY && prime(p))
443 {
444 nrp=p+1-(x+1);
445 ord=trial_divide(nrp);
446 if (prime(ord))
447 {
448 if (bits(ord)>=MIN_P)
449 {
450 mnt++;
451 results(fout,ofile,d,p,nrp,ord,TRUE);
452 }
453 }
454 // else
455 // { // ord has at least a 16-bit factor...
456 // if (bits(ord)>=MIN_P+16)
457 // {
458 // mnt++;
459 // results(fout,ofile,d,p,nrp,ord,FALSE);
460 // }
461 // }
462 }
463 }
464 }
465 if ((-mm*X-C)%K == 0)
466 {
467 x=(-mm*X-C)/K;
468 if (N==3) p=x*x+x+1;
469 if (N==4) p=x*x+1;
470 if (N==6) p=x*x-x+1;
471 if (R==1 || p%R==0)
472 {
473 p=(M*p)/R+x;
474 if (N*bits(p) > MAX_SECURITY) break;
475 if (bits(p)>=MIN_P && N*bits(p)>= MIN_SECURITY && prime(p))
476 {
477 nrp=p+1-(x+1);
478 ord=trial_divide(nrp);
479 if (prime(ord))
480 {
481 if (bits(ord)>=MIN_P)
482 {
483 results(fout,ofile,d,p,nrp,ord,TRUE);
484 mnt++;
485 }
486 }
487 // else
488 // {
489 // if (bits(ord)>=MIN_P+16)
490 // {
491 // results(fout,ofile,d,p,nrp,ord,FALSE);
492 // mnt++;
493 // }
494 // }
495 }
496 }
497 }
498 if (!complete) break; // no more solutions
499 multiply(td,t0,u0,X,Y);
500 }
501 }
502 }
503 }
504 }
505 if (fout) ofile << mnt << " MNT elliptic curves found in this range" << endl;
506 cout << mnt << " MNT elliptic curves found in this range" << endl;
507 return 0;
508 }
509
510