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