Annotation of OpenXM_contrib2/asir2000/builtin/ec.c, Revision 1.5
1.3 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.5 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/ec.c,v 1.4 2000/08/22 05:03:57 noro Exp $
1.3 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52: #include "inline.h"
53:
54: /*
55: * Elliptic curve related functions
56: *
57: * Argument specifications
58: * Point = vector(x,y,z) (projective representation)
59: * ECInfo = vector(a,b,p)
60: * a,b : integer for GF(p), GF2N for GF(2^n)
61: * p : integer for GF(p), polynomial for GF(2^n)
62: *
63: */
64:
65: struct oKeyIndexPair {
66: unsigned int key;
67: int index;
68: };
69:
70: void Psha1();
71: void Psha1_ec();
72: void Pecm_add_ff();
73: void Pecm_chsgn_ff();
74: void Pecm_sub_ff();
75: void Pecm_compute_all_key_homo_ff();
76: void Pnextvect1(),Psort_ktarray(),Pecm_find_match(),Pseparate_vect();
77: void Pecm_set_addcounter();
78: void Pecm_count_order();
79:
80: void ecm_add_ff(VECT,VECT,VECT,VECT *);
81: void ecm_add_gfp(VECT,VECT,VECT,VECT *);
82: void ecm_add_gf2n(VECT,VECT,VECT,VECT *);
83:
84: void ecm_chsgn_ff(VECT,VECT *);
85: void ecm_sub_ff(VECT,VECT,VECT,VECT *);
86:
87: void compute_all_key_homo_gfp(VECT *,int,unsigned int *);
88: void compute_all_key_homo_gf2n(VECT *,int,unsigned int *);
89:
90: unsigned int separate_vect(double *,int);
91: void ecm_find_match(unsigned int *,int,unsigned int *,int,LIST *);
92: int find_match(unsigned int,unsigned int *,int);
93:
94: void sort_ktarray(VECT,VECT,LIST *);
95: int comp_kip(struct oKeyIndexPair *,struct oKeyIndexPair *);
96:
97: int nextvect1(VECT,VECT);
98:
99: unsigned int ecm_count_order_gfp(unsigned int,unsigned int,unsigned int);
100: unsigned int ecm_count_order_gf2n(UP2,GF2N,GF2N);
101:
102: void sha_memory(unsigned char *,unsigned int,unsigned int *);
103:
104: unsigned int ecm_addcounter;
105: extern int current_ff,lm_lazy;
106:
107: struct ftab ec_tab[] = {
108: /* point addition */
109: {"ecm_add_ff",Pecm_add_ff,3},
110:
111: /* point change sign */
112: {"ecm_chsgn_ff",Pecm_chsgn_ff,1},
113:
114: /* point subtraction */
115: {"ecm_sub_ff",Pecm_sub_ff,3},
116:
117: /* key computation for sort and match */
118: {"ecm_compute_all_key_homo_ff",Pecm_compute_all_key_homo_ff,1},
119:
120: /* exhausitve search of rational points */
121: {"ecm_count_order",Pecm_count_order,1},
122:
123: /* common functions */
124: {"nextvect1",Pnextvect1,2},
125: {"sort_ktarray",Psort_ktarray,2},
126: {"separate_vect",Pseparate_vect,1},
127: {"ecm_find_match",Pecm_find_match,2},
128: {"ecm_set_addcounter",Pecm_set_addcounter,-1},
1.2 noro 129: {"sha1",Psha1,-2},
1.1 noro 130: #if 0
131: {"sha1_free",Psha1_free,1},
132: #endif
133: {0,0,0},
134: };
135:
136: void Psha1(arg,rp)
137: NODE arg;
138: Q *rp;
139: {
140: unsigned char *s;
141: unsigned int digest[5];
142: int i,j,l,bl,n;
143: unsigned int t;
144: N z,r;
145:
146: asir_assert(ARG0(arg),O_N,"sha1_free");
147: z = NM((Q)ARG0(arg));
148: n = PL(z);
149: l = n_bits(z);
1.2 noro 150: if ( argc(arg) == 2 )
151: bl = QTOS((Q)ARG1(arg));
152: else
153: bl = (l+7)/8;
1.1 noro 154: s = (unsigned char *)MALLOC(bl);
155: for ( i = 0, j = bl-1; i < n; i++ ) {
156: t = BD(z)[i];
157: if ( j >= 0 ) s[j--] = t&0xff; t>>=8;
158: if ( j >= 0 ) s[j--] = t&0xff; t>>=8;
159: if ( j >= 0 ) s[j--] = t&0xff; t>>=8;
160: if ( j >= 0 ) s[j--] = t;
161: }
162: sha_memory(s,bl,digest);
163: r = NALLOC(5);
164: for ( i = 0; i < 5; i++ )
165: BD(r)[i] = digest[4-i];
166: for ( i = 4; i >= 0 && !BD(r)[i]; i-- );
167: if ( i < 0 )
168: *rp = 0;
169: else {
170: PL(r) = i+1;
171: NTOQ(r,1,*rp);
172: }
173: }
174:
175: #if 0
176: void Psha1_ec(arg,rp)
177: NODE arg;
178: Q *rp;
179: {
180: #include <fj_crypt.h>
181: SHS_CTX context;
182: unsigned char *s;
183: int i,j,l,bl,n;
184: unsigned int t;
185: N z,r;
186: extern int little_endian;
187:
188: asir_assert(ARG0(arg),O_N,"sha1");
189: z = NM((Q)ARG0(arg));
190: n = PL(z);
191: l = n_bits(z);
192: bl = (l+7)/8;
193: s = (unsigned char *)MALLOC(bl);
194: for ( i = 0, j = bl-1; i < n; i++ ) {
195: t = BD(z)[i];
196: if ( j >= 0 ) s[j--] = t&0xff; t>>=8;
197: if ( j >= 0 ) s[j--] = t&0xff; t>>=8;
198: if ( j >= 0 ) s[j--] = t&0xff; t>>=8;
199: if ( j >= 0 ) s[j--] = t;
200: }
201: SHSInit(&context);
202: SHSUpdate(&context,s,bl);
203: SHSFinal(&context);
204: r = NALLOC(5);
205: if ( little_endian )
206: for ( i = 0; i < 5; i++ ) {
207: t = context.digest[4-i];
208: BD(r)[i] = (t>>24)|((t&0xff0000)>>8)|((t&0xff00)<<8)|(t<<24);
209: } else
210: for ( i = 0; i < 5; i++ )
211: BD(r)[i] = context.digest[4-i];
212: for ( i = 4; i >= 0 && !BD(r)[i]; i-- );
213: if ( i < 0 )
214: *rp = 0;
215: else {
216: PL(r) = i+1;
217: NTOQ(r,1,*rp);
218: }
219: }
220: #endif
221:
222: void Pecm_count_order(arg,rp)
223: NODE arg;
224: Q *rp;
225: {
226: N p;
227: UP2 d;
228: Obj a,b;
229: unsigned int p0,a0,b0,ord;
230: VECT ec;
231: Obj *vb;
232:
233: switch ( current_ff ) {
234: case FF_GFP:
235: getmod_lm(&p);
236: if ( n_bits(p) > 32 )
237: error("ecm_count_order : ground field too large");
238: p0 = BD(p)[0];
239: ec = (VECT)ARG0(arg);
240: vb = (Obj *)BDY(ec); simp_ff(vb[0],&a); simp_ff(vb[1],&b);
241: a0 = a?BD(BDY((LM)a))[0]:0;
242: b0 = b?BD(BDY((LM)b))[0]:0;
243: ord = ecm_count_order_gfp(p0,a0,b0);
244: UTOQ(ord,*rp);
245: break;
246: case FF_GF2N:
247: getmod_gf2n(&d);
248: if ( degup2(d) > 10 )
249: error("ecm_count_order : ground field too large");
250: ec = (VECT)ARG0(arg);
251: vb = (Obj *)BDY(ec); simp_ff(vb[0],&a); simp_ff(vb[1],&b);
252: ord = ecm_count_order_gf2n(d,(GF2N)a,(GF2N)b);
253: UTOQ(ord,*rp);
254: break;
255: default:
256: error("ecm_count_order : current_ff is not set");
257: }
258: }
259:
260: void Pecm_set_addcounter(arg,rp)
261: NODE arg;
262: Q *rp;
263: {
264: if ( arg )
265: ecm_addcounter = QTOS((Q)ARG0(arg));
266: UTOQ(ecm_addcounter,*rp);
267: }
268:
269: void Pecm_compute_all_key_homo_ff(arg,rp)
270: NODE arg;
271: VECT *rp;
272: {
273: unsigned int *ka;
274: int len,i;
275: VECT *pa;
276: VECT r,v;
277: LIST *vb;
278: USINT *b;
279:
280: v = (VECT)ARG0(arg);
281: len = v->len;
282: vb = (LIST *)v->body;
283: pa = (VECT *)ALLOCA(len*sizeof(VECT));
284: ka = (unsigned int *)ALLOCA(len*sizeof(unsigned int));
285: for ( i = 0; i < len; i++ )
286: pa[i] = (VECT)BDY(NEXT(BDY(vb[i])));
287: switch ( current_ff ) {
288: case FF_GFP:
289: compute_all_key_homo_gfp(pa,len,ka); break;
290: case FF_GF2N:
291: compute_all_key_homo_gf2n(pa,len,ka); break;
292: default:
293: error("ecm_compute_all_key_homo_ff : current_ff is not set");
294: }
295: MKVECT(r,len); *rp = r;
296: b = (USINT *)r->body;
297: for ( i = 0; i < len; i++ )
298: MKUSINT(b[i],ka[i]);
299: }
300:
301: void Psort_ktarray(arg,rp)
302: NODE arg;
303: LIST *rp;
304: {
305: sort_ktarray((VECT)ARG0(arg),(VECT)ARG1(arg),rp);
306: }
307:
308: void Pecm_add_ff(arg,rp)
309: NODE arg;
310: VECT *rp;
311: {
312: ecm_add_ff(ARG0(arg),ARG1(arg),ARG2(arg),rp);
313: }
314:
315: void Pecm_sub_ff(arg,rp)
316: NODE arg;
317: VECT *rp;
318: {
319: ecm_sub_ff(ARG0(arg),ARG1(arg),ARG2(arg),rp);
320: }
321:
322: void Pecm_chsgn_ff(arg,rp)
323: NODE arg;
324: VECT *rp;
325: {
326: ecm_chsgn_ff(ARG0(arg),rp);
327: }
328:
329: void Pnextvect1(arg,rp)
330: NODE arg;
331: Q *rp;
332: {
333: int index;
334:
335: index = nextvect1(ARG0(arg),ARG1(arg));
336: STOQ(index,*rp);
337: }
338:
339: /* XXX at least n < 32 must hold. What is the strict restriction for n ? */
340:
341: void Pseparate_vect(arg,rp)
342: NODE arg;
343: LIST *rp;
344: {
345: VECT v;
346: int n,i;
347: Q *b;
348: double *w;
349: unsigned int s;
350: NODE ns,nc,t,t1;
351: Q iq;
352: LIST ls,lc;
353:
354: v = (VECT)ARG0(arg);
355: n = v->len; b = (Q *)v->body;
356: w = (double *)ALLOCA(n*sizeof(double));
357: for ( i = 0; i < n; i++ )
358: w[i] = (double)QTOS(b[i]);
359: s = separate_vect(w,n);
360: ns = nc = 0;
361: for ( i = n-1; i >= 0; i-- )
362: if ( s & (1<<i) ) {
363: STOQ(i,iq); MKNODE(t,iq,ns); ns = t;
364: } else {
365: STOQ(i,iq); MKNODE(t,iq,nc); nc = t;
366: }
367: MKLIST(ls,ns); MKLIST(lc,nc);
368: MKNODE(t,lc,0); MKNODE(t1,ls,t);
369: MKLIST(*rp,t1);
370: }
371:
372: void Pecm_find_match(arg,rp)
373: NODE arg;
374: LIST *rp;
375: {
376: VECT g,b;
377: int ng,nb,i;
378: USINT *p;
379: unsigned int *kg,*kb;
380:
381: g = (VECT)ARG0(arg); ng = g->len;
382: kg = (unsigned int *)ALLOCA(ng*sizeof(unsigned int));
383: for ( i = 0, p = (USINT *)g->body; i < ng; i++ )
384: kg[i] = p[i]?BDY(p[i]):0;
385: b = (VECT)ARG1(arg); nb = b->len;
386: kb = (unsigned int *)ALLOCA(nb*sizeof(unsigned int));
387: for ( i = 0, p = (USINT *)b->body; i < nb; i++ )
388: kb[i] = p[i]?BDY(p[i]):0;
389: ecm_find_match(kg,ng,kb,nb,rp);
390: }
391:
392: void ecm_add_ff(p1,p2,ec,pr)
393: VECT p1,p2,ec;
394: VECT *pr;
395: {
396: if ( !p1 )
397: *pr = p2;
398: else if ( !p2 )
399: *pr = p1;
400: else {
401: switch ( current_ff ) {
402: case FF_GFP:
403: ecm_add_gfp(p1,p2,ec,pr); break;
404: case FF_GF2N:
405: ecm_add_gf2n(p1,p2,ec,pr); break;
406: default:
407: error("ecm_add_ff : current_ff is not set");
408: }
409: }
410: }
411:
412: /* ec = [AX,BC] */
413:
414: void ecm_add_gf2n(p1,p2,ec,rp)
415: VECT p1,p2,ec;
416: VECT *rp;
417: {
418: GF2N ax,bc,a0,a1,a2,b0,b1,b2;
419: GF2N a2b0,a0b2,a2b1,a1b2,a02,a04,a22,a24,a0a2,a0a22,a1a2;
1.5 ! noro 420: GF2N t,s,u,r0,r1,r00,r01,r02,r002,r003,r02q;
1.1 noro 421: VECT r;
422: GF2N *vb,*rb;
423:
424: ecm_addcounter++;
425: /* addition with O */
426: if ( !p1 ) {
427: *rp = p2;
428: return;
429: }
430: if ( !p2 ) {
431: *rp = p1;
432: return;
433: }
434: vb = (GF2N *)BDY(ec); ax = vb[0]; bc = vb[1];
435: vb = (GF2N *)BDY(p1); a0 = vb[0]; a1 = vb[1]; a2 = vb[2];
436: vb = (GF2N *)BDY(p2); b0 = vb[0]; b1 = vb[1]; b2 = vb[2];
437:
438: mulgf2n(a2,b0,&a2b0); mulgf2n(a0,b2,&a0b2);
439: if ( !cmpgf2n(a2b0,a0b2) ) {
440: mulgf2n(a2,b1,&a2b1);
441: mulgf2n(a1,b2,&a1b2);
442: if ( !cmpgf2n(a2b1,a1b2) ) {
443: if ( !a0 )
444: *rp = 0;
445: else {
446: squaregf2n(a0,&a02); squaregf2n(a02,&a04);
447: squaregf2n(a2,&a22); squaregf2n(a22,&a24);
448: mulgf2n(a0,a2,&a0a2); squaregf2n(a0a2,&a0a22);
449: mulgf2n(bc,a24,&t); addgf2n(a04,t,&r0);
450: mulgf2n(a04,a0a2,&t); mulgf2n(a1,a2,&a1a2);
451: addgf2n(a02,a1a2,&s); addgf2n(s,a0a2,&u);
452: mulgf2n(u,r0,&s); addgf2n(t,s,&r1);
453:
454: MKVECT(r,3); rb = (GF2N *)r->body;
455: mulgf2n(r0,a0a2,&rb[0]); rb[1] = r1; mulgf2n(a0a22,a0a2,&rb[2]);
456: *rp = r;
457: }
458: } else
459: *rp = 0;
460: } else {
461: mulgf2n(a1,b2,&a1b2); addgf2n(a0b2,a2b0,&r00);
462: mulgf2n(a2,b1,&t); addgf2n(a1b2,t,&r01); mulgf2n(a2,b2,&r02);
463: squaregf2n(r00,&r002); mulgf2n(r002,r00,&r003);
464:
465: addgf2n(r00,r01,&t); mulgf2n(t,r01,&s); mulgf2n(s,r02,&t);
466: if ( ax ) {
467: mulgf2n(r02,ax,&r02q);
468: addgf2n(t,r003,&s); mulgf2n(r02q,r002,&t); addgf2n(s,t,&r0);
469: } else
470: addgf2n(t,r003,&r0);
471:
472: mulgf2n(a0b2,r002,&t); addgf2n(t,r0,&s); mulgf2n(r01,s,&t);
473: mulgf2n(r002,a1b2,&s); addgf2n(r0,s,&u); mulgf2n(r00,u,&s);
474: addgf2n(t,s,&r1);
475:
476: MKVECT(r,3); rb = (GF2N *)r->body;
477: mulgf2n(r0,r00,&rb[0]); rb[1] = r1; mulgf2n(r003,r02,&rb[2]);
478: *rp = r;
479: }
480: }
481:
482: extern LM THREE_LM,FOUR_LM,EIGHT_LM;
483:
484: /* 0 < p < 2^16, 0 <= a,b < p */
485:
486: unsigned int ecm_count_order_gfp(p,a,b)
487: unsigned int p,a,b;
488: {
1.5 ! noro 489: unsigned int x,rhs,ord,t;
1.1 noro 490:
491: for ( x = 0, ord = 1; x < p; x++ ) {
492: DMAR(x,x,a,p,t) /* t = x^2+a mod p */
493: DMAR(t,x,b,p,rhs) /* rhs = x*(x^2+a)+b mod p */
494: if ( !rhs )
495: ord++;
496: else if ( small_jacobi(rhs,p)==1 )
497: ord+=2;
498: }
499: return ord;
500: }
501:
502: unsigned int ecm_count_order_gf2n(d,a,b)
503: UP2 d;
504: GF2N a,b;
505: {
506: error("ecm_count_order_gf2n : not implemented yet");
1.5 ! noro 507: /* NOTREACHED */
! 508: return 0;
1.1 noro 509: }
510:
511: /* ec = [AX,BC] */
512:
513: void ecm_add_gfp(p1,p2,ec,pr)
514: VECT p1,p2,ec;
515: VECT *pr;
516: {
517: LM aa,bb,x1,y1,z1,x2,y2,z2,x1z2,v1,y1z2,u1,u2,v2,v3,z1z2;
518: LM v2x1z2,a1,x3,y3,z3,w1,s1,s2,s3,s1y1,b1,h1;
519: LM t,s,u;
520: LM *vb;
521: VECT r;
522:
523: ecm_addcounter++;
524: /* addition with O */
525: if( !p1 ) {
526: *pr = p2;
527: return;
528: }
529: if( !p2 ) {
530: *pr = p1;
531: return;
532: }
533:
534: /* set parameters */
535:
536: /* aa = ec[0]; bb = ec[1]; */
537: vb = (LM *)BDY(ec); aa = vb[0]; bb = vb[1];
538:
539: /* x1 = p1[0]; y1 = p1[1]; z1 = p1[2]; */
540: vb = (LM *)BDY(p1); x1 = vb[0]; y1 = vb[1]; z1 = vb[2];
541:
542: /* x2 = p2[0]; y2 = p2[1]; z2 = p2[2]; */
543: vb = (LM *)BDY(p2); x2 = vb[0]; y2 = vb[1]; z2 = vb[2];
544:
545: /* addition */
546:
547: /* x1z2 = (x1*z2) %p; */
548: mullm(x1,z2,&x1z2);
549:
550: /* v1 = (x2*z1-x1z2) %p; */
551: lm_lazy = 1;
552: mullm(x2,z1,&t); sublm(t,x1z2,&s);
553: lm_lazy = 0; simplm(s,&v1);
554:
555: /* y1z2 = (y1*z2) %p; */
556: mullm(y1,z2,&y1z2);
557:
558: /* u1 = (y2*z1-y1z2) %p; */
559: lm_lazy = 1;
560: mullm(y2,z1,&t); sublm(t,y1z2,&s);
561: lm_lazy = 0; simplm(s,&u1);
562:
563: if( v1 != 0 ) {
564: /* u2 = (u1*u1) %p; */
565: mullm(u1,u1,&u2);
566:
567: /* v2 = (v1*v1) %p; */
568: mullm(v1,v1,&v2);
569:
570: /* v3 = (v1*v2) %p; */
571: mullm(v1,v2,&v3);
572:
573: /* z1z2 = (z1*z2) %p; */
574: mullm(z1,z2,&z1z2);
575:
576: /* v2x1z2 = (v2*x1z2) %p; */
577: mullm(v2,x1z2,&v2x1z2);
578:
579: /* a1 = (u2*z1z2-v3-2*v2x1z2) %p; */
580: lm_lazy = 1;
581: mullm(u2,z1z2,&t); addlm(v2x1z2,v2x1z2,&s);
582: addlm(v3,s,&u); sublm(t,u,&s);
583: lm_lazy = 0; simplm(s,&a1);
584:
585: /* x3 = ( v1 * a1 ) %p; */
586: mullm(v1,a1,&x3);
587:
588: /* y3 = ( u1 * ( v2x1z2 - a1 ) - v3 * y1z2 ) %p; */
589: lm_lazy = 1;
590: sublm(v2x1z2,a1,&t); mullm(u1,t,&s); mullm(v3,y1z2,&u); sublm(s,u,&t);
591: lm_lazy = 0; simplm(t,&y3);
592:
593: /* z3 = ( v3 * z1z2 ) %p; */
594: mullm(v3,z1z2,&z3);
595: } else if( u1 == 0 ) {
596: /* w1 = (aa*z1*z1+3*x1*x1) %p; */
597: lm_lazy = 1;
598: mullm(z1,z1,&t); mullm(aa,t,&s);
599: mullm(x1,x1,&t); mullm(THREE_LM,t,&u); addlm(s,u,&t);
600: lm_lazy = 0; simplm(t,&w1);
601:
602: /* s1 = (y1*z1) %p; */
603: mullm(y1,z1,&s1);
604:
605: /* s2 = (s1*s1) %p; */
606: mullm(s1,s1,&s2);
607:
608: /* s3 = (s1*s2) %p; */
609: mullm(s1,s2,&s3);
610:
611: /* s1y1 = (s1*y1) %p; */
612: mullm(s1,y1,&s1y1);
613:
614: /* b1 = (s1y1*x1) %p; */
615: mullm(s1y1,x1,&b1);
616:
617: /* h1 = (w1*w1-8*b1) %p; */
618: lm_lazy = 1;
619: mullm(w1,w1,&t); mullm(EIGHT_LM,b1,&s); sublm(t,s,&u);
620: lm_lazy = 0; simplm(u,&h1);
621:
622: /* x3 = ( 2 * h1 * s1 ) %p; */
623: lm_lazy = 1;
624: mullm(h1,s1,&t); addlm(t,t,&s);
625: lm_lazy = 0; simplm(s,&x3);
626:
627: /* y3 = ( w1 * ( 4 * b1 - h1 ) - 8 * s1y1 * s1y1 ) %p; */
628: lm_lazy = 1;
629: mullm(FOUR_LM,b1,&t); sublm(t,h1,&s); mullm(w1,s,&u);
630: mullm(s1y1,s1y1,&t); mullm(EIGHT_LM,t,&s); sublm(u,s,&t);
631: lm_lazy = 0; simplm(t,&y3);
632:
633: /* z3 = ( 8 * s3 ) %p; */
634: mullm(EIGHT_LM,s3,&z3);
635: } else {
636: *pr = 0;
637: return;
638: }
639: if ( !z3 )
640: *pr = 0;
641: else {
642: MKVECT(r,3); *pr = r;
643: vb = (LM *)BDY(r); vb[0] = x3; vb[1] = y3; vb[2] = z3;
644: }
645: }
646:
647: void ecm_chsgn_ff(p,pr)
648: VECT p;
649: VECT *pr;
650: {
1.5 ! noro 651: Obj x,y,z;
1.1 noro 652: LM tl;
653: GF2N tg;
654: Obj *vb;
655: VECT r;
656:
657: if( !p ) {
658: *pr = 0;
659: return;
660: }
661:
662: /* x = p[0]; y = p[1]; z = p[2]; */
663: vb = (Obj *)BDY(p); x = vb[0]; y = vb[1]; z = vb[2];
664: switch ( current_ff ) {
665: case FF_GFP:
666: if ( !y )
667: *pr = p;
668: else {
669: chsgnlm((LM)y,&tl);
670: MKVECT(r,3); *pr = r;
671: vb = (Obj *)BDY(r); vb[0] = x; vb[1] = (Obj)tl; vb[2] = z;
672: }
673: break;
674: case FF_GF2N:
675: addgf2n((GF2N)x,(GF2N)y,&tg);
676: MKVECT(r,3); *pr = r;
677: vb = (Obj *)BDY(r); vb[0] = x; vb[1] = (Obj)tg; vb[2] = z;
678: break;
679: default:
680: error("ecm_chsgn_ff : current_ff is not set");
681: }
682: }
683:
684: void ecm_sub_ff(p1,p2,ec,pr)
685: VECT p1,p2,ec;
686: VECT *pr;
687: {
688: VECT mp2;
689:
690: ecm_chsgn_ff(p2,&mp2);
691: ecm_add_ff(p1,mp2,ec,pr);
692: }
693:
694: /* tplist = [[t,p],...]; t:interger, p=[p0,p1]:point (vector) */
695:
696: int comp_kip(a,b)
697: struct oKeyIndexPair *a,*b;
698: {
699: unsigned int ka,kb;
700:
701: ka = a->key; kb = b->key;
702: if ( ka > kb )
703: return 1;
704: else if ( ka < kb )
705: return -1;
706: else
707: return 0;
708: }
709:
710: #define EC_GET_XZ(p,x,z) \
711: if ( !(p) ) {\
712: (x)=0; (z)=(LM)ONE;\
713: } else { \
714: LM *vb;\
715: vb = (LM *)BDY((VECT)(p));\
716: (x) = vb[0]; (z) = vb[2];\
717: }
718:
719: #define EC_GET_XZ_GF2N(p,x,z) \
720: if ( !(p) ) {\
721: (x)=0; (z)=(GF2N)ONE;\
722: } else { \
723: GF2N *vb;\
724: vb = (GF2N *)BDY((VECT)(p));\
725: (x) = vb[0]; (z) = vb[2];\
726: }
727:
728: void compute_all_key_homo_gfp(pa,len,ka)
729: VECT *pa;
730: int len;
731: unsigned int *ka;
732: {
733: LM *b,*x,*z;
734: int i;
735: LM t,s,m;
736:
737: b = (LM *)ALLOCA((len+1)*sizeof(LM));
738: x = (LM *)ALLOCA(len*sizeof(LM));
739: z = (LM *)ALLOCA(len*sizeof(LM));
740: MKLM(ONEN,b[0]);
741: for ( i = 1; i <= len; i++ ) {
742: EC_GET_XZ(pa[i-1],x[i-1],z[i-1]);
743: mullm(b[i-1],z[i-1],&b[i]);
744: }
745: /* b[0] = 1 */
746: divlm(b[0],b[len],&m);
747: for ( i = len-1; i >= 0; i-- ) {
748: mullm(m,b[i],&s); mullm(s,x[i],&t); s = t;
749: ka[i] = s ? s->body->b[0] : 0; ka[i] |= 0x80000000;
750: mullm(m,z[i],&s); m = s;
751: }
752: }
753:
754: void compute_all_key_homo_gf2n(pa,len,ka)
755: VECT *pa;
756: int len;
757: unsigned int *ka;
758: {
759: GF2N *b,*x,*z;
760: int i;
761: GF2N t,s,m;
762:
763: b = (GF2N *)ALLOCA((len+1)*sizeof(Q));
764: x = (GF2N *)ALLOCA(len*sizeof(Q));
765: z = (GF2N *)ALLOCA(len*sizeof(Q));
766: MKGF2N(ONEUP2,b[0]);
767: for ( i = 1; i <= len; i++ ) {
768: EC_GET_XZ_GF2N(pa[i-1],x[i-1],z[i-1]);
769: mulgf2n(b[i-1],z[i-1],&b[i]);
770: }
771: invgf2n(b[len],&m);
772: for ( i = len-1; i >= 0; i-- ) {
773: mulgf2n(m,b[i],&s); mulgf2n(s,x[i],&t); s = t;
774: ka[i] = s ? s->body->b[0] : 0; ka[i] |= 0x80000000;
775: mulgf2n(m,z[i],&s); m = s;
776: }
777: }
778:
779: unsigned int separate_vect(v,n)
780: double *v;
781: int n;
782: {
783: unsigned int max = 1<<n;
784: unsigned int i,j,i0;
785: double all,a,total,m;
786:
787: for ( i = 0, all = 1; i < (unsigned int)n; i++ )
788: all *= v[i];
789:
790: for ( i = 0, m = 0; i < max; i++ ) {
791: for ( a = 1, j = 0; j < (unsigned int)n; j++ )
792: if ( i & (1<<j) )
793: a *= v[j];
794: total = a+(all/a)*2;
795: if ( !m || total < m ) {
796: m = total;
797: i0 = i;
798: }
799: }
800: return i0;
801: }
802:
803: void ecm_find_match(g,ng,b,nb,r)
804: unsigned int *g;
805: int ng;
806: unsigned int *b;
807: int nb;
808: LIST *r;
809: {
810: int i,j;
811: Q iq,jq;
812: NODE n0,n1,c0,c;
813: LIST l;
814:
815: for ( i = 0, c0 = 0; i < ng; i++ ) {
816: j = find_match(g[i],b,nb);
817: if ( j >= 0 ) {
818: STOQ(i,iq); STOQ(j,jq);
819: MKNODE(n1,jq,0); MKNODE(n0,iq,n1); MKLIST(l,n0);
820: NEXTNODE(c0,c);
821: BDY(c) = (pointer)l;
822: }
823: }
824: if ( c0 )
825: NEXT(c) = 0;
826: MKLIST(*r,c0);
827: }
828:
829: int find_match(k,key,n)
830: unsigned int k;
831: unsigned int *key;
832: int n;
833: {
834: int s,e,m;
835:
836: for ( s = 0, e = n; (e-s) > 1; ) {
837: m = (s+e)/2;
838: if ( k==key[m] )
839: return m;
840: else if ( k > key[m] )
841: s = m;
842: else
843: e = m;
844: }
845: if ( k == key[s] )
846: return s;
847: else
848: return -1;
849: }
850:
851: int nextvect1(vect,bound)
852: VECT vect,bound;
853: {
854: int size,i,a;
855: Q *vb,*bb;
856:
857: size = vect->len;
858: vb = (Q *)vect->body;
859: bb = (Q *)bound->body;
860: for ( i = size-1; i >= 0; i-- )
861: if ( (a=QTOS(vb[i])) < QTOS(bb[i]) ) {
862: a++; STOQ(a,vb[i]);
863: break;
864: } else
865: vb[i] = 0;
866: return i;
867: }
868:
869: void sort_ktarray(karray,tarray,rp)
870: VECT karray,tarray;
871: LIST *rp;
872: {
873: NODE r,r1;
874: int i,i0,k,len,same,tsame;
875: struct oKeyIndexPair *kip;
876: VECT key,value,v;
877: Q *tb,*samebuf;
878: USINT *kb;
879: Obj *svb;
880: USINT *skb;
881:
882: len = karray->len;
883: kb = (USINT *)karray->body;
884:
885: kip = (struct oKeyIndexPair *)ALLOCA(len*sizeof(struct oKeyIndexPair));
886: for ( i = 0; i < len; i++ ) {
887: kip[i].key = BDY(kb[i]); kip[i].index = i;
888: }
889: qsort((void *)kip,len,sizeof(struct oKeyIndexPair),
890: (int (*)(const void *,const void *))comp_kip);
891:
892: for ( same = tsame = i = i0 = 0, k = 1; i < len; i++, tsame++ )
893: if ( kip[i0].key != kip[i].key ) {
894: i0 = i; k++;
895: same = MAX(tsame,same);
896: tsame = 0;
897: }
898: same = MAX(tsame,same);
899: samebuf = (Q *)ALLOCA(same*sizeof(Q));
900:
901: MKVECT(key,k); skb = (USINT *)BDY(key);
902: MKVECT(value,k); svb = (Obj *)BDY(value);
903:
904: tb = (Q *)tarray->body;
905: for ( same = i = i0 = k = 0; i <= len; i++ ) {
906: if ( i == len || kip[i0].key != kip[i].key ) {
907: skb[k] = kb[kip[i0].index];
908: if ( same > 1 ) {
909: MKVECT(v,same);
910: bcopy((char *)samebuf,(char *)v->body,same*sizeof(Q));
911: svb[k] = (Obj)v;
912: } else
913: svb[k] = (Obj)samebuf[0];
914: i0 = i;
915: k++;
916: same = 0;
917: if ( i == len )
918: break;
919: }
920: samebuf[same++] = tb[kip[i].index];
921: }
922: MKNODE(r1,value,0); MKNODE(r,key,r1); MKLIST(*rp,r);
923: }
924:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>