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