]> git.saurik.com Git - apple/security.git/blob - OSX/libsecurity_cryptkit/lib/CurveParamDocs/schoofs.c
projects / apple / security.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Security-59754.80.3.tar.gz
[apple/security.git] / OSX / libsecurity_cryptkit / lib / CurveParamDocs / schoofs.c
1 /* schoofs.c
2
3 Elliptic curve order calculator, for
4
5 y^2 = x^3 + a x + b (mod p)
6
7 (NOTE:
8 This version has order sieving, triggering on the
9 initial b parameter and incrementing same.
10 Parameter details are described in schoof.c)
11
12 Compile with:
13
14 % cc -O schoofs.c giants.c tools.c -lm -o schoofs
15
16 and run, entering the a,b parameters.
17
18 * Change history:
19
20 20 Mar 01 (REC) Added binarysegsquare() and base-4 ladder
21 20 Mar 01 (REC) Bumped MAX_DIGS as in schoof.c
22 4 Feb 99 (REC) Sieving invoked.
23 2 Feb 99 (REC) Added explicit CRT result
24 12 Jan 99 (REC) Removed (hopefully) last of memory crashes
25 20 Jan 98 (REC) Creation
26
27 * c. 1998 Perfectly Scientific, Inc.
28 * All Rights Reserved.
29 *
30 *
31 *************************************************************/
32
33 #include <stdio.h>
34 #include <math.h>
35 #include"giants.h"
36 #include "tools.h"
37
38 #define P_BREAK 32
39
40
41 #define Q 256 /* See schoof.c for explanation. */
42 #define L_MAX 100
43 #define MAX_DIGS (2 + (Q+15)/8) /* Size of (squared) coefficients. */
44 #define MAX_COEFFS ((L_MAX+1)*(L_MAX+2))
45
46 typedef struct
47 {
48 int deg;
49 giant *coe;
50 } polystruct;
51 typedef polystruct *poly;
52
53
54 static giant *mcand, coe, tmp, err, aux, aux2, globx, globy, t1, t2,
55 t3, t4, t5;
56 static poly qscratch, rscratch, sscratch, tscratch, pbuff, pbuff2,
57 precip, cubic, powx, powy, kxn, kyn, kxd, kyd,
58 txn, txd, tyn, tyd, txn1, txd1, tyn1, tyd1,
59 nx, dx, ny, dy, mn, md, tmp1, tmp2, tmp3, tmp4;
60 static poly s[L_MAX+2], smonic;
61 static giant p, a, b, magcheck;
62 int L;
63
64 void quickmodg(giant g, giant x)
65 { int sgn = x->sign;
66
67 if(sgn == 0) return;
68 if(sgn > 0) {
69 if(gcompg(x, g) >= 0) subg(g, x);
70 return;
71 }
72 addg(g,x);
73 return;
74 }
75
76 int
77 log_2(int n) {
78 int c = 1, d = -1;
79 while(c <= n) {
80 c <<= 1;
81 ++d;
82 }
83 return(d);
84 }
85
86 void
87 justifyp(poly x) {
88 int j;
89 for(j = x->deg; j >= 0; j--) {
90 if(!is0(x->coe[j])) break;
91 }
92 x->deg = (j>0)? j : 0;
93 }
94
95 void
96 polyrem(poly x) {
97 int j;
98 for(j=0; j <= x->deg; j++) {
99 modg(p, x->coe[j]);
100 }
101 justifyp(x);
102 }
103
104
105 giant *
106 newa(int n) {
107 giant *p = (giant *)malloc(n*sizeof(giant));
108 int j;
109 for(j=0; j<n; j++) {
110 p[j] = newgiant(MAX_DIGS);
111 }
112 return(p);
113 }
114
115 poly
116 newpoly(int coeffs) {
117 poly pol;
118 pol = (poly) malloc(sizeof(polystruct));
119 pol->coe = (giant *)newa(coeffs);
120 return(pol);
121 }
122
123 void
124 init_all() {
125 int j;
126
127 j = (2*Q + log_2(MAX_COEFFS) + 32 + 15)/16;
128 globx = newgiant(MAX_COEFFS * j);
129 globy = newgiant(MAX_COEFFS * j);
130
131 init_tools(MAX_DIGS);
132 powx = newpoly(MAX_COEFFS);
133 powy = newpoly(MAX_COEFFS);
134 kxn = newpoly(MAX_COEFFS);
135 kxd = newpoly(MAX_COEFFS);
136 kyn = newpoly(MAX_COEFFS);
137 kyd = newpoly(MAX_COEFFS);
138 txn = newpoly(MAX_COEFFS);
139 txd = newpoly(MAX_COEFFS);
140 tyn = newpoly(MAX_COEFFS);
141 tyd = newpoly(MAX_COEFFS);
142 txn1 = newpoly(MAX_COEFFS);
143 txd1 = newpoly(MAX_COEFFS);
144 tyn1 = newpoly(MAX_COEFFS);
145 tyd1 = newpoly(MAX_COEFFS);
146 nx = newpoly(MAX_COEFFS);
147 ny = newpoly(MAX_COEFFS);
148 dx = newpoly(MAX_COEFFS);
149 dy = newpoly(MAX_COEFFS);
150 mn = newpoly(MAX_COEFFS);
151 md = newpoly(MAX_COEFFS);
152 tmp1 = newpoly(MAX_COEFFS);
153 tmp2 = newpoly(MAX_COEFFS);
154 tmp3 = newpoly(MAX_COEFFS);
155 tmp4 = newpoly(MAX_COEFFS);
156 mcand = (giant *)newa(MAX_COEFFS);
157
158 s[0] = newpoly(MAX_COEFFS);
159
160 for(j=1; j<=L_MAX+1; j++) {
161 s[j] = newpoly(j*(j+1));
162 }
163 smonic = newpoly(MAX_COEFFS);
164 /* The next three need extra space for remaindering routine. */
165 qscratch = newpoly(2*MAX_COEFFS);
166 rscratch = newpoly(2*MAX_COEFFS);
167 pbuff = newpoly(2*MAX_COEFFS);
168 pbuff2 = newpoly(MAX_COEFFS);
169 sscratch = newpoly(MAX_COEFFS);
170 tscratch = newpoly(MAX_COEFFS);
171 tmp = newgiant(MAX_DIGS);
172 err = newgiant(MAX_DIGS);
173 coe = newgiant(MAX_DIGS);
174 aux = newgiant(MAX_DIGS);
175 aux2 = newgiant(MAX_DIGS);
176 t1 = newgiant(MAX_DIGS);
177 t2 = newgiant(MAX_DIGS);
178 t3 = newgiant(MAX_DIGS);
179 t4 = newgiant(MAX_DIGS);
180 t5 = newgiant(MAX_DIGS);
181 cubic = newpoly(4);
182 p = newgiant(MAX_DIGS);
183 a = newgiant(MAX_DIGS);
184 b = newgiant(MAX_DIGS);
185 magcheck = newgiant(MAX_DIGS);
186 precip = newpoly(MAX_COEFFS);
187 }
188
189 void
190 atoa(giant *a, giant *b, int n) {
191 int j;
192 for(j=0; j<n; j++) gtog(a[j], b[j]);
193 }
194
195 void
196 ptop(poly x, poly y)
197 /* y := x. */
198 {
199 y->deg = x->deg;
200 atoa(x->coe, y->coe, y->deg+1);
201 }
202
203 void
204 negp(poly y)
205 /* y := -y. */
206 { int j;
207 for(j=0; j <= y->deg; j++) {
208 negg(y->coe[j]);
209 quickmodg(p, y->coe[j]);
210 }
211 }
212
213 int
214 iszer(giant a) {
215
216 if(a->sign == 0) return(1);
217 return(0);
218
219 }
220
221 int
222 iszerop(poly x) {
223 if(x->deg == 0 && iszer(x->coe[0])) return 1;
224 return 0;
225 }
226
227 int
228 is0(giant a) {
229 return(iszer(a));
230 }
231
232 int
233 is1(giant a) {
234 return(isone(a));
235 }
236
237
238 void
239 addp(poly x, poly y)
240 /* y += x. */
241 {
242 int d = x->deg, j;
243
244 if(y->deg > d) d = y->deg;
245 for(j = 0; j <= d; j++) {
246 if((j <= x->deg) && (j <= y->deg)) {
247 addg(x->coe[j], y->coe[j]);
248 quickmodg(p, y->coe[j]);
249 continue;
250 }
251 if((j <= x->deg) && (j > y->deg)) {
252 gtog(x->coe[j], y->coe[j]);
253 quickmodg(p, y->coe[j]);
254 continue;
255 }
256 }
257 y->deg = d;
258 justifyp(y);
259 }
260
261 void
262 subp(poly x, poly y)
263 /* y -= x. */
264 {
265 int d = x->deg, j;
266
267 if(y->deg > d) d = y->deg;
268 for(j = 0; j <= d; j++) {
269 if((j <= x->deg) && (j <= y->deg)) {
270 subg(x->coe[j], y->coe[j]);
271 quickmodg(p, y->coe[j]);
272 continue;
273 }
274 if((j <= x->deg) && (j > y->deg)) {
275 gtog(x->coe[j], y->coe[j]);
276 negg(y->coe[j]);
277 quickmodg(p, y->coe[j]);
278 continue;
279 }
280 }
281 y->deg = d;
282 justifyp(y);
283 }
284
285
286 void
287 grammarmulp(poly a, poly b)
288 /* b *= a. */
289 {
290 int dega = a->deg, degb = b->deg, deg = dega + degb;
291 register int d, kk, first, diffa;
292
293 for(d=deg; d>=0; d--) {
294 diffa = d-dega;
295 itog(0, coe);
296 for(kk=0; kk<=d; kk++) {
297 if((kk>degb)||(kk<diffa)) continue;
298 gtog(b->coe[kk], tmp);
299 mulg(a->coe[d-kk], tmp);
300 modg(p, tmp);
301 addg(tmp, coe);
302 quickmodg(p, coe);
303 }
304 gtog(coe, mcand[d]);
305 }
306 atoa(mcand, b->coe, deg+1);
307 b->deg = deg;
308 justifyp(b);
309 }
310
311 void
312 atop(giant *a, poly x, int deg)
313 /* Copy array to poly, with monic option. */
314 {
315 int adeg = abs(deg);
316 x->deg = adeg;
317 atoa(a, x->coe, adeg);
318 if(deg < 0) {
319 itog(1, x->coe[adeg]);
320 } else {
321 gtog(a[adeg], x->coe[adeg]);
322 }
323 }
324
325 void
326 just(giant g) {
327 while((g->n[g->sign-1] == 0) && (g->sign > 0)) --g->sign;
328 }
329
330 void
331 unstuff_partial(giant g, poly y, int words){
332 int j;
333 for(j=0; j < y->deg; j++) {
334 bcopy((g->n) + j*words, y->coe[j]->n, words*sizeof(short));
335 y->coe[j]->sign = words;
336 just(y->coe[j]);
337 }
338 }
339
340 void
341 stuff(poly x, giant g, int words) {
342 int deg = x->deg, j, coedigs;
343
344 g->sign = words*(1 + deg);
345 for(j=0; j <= deg; j++) {
346 coedigs = (x->coe[j])->sign;
347 bcopy(x->coe[j]->n, (g->n) + j*words, coedigs*sizeof(short));
348 bzero((g->n) + (j*words+coedigs),
349 sizeof(short)*(words-coedigs));
350 }
351 just(g);
352 }
353
354 int maxwords = 0;
355 void
356
357 binarysegmul(poly x, poly y) {
358 int bits = bitlen(p), xwords, ywords, words;
359
360 xwords = (2*bits + log_2(x->deg+1) + 32 + 15)/16;
361 ywords = (2*bits + log_2(y->deg+1) + 32 + 15)/16;
362 if(xwords > ywords) words = xwords; else words = ywords;
363 if(words > maxwords) {
364 maxwords = words;
365 // printf("m: %d\n", words);
366 fflush(stdout);
367 }
368 stuff(x, globx, words);
369 stuff(y, globy, words);
370 mulg(globx, globy);
371 gtog(y->coe[y->deg], globx); /* Save high coeff. */
372 y->deg += x->deg;
373 gtog(globx, y->coe[y->deg]); /* Move high coeff. */
374 unstuff_partial(globy, y, words);
375 mulg(x->coe[x->deg], y->coe[y->deg]); /* resolve high coeff. */
376 justifyp(y);
377 }
378
379 binarysegsquare(poly y) {
380 int bits = bitlen(p), words;
381 words = (2*bits + log_2(y->deg+1) + 32 + 15)/16;
382 stuff(y, globy, words);
383 squareg(globy);
384 gtog(y->coe[y->deg], globx); /* Save high coeff. */
385 y->deg += y->deg;
386 gtog(globx, y->coe[y->deg]); /* Move high coeff. */
387 unstuff_partial(globy, y, words);
388 mulg(y->coe[y->deg], y->coe[y->deg]); /* resolve high coeff. */
389 justifyp(y);
390 }
391
392
393 void
394 assess(poly x, poly y){
395 int max = 0, j;
396 for(j=0; j <= x->deg; j++) {
397 if(bitlen(x->coe[j]) > max) max = bitlen(x->coe[j]);
398 }
399 // printf("max: %d ", max);
400 max = 0;
401 for(j=0; j <= y->deg; j++) {
402 if(bitlen(y->coe[j]) > max) max = bitlen(y->coe[j]);
403 }
404 // printf("%d\n", max);
405
406
407 }
408
409 int
410 pcompp(poly x, poly y) {
411 int j;
412 if(x->deg != y->deg) return 1;
413 for(j=0; j <= x->deg; j++) {
414 if(gcompg(x->coe[j], y->coe[j])) return 1;
415 }
416 return 0;
417 }
418
419 int maxdeg = 0;
420
421 void
422 mulp(poly x, poly y)
423 /* y *= x. */
424 {
425 int n, degx = x->deg, degy = y->deg;
426
427 /*
428 if(degx > max_deg) {
429 max_deg = degx; printf("xdeg: %d\n", degx);
430 }
431
432 if(degy > max_deg) {
433 max_deg = degy; printf("ydeg: %d\n", degy);
434 }
435 */
436 if((degx < P_BREAK) || (degy < P_BREAK)) {
437 grammarmulp(x,y);
438 justifyp(y);
439 return;
440 }
441 if(x==y) binarysegsquare(y);
442 else binarysegmul(x, y);
443 }
444
445 void
446 revp(poly x)
447 /* Reverse the coefficients of x. */
448 { int j, deg = x->deg;
449
450 for(j=0; j <= deg/2; j++) {
451 gtog(x->coe[deg-j], tmp);
452 gtog(x->coe[j], x->coe[deg-j]);
453 gtog(tmp, x->coe[j]);
454 }
455 justifyp(x);
456 }
457
458 void
459 recipp(poly f, int deg)
460 /* f := 1/f, via newton method. */
461 {
462 int lim = deg + 1, prec = 1;
463
464 sscratch->deg = 0; itog(1, sscratch->coe[0]);
465 itog(1, aux);
466 while(prec < lim) {
467 prec <<= 1;
468 if(prec > lim) prec = lim;
469 f->deg = prec-1;
470 ptop(sscratch, tscratch);
471 mulp(f, tscratch);
472 tscratch->deg = prec-1;
473 polyrem(tscratch);
474 subg(aux, tscratch->coe[0]);
475 quickmodg(p, tscratch->coe[0]);
476 mulp(sscratch, tscratch);
477 tscratch->deg = prec-1;
478 polyrem(tscratch);
479 subp(tscratch, sscratch);
480 sscratch->deg = prec-1;
481 polyrem(sscratch);
482 }
483 justifyp(sscratch);
484 ptop(sscratch, f);
485 }
486
487 int
488 left_justifyp(poly x, int start)
489 /* Left-justify the polynomial, checking from position "start." */
490 {
491 int j, shift = 0;
492
493 for(j = start; j <= x->deg; j++) {
494 if(!is0(x->coe[j])) {
495 shift = start;
496 break;
497 }
498 }
499 x->deg -= shift;
500 for(j=0; j<= x->deg; j++) {
501 gtog(x->coe[j+shift], x->coe[j]);
502 }
503 return(shift);
504 }
505
506 void
507 remp(poly x, poly y, int mode)
508 /* y := x (mod y).
509 mode = 0 is normal operation,
510 = 1 jams a fixed reciprocal,
511 = 2 uses the fixed reciprocal.
512 */
513 {
514 int degx = x->deg, degy = y->deg, d, shift;
515
516 if(degy == 0) {
517 y->deg = 0;
518 itog(0, y->coe[0]);
519 return;
520 }
521 d = degx - degy;
522 if(d < 0) {
523 ptop(x, y);
524 return;
525 }
526 revp(x); revp(y);
527 ptop(y, rscratch);
528 switch(mode) {
529 case 0: recipp(rscratch, d);
530 break;
531 case 1: recipp(rscratch, degy); /* degy -1. */
532 ptop(rscratch, precip);
533 rscratch->deg = d; justifyp(rscratch);
534 break;
535 case 2: ptop(precip, rscratch);
536 rscratch->deg = d; justifyp(rscratch);
537 break;
538 }
539 /* Next, a limited-precision multiply. */
540 if(d < degx) { x->deg = d; justifyp(x);}
541 mulp(x, rscratch);
542 rscratch->deg = d;
543 polyrem(rscratch);
544 justifyp(rscratch);
545 x->deg = degx; justifyp(x);
546 mulp(rscratch, y);
547 subp(x, y);
548 negp(y); polyrem(y);
549 shift = left_justifyp(y, d+1);
550 for(d = y->deg+1; d <= degx-shift; d++) itog(0, y->coe[d]);
551 y->deg = degx - shift;
552 revp(y);
553 revp(x);
554 }
555
556 fullmod(poly x) {
557 polyrem(x);
558 ptop(smonic, s[0]);
559 remp(x, s[0], 2);
560 ptop(s[0], x);
561 polyrem(x);
562 }
563
564 scalarmul(giant s, poly x) {
565 int j;
566 for(j=0; j <= x->deg; j++) {
567 mulg(s, x->coe[j]);
568 modg(p, x->coe[j]);
569 }
570 }
571
572
573 schain(int el) {
574 int j;
575
576 s[0]->deg = 0;
577 itog(0, s[0]->coe[0]);
578
579 s[1]->deg = 0;
580 itog(1, s[1]->coe[0]);
581 s[2]->deg = 0;
582 itog(2, s[2]->coe[0]);
583
584 s[3]->deg = 4;
585 gtog(a, aux); mulg(a, aux); negg(aux);
586 gtog(aux, s[3]->coe[0]);
587 gtog(b, aux); smulg(12, aux);
588 gtog(aux, s[3]->coe[1]);
589 gtog(a, aux); smulg(6, aux);
590 gtog(aux, s[3]->coe[2]);
591 itog(0, s[3]->coe[3]);
592 itog(3, s[3]->coe[4]);
593
594 s[4]->deg = 6;
595 gtog(a, aux); mulg(a, aux); mulg(a, aux);
596 gtog(b, tmp); mulg(b, tmp); smulg(8, tmp); addg(tmp, aux);
597 negg(aux);
598 gtog(aux, s[4]->coe[0]);
599 gtog(b, aux); mulg(a, aux); smulg(4, aux); negg(aux);
600 gtog(aux, s[4]->coe[1]);
601 gtog(a, aux); mulg(a, aux); smulg(5, aux); negg(aux);
602 gtog(aux, s[4]->coe[2]);
603 gtog(b, aux); smulg(20, aux);
604 gtog(aux, s[4]->coe[3]);
605 gtog(a, aux); smulg(5, aux);
606 gtog(aux, s[4]->coe[4]);
607 itog(0, s[4]->coe[5]);
608 itog(1, s[4]->coe[6]);
609 itog(4, aux);
610 scalarmul(aux, s[4]);
611 cubic->deg = 3;
612 itog(1, cubic->coe[3]);
613 itog(0, cubic->coe[2]);
614 gtog(a, cubic->coe[1]);
615 gtog(b, cubic->coe[0]);
616 for(j=5; j <= el; j++) {
617 if(j % 2 == 0) {
618 ptop(s[j/2 + 2], s[j]); mulp(s[j/2-1], s[j]);
619 polyrem(s[j]); mulp(s[j/2-1], s[j]); polyrem(s[j]);
620 ptop(s[j/2-2], s[0]); mulp(s[j/2+1], s[0]); polyrem(s[0]);
621 mulp(s[j/2+1], s[0]); polyrem(s[0]);
622 subp(s[0], s[j]); mulp(s[j/2], s[j]); polyrem(s[j]);
623 gtog(p, aux); iaddg(1, aux); gshiftright(1, aux);
624 scalarmul(aux, s[j]);
625 } else {
626 ptop(s[(j-1)/2+2], s[j]);
627 mulp(s[(j-1)/2], s[j]); polyrem(s[j]);
628 mulp(s[(j-1)/2], s[j]); polyrem(s[j]);
629 mulp(s[(j-1)/2], s[j]); polyrem(s[j]);
630 ptop(s[(j-1)/2-1], s[0]);
631 mulp(s[(j-1)/2 + 1], s[0]); polyrem(s[0]);
632 mulp(s[(j-1)/2 + 1], s[0]); polyrem(s[0]);
633 mulp(s[(j-1)/2 + 1], s[0]); polyrem(s[0]);
634 if(((j-1)/2) % 2 == 1) {
635 mulp(cubic, s[0]); polyrem(s[0]);
636 mulp(cubic, s[0]); polyrem(s[0]);
637 } else {
638 mulp(cubic, s[j]); polyrem(s[j]);
639 mulp(cubic, s[j]); polyrem(s[j]);
640 }
641 // polyout(s[1]); polyout(s[3]); polyout(s[0]); polyout(s[j]);
642 subp(s[0], s[j]);
643 polyrem(s[j]);
644 }
645 }
646 }
647
648 init_recip(int el) {
649 int j;
650 ptop(s[el], smonic);
651 if(el == 2) {
652 mulp(cubic, smonic); polyrem(smonic);
653 }
654 gtog(smonic->coe[smonic->deg], aux); /* High coeff. */
655 binvg(p, aux);
656 scalarmul(aux, smonic); /* s is now monic. */
657 s[0]->deg = smonic->deg + 1;
658 for(j=0; j <= s[0]->deg; j++) itog(1, s[0]->coe[j]);
659 ptop(smonic, pbuff);
660 remp(s[0], pbuff, 1); /* Initialize reciprocal of s as precip. */
661 }
662
663 void powerpoly(poly x, giant n)
664 /* Base-4 window. */
665 { int pos, code;
666 ptop(x, pbuff); /* x. */
667 ptop(pbuff, pbuff2);
668 mulmod(pbuff2, pbuff2); mulmod(pbuff, pbuff2); /* x^3. */
669 pos = bitlen(n)-2;
670 while(pos >= 0) {
671 mulmod(x, x);
672 if(pos==0) {
673 if(bitval(n, pos) != 0) {
674 mulmod(pbuff, x);
675 }
676 break;
677 }
678 code = (bitval(n, pos) != 0) * 2 + (bitval(n, pos-1) != 0);
679 switch(code) {
680 case 0: mulmod(x,x); break;
681 case 1: mulmod(x,x);
682 mulmod(pbuff, x);
683 break;
684 case 2: mulmod(pbuff, x);
685 mulmod(x,x); break;
686 case 3: mulmod(x,x); mulmod(pbuff2, x); break;
687 }
688 pos -= 2;
689 }
690 }
691
692 mulmod(poly x, poly y) {
693 mulp(x, y); fullmod(y);
694 }
695
696 elldoublep(poly n1, poly d1, poly m1, poly c1, poly n0, poly d0,
697 poly m0, poly c0) {
698
699 ptop(n1, mn); mulmod(n1, mn);
700 ptop(mn, pbuff); addp(mn, mn); addp(pbuff, mn);
701 fullmod(mn);
702 ptop(d1, pbuff); mulmod(d1, pbuff);
703 scalarmul(a, pbuff); addp(pbuff, mn);
704 fullmod(mn);
705 mulmod(c1, mn);
706 ptop(m1, md); addp(md, md);
707 mulmod(d1, md); mulmod(d1, md); mulmod(cubic, md);
708
709 ptop(d1, n0); mulmod(mn, n0); mulmod(mn, n0);
710 mulmod(cubic, n0);
711 ptop(n1, pbuff); addp(pbuff, pbuff); fullmod(pbuff);
712 mulmod(md, pbuff); mulmod(md, pbuff);
713 subp(pbuff, n0); fullmod(n0);
714 ptop(md, d0); mulmod(md, d0); mulmod(d1, d0);
715
716 ptop(mn, m0); mulmod(c1, m0);
717 ptop(d0, pbuff); mulmod(n1, pbuff);
718 ptop(n0, c0); mulmod(d1, c0); subp(c0, pbuff);
719 fullmod(pbuff);
720 mulmod(pbuff, m0);
721 ptop(m1, pbuff); mulmod(md, pbuff);
722 mulmod(d1, pbuff); mulmod(d0, pbuff);
723 subp(pbuff, m0); fullmod(m0);
724
725 ptop(c1, c0); mulmod(md, c0); mulmod(d1, c0); mulmod(d0, c0);
726 }
727
728 elladdp(poly n1, poly d1, poly m1, poly c1, poly n2, poly d2, poly m2, poly c2, poly n0, poly d0, poly m0, poly c0) {
729 ptop(m2, mn); mulmod(c1, mn);
730 ptop(m1, pbuff); mulmod(c2, pbuff);
731 subp(pbuff, mn); fullmod(mn);
732 mulmod(d1, mn); mulmod(d2, mn);
733
734 ptop(n2, md); mulmod(d1, md);
735 ptop(n1, pbuff); mulmod(d2, pbuff);
736 subp(pbuff, md); fullmod(md);
737 mulmod(c1, md); mulmod(c2, md);
738
739 ptop(cubic, n0); mulmod(mn, n0); mulmod(mn, n0);
740 mulmod(d1, n0); mulmod(d2, n0);
741 ptop(n1, pbuff); mulmod(d2, pbuff);
742 ptop(n2, d0); mulmod(d1, d0);
743 addp(d0, pbuff); mulmod(md, pbuff); mulmod(md, pbuff);
744 subp(pbuff, n0); fullmod(n0);
745
746 ptop(md, d0); mulmod(md, d0); mulmod(d1, d0); mulmod(d2, d0);
747
748 ptop(mn, m0); mulmod(c1, m0);
749 ptop(d0, pbuff); mulmod(n1, pbuff);
750 ptop(d1, c0); mulmod(n0, c0);
751 subp(c0, pbuff); fullmod(pbuff);
752 mulmod(pbuff, m0);
753 ptop(m1, pbuff); mulmod(md, pbuff);
754 mulmod(d0, pbuff); mulmod(d1, pbuff);
755 subp(pbuff, m0); fullmod(m0);
756
757 ptop(md, c0); mulmod(c1, c0); mulmod(d0, c0); mulmod(d1, c0);
758
759 }
760
761 polyout(poly x) {
762 int j;
763 for(j=0; j <= x->deg; j++) {printf("%d: ",j); gout(x->coe[j]);}
764 }
765
766 #define SIEVE_CUT 7
767
768 main(int argc, char **argv) {
769 int j, ct, el, xmatch, ymatch;
770 int k, k2, t, linit;
771 int T[L_MAX], P[L_MAX];
772 giant ss[L_MAX];
773 unsigned int ord, ordminus;
774
775 printf("Initializing...\n"); fflush(stdout);
776 init_all();
777
778 for(j=0; j < L_MAX; j++) { /* Array to be used for CRT reconstruction. */
779 ss[j] = NULL;
780 }
781
782 printf("Give p, a, b on separate lines:\n"); fflush(stdout);
783 gin(p); /* Field prime. */
784 if(bitlen(p) > Q) {
785 fprintf(stderr, "p too large, larger than %d bits.\n", Q);
786 exit(0);
787 }
788
789 gin(a); gin(b); /* Curve parameters. */
790
791 newb:
792 printf("Checking discriminant for:\nb = ");
793 gout(b);
794 gtog(a, t1); squareg(t1); modg(p, t1); mulg(a, t1); modg(p, t1);
795 gshiftleft(2, t1); /* 4 a^3. */
796 gtog(b, t2); squareg(t2); modg(p, t2); smulg(27, t2);
797 addg(t2, t1); modg(p, t1);
798 if(isZero(t1)) {
799 printf("Discriminant FAILED\n");
800 iaddg(1, b);
801 goto newb;
802 }
803 printf("Discriminant PASSED\n");
804 printf("Starting sieve process:\n");
805
806
807 schain(SIEVE_CUT + 1); /* Do first piece of chain, for small-sieving. */
808 linit = 0;
809 ct = 0;
810 itog(1, magcheck);
811 for(el = 2; el <= L_MAX; el += 1 + (el%2)) {
812 if(!primeq(el)) continue;
813 L = el; while(L*el <= 4) L *= el;
814 printf("Resolving Schoof L = %d...\n", L);
815 P[ct] = L; /* Stuff another prime power. */
816 smulg(L, magcheck);
817 gtog(magcheck, tmp); squareg(tmp); gshiftright(2, tmp);
818 if(gcompg(tmp, p) > 0) break; /* Done...go to CRT phase. */
819 if((L > SIEVE_CUT) && (linit == 0)) {
820 schain(L_MAX+1);
821 linit = 1;
822 }
823 init_recip(L);
824 // printf("s: "); polyout(s[L]);
825 gtog(p, aux2);
826 k = idivg(L, aux2); /* p (mod L). */
827 gtog(p, aux2);
828 k2 = idivg(el, aux2);
829 // printf("power...\n");
830 txd->deg = 0; itog(1, txd->coe[0]);
831 tyd->deg = 0; itog(1, tyd->coe[0]);
832 txn->deg = 1; itog(0, txn->coe[0]); itog(1, txn->coe[1]);
833 ptop(txn, kxn);
834
835 gtog(p, aux2);
836 powerpoly(txn, aux2); /* x^p. */
837 // printf("x^p done...\n");
838 ptop(txn, powx);
839 powerpoly(powx, aux2);
840 // printf("x^p^2 done...\n");
841 ptop(cubic, tyn);
842 gtog(p, aux2); itog(1, aux); subg(aux, aux2);
843 gshiftright(1, aux2); /* aux2 := (p-1)/2. */
844 powerpoly(tyn, aux2); /* y^p. */
845 // printf("y^p done...\n");
846 ptop(tyn, powy); mulp(tyn, powy); fullmod(powy);
847 mulp(cubic, powy); fullmod(powy);
848 powerpoly(powy, aux2);
849 mulp(tyn, powy); fullmod(powy);
850 // printf("Powers done...\n");
851
852 // printf("pow" ); polyout(powx); polyout(powy);
853 ptop(txn, txn1); ptop(txd, txd1); /* Save t = 1 case. */
854 ptop(tyn, tyn1); ptop(tyd, tyd1);
855 /* We now shall test
856 (powx, y powy) + k(kxn/kxd, y kyn/kyd) = t(txn/txd, y tyn/tyd)
857 */
858
859 if(k==1) { ptop(txd, kxd); ptop(txd, kyd);
860 ptop(txd, kyn);
861 } else {
862 ptop(s[k], kxd); mulp(s[k], kxd); fullmod(kxd);
863 if(k%2==0) { mulp(cubic, kxd); fullmod(kxd); }
864 mulp(kxd, kxn); fullmod(kxn);
865 ptop(s[k-1], pbuff); mulp(s[k+1], pbuff); fullmod(pbuff);
866 if(k%2==1) {mulp(cubic, pbuff); fullmod(pbuff); }
867 subp(pbuff, kxn); fullmod(kxn);
868
869 ptop(s[k+2], kyn); mulp(s[k-1], kyn); fullmod(kyn);
870 mulp(s[k-1], kyn); fullmod(kyn);
871 if(k > 2) {
872 ptop(s[k-2], pbuff); mulp(s[k+1], pbuff); fullmod(pbuff);
873 mulp(s[k+1], pbuff); fullmod(pbuff);
874 subp(pbuff, kyn); fullmod(kyn);
875 }
876 ptop(s[k], kyd); mulp(s[k], kyd); fullmod(kyd);
877 mulp(s[k], kyd); fullmod(kyd);
878 if(k%2==0) { mulp(cubic, kyd); fullmod(kyd);
879 mulp(cubic, kyd); fullmod(kyd);}
880 itog(4, aux2); scalarmul(aux2, kyd);
881 }
882 //printf("kP: "); polyout(kxn); polyout(kxd); polyout(kyn); polyout(kyd);
883 /* Commence t = 0 check. */
884 // printf("Checking t = %d ...\n", 0);
885 fflush(stdout);
886
887 ptop(powx, pbuff); mulp(kxd, pbuff);
888 subp(kxn, pbuff);
889 fullmod(pbuff);
890
891 xmatch = ymatch = 0;
892 if(iszerop(pbuff)) {
893 xmatch = 1;
894 /* Now check y coords. */
895 if(L == 2) goto resolve;
896 ptop(powy, pbuff); mulp(kyd, pbuff);
897 addp(kyn, pbuff);
898 fullmod(pbuff);
899 if(iszerop(pbuff)) {
900 resolve:
901 printf("%d %d\n", L, 0);
902 T[ct++] = 0;
903 if((k2 + 1 - T[ct-1]) % el == 0) {
904 printf("TILT: %d\n", el);
905 el = 2;
906 iaddg(1, b);
907 goto newb;
908 }
909 continue;
910 } else ymatch = 1;
911 }
912 /* Combine pt1 and pt2. */
913 if((xmatch == 1) && (ymatch == 1))
914 elldoublep(powx, txd, powy, txd, nx, dx, ny, dy);
915 else
916 elladdp(powx, txd, powy, txd, kxn, kxd, kyn, kyd, nx, dx, ny, dy);
917
918 /* Now {nx/dx, ny/dy} is (fixed) LHS. */
919 // printf("add12: "); polyout(nx); polyout(dx); polyout(ny); polyout(dy);
920 /* Commence t > 0 check. */
921 for(t=1; t <= L/2; t++) {
922 // printf("Checking t = %d ...\n", t);
923 if(t > 1) { /* Add (tx1, ty1) to (txn, tyn). */
924 ptop(txn1, pbuff); mulmod(txd, pbuff);
925 ptop(txn, powx); mulmod(txd1, powx);
926 subp(powx, pbuff); fullmod(pbuff);
927 if(!iszerop(pbuff))
928 elladdp(txn1, txd1, tyn1, tyd1, txn, txd, tyn, tyd,
929 tmp1, tmp2, tmp3, tmp4);
930 else elldoublep(txn, txd, tyn, tyd,
931 tmp1, tmp2, tmp3, tmp4);
932 ptop(tmp1, txn); ptop(tmp2, txd);
933 ptop(tmp3, tyn); ptop(tmp4, tyd);
934 }
935 // printf("tQ: "); polyout(txn); polyout(txd); polyout(tyn); polyout(tyd);
936 /* Next, check {nx/dx, ny/dy} =? {txn/txd, tyn/tyd}. */
937 ptop(nx, pbuff); mulmod(txd, pbuff);
938 ptop(dx, powx); mulmod(txn, powx);
939 subp(powx, pbuff); fullmod(pbuff);
940 if(!iszerop(pbuff)) continue;
941 /* Next, check y. */
942 // printf("y check!\n");
943 ptop(ny, pbuff); mulmod(tyd, pbuff);
944 ptop(dy, powx); mulmod(tyn, powx);
945 subp(powx, pbuff); fullmod(pbuff);
946 if(iszerop(pbuff)) {
947 printf("%d %d\n", L, t);
948 T[ct++] = t;
949 } else {
950 printf("%d %d\n", L, L-t);
951 T[ct++] = L-t;
952 }
953 if((k2 + 1 - T[ct-1]) % el == 0) {
954 printf("TILT: %d\n", el);
955 el = 2;
956 iaddg(1, b);
957 goto newb;
958 }
959
960 fflush(stdout);
961 break;
962 }
963 }
964
965 /* Now, prime powers P[] and CRT residues T[] are intact. */
966 printf("Prime powers L:\n");
967 printf("{");
968 for(j=0; j < ct-1; j++) {
969 printf("%d, ", P[j]);
970 }
971 printf("%d }\n", P[ct-1]);
972
973 printf("Residues t (mod L):\n");
974 printf("{");
975 for(j=0; j < ct-1; j++) {
976 printf("%d, ", T[j]);
977 }
978 printf("%d }\n", T[ct-1]);
979
980 /* Mathematica algorithm for order:
981 plis = {2^5, 3^3, 5^2, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47};
982 tlis = {1, 26, 4, 2, 4, 11, 6, 5, 19, 22, 10, 16, 7, 22, 11};
983 prod = Apply[Times, plis];
984 prlis = prod/plis;
985 invlis = Table[PowerMod[prlis[[q]], -1, plis[[q]]],{q,1,Length[plis]}];
986 p = 2^127 - 1;
987 t = Mod[tlis . (prlis * invlis), prod];
988 ord = p + 1 - If[t^2 > 4p, t - prod, t]
989 */
990
991 itog(1, t1);
992 for(j=0; j < ct; j++) {
993 smulg(P[j], t1);
994 }
995
996 for(j=0; j < 2*ct; j++) {
997 if(!ss[j]) ss[j] = newgiant(MAX_DIGS);
998 }
999
1000 for(j=0; j < ct; j++) {
1001 gtog(t1, ss[j]);
1002 itog(P[j], t2);
1003 divg(t2, ss[j]);
1004 }
1005
1006 for(j=0; j < ct; j++) {
1007 gtog(ss[j], ss[j+ct]);
1008 itog(P[j], t2);
1009 invg(t2, ss[j+ct]);
1010 }
1011
1012 itog(0, t4);
1013 for(j=0; j < ct; j++) {
1014 itog(T[j], t5);
1015 mulg(ss[j], t5);
1016 mulg(ss[j+ct], t5);
1017 addg(t5, t4);
1018 }
1019 modg(t1, t4);
1020 gtog(p, t5);
1021 iaddg(1, t5);
1022 gtog(t4, t2);
1023 squareg(t4);
1024 gtog(p, t3); gshiftleft(2, t3);
1025 if(gcompg(t4, t3) > 0) subg(t1, t2);
1026 subg(t2, t5);
1027 printf("Parameters:\n");
1028 printf("p = "); gout(p);
1029 printf("a = "); gout(a);
1030 printf("b = "); gout(b);
1031 printf("Curve order:\n");
1032 printf("o = "); gout(t5);
1033 printf("pprob: %d\n", prime_probable(t5));
1034 printf("Twist order:\n");
1035 printf("o' = ");
1036 addg(t2, t5);
1037 addg(t2, t5);
1038 gout(t5);
1039 printf("pprob: %d\n", prime_probable(t5));
1040
1041 iaddg(1,b);
1042 goto newb;
1043 }
mirror of Security