![[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: