![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to upo.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to upo.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/upo.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3: #include "base.h"
! 4:
! 5: void find_root_up2();
! 6:
! 7: #define INLINE
! 8:
! 9: #if defined(VISUAL)
! 10: #undef INLINE
! 11: #define INLINE __inline
! 12: #endif
! 13:
! 14: #if defined(linux)
! 15: #undef INLINE
! 16: #define INLINE inline
! 17: #endif
! 18:
! 19: static int up2_bbtab_initialized;
! 20: static unsigned short up2_sqtab[256];
! 21: static unsigned short up2_bbtab[65536];
! 22:
! 23: int up2_kara_mag = 10;
! 24: V up2_var;
! 25:
! 26: extern GEN_UP2 current_mod_gf2n;
! 27: extern int lm_lazy;
! 28:
! 29: #define GCD_COEF(q,p1,p2,p3,q1,q2,q3)\
! 30: switch (q){\
! 31: case 2:(p3)=(p1)^((p2)<<1);(q3)=(q1)^((q2)<<1);break;\
! 32: case 3:(p3)=(p1)^((p2)<<1)^(p2);(q3)=(q1)^((q2)<<1)^(q2);break;\
! 33: case 4:(p3)=(p1)^((p2)<<2);(q3)=(q1)^((q2)<<2);break;\
! 34: case 5:(p3)=(p1)^((p2)<<2)^(p2);(q3)=(q1)^((q2)<<2)^(q2);break;\
! 35: case 6:(p3)=(p1)^((p2)<<2)^((p2)<<1);(q3)=(q1)^((q2)<<2)^((q2)<<1);break;\
! 36: case 7:(p3)=(p1)^((p2)<<2)^((p2)<<1)^(p2);(q3)=(q1)^((q2)<<2)^((q2)<<1)^(q2);break;\
! 37: default:(p3)=(p1)^mulup2_bb((p2),(q));(q3)=(q1)^mulup2_bb((q2),(q));break;\
! 38: }\
! 39:
! 40: #if defined(P5)
! 41:
! 42: #define GF2M_MUL_1_HALF(high,low,a,b) (low)=NNgf2m_mul_1h(&(high),a,b);
! 43: #define GF2M_MUL_1(high,low,a,b) (low)=NNgf2m_mul_11(&(high),a,b);
! 44:
! 45: #else
! 46:
! 47: /*
! 48: * 32bit x 16bit -> 64bit (48bit); for b <= 0xffff
! 49: * a short-cut version of GF2M_MUL_1
! 50: */
! 51:
! 52: #define GF2M_MUL_1_HALF(high,low,a,b)\
! 53: {\
! 54: unsigned int _a,_b,_c,_d,_t,_s;\
! 55: unsigned int _ll,_cc;\
! 56: _a = (a);\
! 57: _b = (b);\
! 58: _c = _a&0xff00ff00; /* c = a4 0 a2 0 */\
! 59: _a ^= _c; /* a = 0 a3 0 a1 */\
! 60: _d = _b&0xff00ff00; /* d = 0 0 b2 0 */\
! 61: _b ^= _d; /* b = 0 0 0 b1 */\
! 62: _c |= (_d>>8); /* c = a4 00 a2 b2 */\
! 63: _b |= (_a<<8); /* b = a3 00 a1 b1 */\
! 64: _a = (_c>>16); /* a = a4 00 */\
! 65: _c &= 0xffff; /* c = a2 b2 */\
! 66: _d = (_b>>16); /* d = a3 00 */\
! 67: _b &= 0xffff; /* b = a1 b1 */\
! 68: _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\
! 69: _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\
! 70: _a ^= _c; _d^= _b;\
! 71: _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
! 72: _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
! 73: _cc ^= _ll;\
! 74: (low) = (unsigned int)(_ll^(_cc<<16));\
! 75: (high) = (unsigned int)(_cc>>16);\
! 76: }
! 77:
! 78: /*
! 79: * 32bit x 32bit -> 64bit
! 80: * This is an inline expansion of 4byte x 4byte Karatsuba multiplication
! 81: * Necessary indices for up2_bbtab are efficiently generated.
! 82: */
! 83:
! 84: #define GF2M_MUL_1(high,low,a,b)\
! 85: {\
! 86: unsigned int _a,_b,_c,_d,_t,_s;\
! 87: unsigned int _uu,_ll,_cc;\
! 88: _a = (a);\
! 89: _b = (b);\
! 90: _c = _a&0xff00ff00; /* _c = a4 0 a2 0 */\
! 91: _a ^= _c; /* a = 0 a3 0 a1 */\
! 92: _d = _b&0xff00ff00; /* _d = b4 0 b2 0 */\
! 93: _b ^= _d; /* b = 0 b3 0 b1 */\
! 94: _c |= (_d>>8); /* _c = a4 b4 a2 b2 */\
! 95: _b |= (_a<<8); /* b = a3 b3 a1 b1 */\
! 96: _a = (_c>>16); /* a = a4 b4 */\
! 97: _c &= 0xffff; /* _c = a2 b2 */\
! 98: _d = (_b>>16); /* _d = a3 b3 */\
! 99: _b &= 0xffff; /* b = a1 b1 */\
! 100: _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
! 101: _uu = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
! 102: _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\
! 103: _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\
! 104: _a ^= _c; _d^= _b;\
! 105: _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
! 106: _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
! 107: _cc ^= (_ll^_uu);\
! 108: (low) = (unsigned int)(_ll^(_cc<<16));\
! 109: (high) = (unsigned int)(_uu^(_cc>>16));\
! 110: }
! 111: #endif
! 112:
! 113: #define REMUP2_ONESTEP(a,q,s,n,w)\
! 114: if ( !s ) a[q] ^= w;\
! 115: else {\
! 116: a[q] ^= (w<<s);\
! 117: if ( q+1 <= n )\
! 118: a[q+1] ^= (w>>(32-s));\
! 119: }
! 120:
! 121: void ptoup2(n,nr)
! 122: P n;
! 123: UP2 *nr;
! 124: {
! 125: DCP dc;
! 126: UP2 r,s;
! 127: int d,w,i;
! 128:
! 129: if ( !n )
! 130: *nr = 0;
! 131: else if ( OID(n) == O_N )
! 132: if ( EVENN(NM((Q)n)) )
! 133: *nr = 0;
! 134: else
! 135: *nr = ONEUP2;
! 136: else {
! 137: d = UDEG(n); w = (d+1)/BSH + ((d+1)%BSH?1:0);
! 138: up2_var = VR(n);
! 139: W_NEWUP2(r,w);
! 140: for ( dc = DC(n); dc; dc = NEXT(dc) )
! 141: if ( !EVENN(NM((Q)COEF(dc))) ) {
! 142: i = QTOS(DEG(dc));
! 143: r->b[i/BSH] |= 1<<(i%BSH);
! 144: }
! 145: for ( i = w-1; i >= 0 && !r->b[i]; i-- );
! 146: r->w = i+1;
! 147: NEWUP2(s,r->w); *nr = s;
! 148: _copyup2(r,s);
! 149: }
! 150: }
! 151:
! 152: void ptoup2_sparse(n,nr)
! 153: P n;
! 154: UP2 *nr;
! 155: {
! 156: DCP dc;
! 157: UP2 s;
! 158: int d,i;
! 159:
! 160: if ( !n )
! 161: *nr = 0;
! 162: else if ( OID(n) == O_N )
! 163: if ( EVENN(NM((Q)n)) )
! 164: *nr = 0;
! 165: else {
! 166: NEWUP2(s,1); s->w = 1; s->b[0] = 0;
! 167: *nr = s;
! 168: }
! 169: else {
! 170: d = UDEG(n);
! 171: for ( dc = DC(n), i = 0; dc; dc = NEXT(dc) )
! 172: if ( !EVENN(NM((Q)COEF(dc))) )
! 173: i++;
! 174: NEWUP2(s,i); s->w = i; *nr = s;
! 175: for ( dc = DC(n), i = 0; dc; dc = NEXT(dc) )
! 176: if ( !EVENN(NM((Q)COEF(dc))) )
! 177: s->b[i++] = QTOS(DEG(dc));
! 178: }
! 179: }
! 180:
! 181: void up2top(n,nr)
! 182: UP2 n;
! 183: P *nr;
! 184: {
! 185: int i,d;
! 186: DCP dc0,dc;
! 187:
! 188: if ( !n )
! 189: *nr = 0;
! 190: else if ( !(d = degup2(n)) )
! 191: *nr = (P)ONE;
! 192: else {
! 193: for ( i = d, dc0 = 0; i >= 0; i-- )
! 194: if ( n->b[i/BSH] & (1<<(i%BSH)) ) {
! 195: NEXTDC(dc0,dc); STOQ(i,DEG(dc)); COEF(dc) = (P)ONE;
! 196: }
! 197: if ( !up2_var )
! 198: up2_var = CO->v;
! 199: MKP(up2_var,dc0,*nr);
! 200: }
! 201: }
! 202:
! 203: void up2tovect(n,nr)
! 204: UP2 n;
! 205: VECT *nr;
! 206: {
! 207: int i,d;
! 208: VECT v;
! 209:
! 210: if ( !n )
! 211: *nr = 0;
! 212: else {
! 213: d = degup2(n);
! 214: MKVECT(v,d+1); *nr = v;
! 215: for ( i = d; i >= 0; i-- )
! 216: if ( n->b[i/BSH] & (1<<(i%BSH)) )
! 217: v->body[i] = (pointer)ONE;
! 218: }
! 219: }
! 220:
! 221: void up2ton(p,n)
! 222: UP2 p;
! 223: Q *n;
! 224: {
! 225: N nm;
! 226: int w;
! 227:
! 228: if ( !p )
! 229: *n = 0;
! 230: else {
! 231: w = p->w;
! 232: nm = NALLOC(w);
! 233: nm->p = w;
! 234: bcopy(p->b,nm->b,w*sizeof(unsigned int));
! 235: NTOQ(nm,1,*n);
! 236: }
! 237: }
! 238:
! 239: void ntoup2(n,p)
! 240: Q n;
! 241: UP2 *p;
! 242: {
! 243: N nm;
! 244: UP2 t;
! 245: int w;
! 246:
! 247: if ( !n )
! 248: *p = 0;
! 249: else {
! 250: nm = NM(n); w = nm->p;
! 251: NEWUP2(t,w); *p = t; t->w = w;
! 252: bcopy(nm->b,t->b,w*sizeof(unsigned int));
! 253: }
! 254: }
! 255:
! 256: void gen_simpup2(p,m,r)
! 257: UP2 p;
! 258: GEN_UP2 m;
! 259: UP2 *r;
! 260: {
! 261: if ( lm_lazy || !m )
! 262: *r = p;
! 263: else if ( m->id == UP2_DENSE )
! 264: remup2(p,m->dense,r);
! 265: else
! 266: remup2_sparse(p,m->sparse,r);
! 267: if ( *r && !((*r)->w) )
! 268: *r = 0;
! 269: }
! 270:
! 271: void gen_simpup2_destructive(p,m)
! 272: UP2 p;
! 273: GEN_UP2 m;
! 274: {
! 275: UP2 t;
! 276:
! 277: if ( lm_lazy || !m )
! 278: return;
! 279: else if ( m->id == UP2_DENSE ) {
! 280: remup2(p,m->dense,&t);
! 281: _copyup2(t,p);
! 282: } else
! 283: remup2_sparse_destructive(p,m->sparse);
! 284: }
! 285:
! 286: void gen_invup2(p,m,r)
! 287: UP2 p;
! 288: GEN_UP2 m;
! 289: UP2 *r;
! 290: {
! 291: if ( !m )
! 292: error("gen_invup2 : invalid modulus");
! 293: else
! 294: invup2(p,m->dense,r);
! 295: }
! 296:
! 297: void gen_pwrmodup2(a,b,m,c)
! 298: UP2 a;
! 299: Q b;
! 300: GEN_UP2 m;
! 301: UP2 *c;
! 302: {
! 303: if ( !m )
! 304: pwrmodup2(a,b,0,c);
! 305: else if ( m->id == UP2_DENSE )
! 306: pwrmodup2(a,b,m->dense,c);
! 307: else
! 308: pwrmodup2_sparse(a,b,m->sparse,c);
! 309: }
! 310:
! 311: void simpup2(p,m,r)
! 312: UP2 p,m;
! 313: UP2 *r;
! 314: {
! 315: if ( !lm_lazy && m )
! 316: remup2(p,m,r);
! 317: else
! 318: *r = p;
! 319: }
! 320:
! 321: int degup2(a)
! 322: UP2 a;
! 323: {
! 324: unsigned int l,i,t;
! 325:
! 326: if ( !a || !a->w )
! 327: return -1;
! 328: else {
! 329: l = a->w; t = a->b[l-1];
! 330: for ( i = 0; t; t>>=1, i++);
! 331: return i+(l-1)*BSH-1;
! 332: }
! 333: }
! 334:
! 335: int degup2_sparse(a)
! 336: UP2 a;
! 337: {
! 338: if ( !a || !a->w )
! 339: return -1;
! 340: else
! 341: return a->b[0];
! 342: }
! 343:
! 344: int degup2_1(a)
! 345: unsigned int a;
! 346: {
! 347: int i;
! 348:
! 349: for ( i = 0; a; a>>=1, i++);
! 350: return i-1;
! 351: }
! 352:
! 353: void addup2(a,b,c)
! 354: UP2 a,b;
! 355: UP2 *c;
! 356: {
! 357: int i;
! 358: UP2 t;
! 359: int w,wa,wb;
! 360:
! 361: if ( !a )
! 362: *c = b;
! 363: else if ( !b )
! 364: *c = a;
! 365: else {
! 366: wa = a->w; wb = b->w;
! 367: if ( wa < wb ) {
! 368: t = a; a = b; b = t;
! 369: w = wa; wa = wb; wb = w;
! 370: }
! 371: NEWUP2(t,wa);
! 372: for ( i = 0; i < wb; i++ )
! 373: t->b[i] = a->b[i]^b->b[i];
! 374: for ( ; i < wa; i++ )
! 375: t->b[i] = a->b[i];
! 376: for ( i = wa-1; i >= 0 && !t->b[i]; i-- );
! 377: if ( i < 0 )
! 378: *c = 0;
! 379: else {
! 380: t->w = i+1;
! 381: *c = t;
! 382: }
! 383: }
! 384: }
! 385:
! 386: void subup2(a,b,c)
! 387: UP2 a,b;
! 388: UP2 *c;
! 389: {
! 390: addup2(a,b,c);
! 391: }
! 392:
! 393:
! 394: /*
! 395: s[w-1] ... s[0]
! 396: x m
! 397: ---------------------
! 398: t[w]t[w-1] ... t[0]
! 399: + r[w]r[w-1] ... r[0]
! 400: ---------------------
! 401: */
! 402:
! 403: INLINE void mulup2_n1(s,w,m,r)
! 404: unsigned int *s,*r;
! 405: int w;
! 406: unsigned int m;
! 407: {
! 408: int i;
! 409: unsigned int _u,_l,t;
! 410:
! 411: for ( i = 0; i < w; i++ )
! 412: if ( t = s[i] ) {
! 413: GF2M_MUL_1(_u,_l,t,m);
! 414: r[i] ^= _l; r[i+1] ^= _u;
! 415: }
! 416: }
! 417:
! 418: INLINE void mulup2_nh(s,w,m,r)
! 419: unsigned int *s,*r;
! 420: int w;
! 421: unsigned int m; /* 0 <= b <= 0xffff */
! 422: {
! 423: int i;
! 424: unsigned int u,l;
! 425:
! 426: for ( i = 0, r[0] = 0; i < w; i++ ) {
! 427: GF2M_MUL_1_HALF(u,l,s[i],m);
! 428: r[i] ^= l; r[i+1] = u;
! 429: }
! 430: }
! 431:
! 432: void _mulup2_1(a,b,c)
! 433: UP2 a,c;
! 434: unsigned int b;
! 435: {
! 436: int w;
! 437:
! 438: if ( !a || !a->w || !b )
! 439: c->w = 0;
! 440: else if ( b <= 0xffff )
! 441: _mulup2_h(a,b,c);
! 442: else {
! 443: w = a->w;
! 444: bzero((char *)c->b,(w+1)*sizeof(unsigned int));
! 445: mulup2_n1(a->b,w,b,c->b);
! 446: c->w = c->b[w] ? w+1 : w;
! 447: }
! 448: }
! 449:
! 450: void _mulup2_h(a,b,c)
! 451: UP2 a,c;
! 452: unsigned int b; /* 0 <= b <= 0xffff */
! 453: {
! 454: int w;
! 455:
! 456: if ( !a || !a->w || !b )
! 457: c->w = 0;
! 458: else {
! 459: w = a->w;
! 460: mulup2_nh(a->b,w,b,c->b);
! 461: c->w = c->b[w] ? w+1 : w;
! 462: }
! 463: }
! 464:
! 465: void mulup2(a,b,c)
! 466: UP2 a,b;
! 467: UP2 *c;
! 468: {
! 469: UP2 t;
! 470: int wa,wb,w;
! 471:
! 472: if ( !a || !b )
! 473: *c = 0;
! 474: else if ( a->w==1 && a->b[0]==1 ) {
! 475: NEWUP2(t,b->w); *c = t; _copyup2(b,t);
! 476: } else if ( b->w==1 && b->b[0]==1 ) {
! 477: NEWUP2(t,a->w); *c = t; _copyup2(a,t);
! 478: } else {
! 479: wa = a->w; wb = b->w;
! 480: if ( wa != wb ) {
! 481: if ( wb > wa ) {
! 482: t = a; a = b; b = t;
! 483: w = wa; wa = wb; wb = w;
! 484: }
! 485: W_NEWUP2(t,wa);
! 486: _copyup2(b,t);
! 487: bzero((char *)(t->b+wb),(wa-wb)*sizeof(unsigned int));
! 488: t->w = wa;
! 489: b = t;
! 490: }
! 491: NEWUP2(t,2*wa);
! 492: _kmulup2_(a->b,b->b,wa,t->b);
! 493: t->w = 2*wa; _adjup2(t);
! 494: *c = t;
! 495: }
! 496: }
! 497:
! 498: void _kmulup2_(a,b,w,c)
! 499: unsigned int *a,*b,*c;
! 500: int w;
! 501: {
! 502: switch ( w ) {
! 503: case 1: GF2M_MUL_1(c[1],c[0],*a,*b); break;
! 504: case 2: _mulup2_22(a,b,c); break;
! 505: case 3: _mulup2_33(a,b,c); break;
! 506: case 4: _mulup2_44(a,b,c); break;
! 507: case 5: _mulup2_55(a,b,c); break;
! 508: case 6: _mulup2_66(a,b,c); break;
! 509: default: _mulup2_nn(a,b,w,c); break;
! 510: }
! 511: }
! 512:
! 513: void _mulup2_nn(a,b,w,c)
! 514: unsigned int *a,*b,*c;
! 515: int w;
! 516: {
! 517: int wlow,whigh;
! 518: struct _oUP2 ablow,abhigh,alow,ahigh,blow,bhigh,aa,bb,mid,cmid;
! 519:
! 520: /* wlow >= whigh */
! 521: wlow = (w+1)/2; whigh = w-wlow;
! 522:
! 523: alow.w = wlow; alow.b = a;
! 524: ahigh.w = whigh; ahigh.b = a+wlow;
! 525:
! 526: blow.w = wlow; blow.b = b;
! 527: bhigh.w = whigh; bhigh.b = b+wlow;
! 528:
! 529: ablow.b = c;
! 530: abhigh.b = c+wlow*2;
! 531:
! 532: _kmulup2_(a,b,wlow,ablow.b); ablow.w = 2*wlow;
! 533: _kmulup2_(a+wlow,b+wlow,whigh,abhigh.b); abhigh.w = 2*whigh;
! 534:
! 535: W_NEW_UP2(aa,wlow);
! 536: _addup2_(&alow,&ahigh,&aa); aa.w = wlow;
! 537:
! 538: W_NEW_UP2(bb,wlow);
! 539: _addup2_(&blow,&bhigh,&bb); bb.w = wlow;
! 540:
! 541: W_NEW_UP2(mid,2*wlow);
! 542: _kmulup2_(aa.b,bb.b,wlow,mid.b); mid.w = 2*wlow;
! 543:
! 544: _addtoup2_(&ablow,&mid); _addtoup2_(&abhigh,&mid);
! 545:
! 546: cmid.w = 2*wlow; cmid.b = c+wlow;
! 547: _addtoup2_(&mid,&cmid);
! 548: }
! 549:
! 550: void _mulup2(a,b,c)
! 551: UP2 a,b,c;
! 552: {
! 553: int wa,wb,w;
! 554: int i;
! 555: unsigned int *ba,*bb;
! 556: unsigned int mul;
! 557:
! 558: if ( !a || !b || !a->w || !b->w ) {
! 559: c->w = 0; return;
! 560: }
! 561: wa = a->w; wb = b->w;
! 562: w = wa + wb;
! 563: bzero((char *)c->b,w*sizeof(unsigned int));
! 564: for ( i = 0, ba = a->b, bb = b->b; i < b->w; i++, bb++ )
! 565: if ( mul = *bb )
! 566: mulup2_n1(ba,wa,mul,c->b+i);
! 567: c->w = w;
! 568: _adjup2(c);
! 569: }
! 570:
! 571: void _mulup2_(a,b,c)
! 572: _UP2 a,b,c;
! 573: {
! 574: int wa,wb,w;
! 575: int i;
! 576: unsigned int *ba,*bb;
! 577: unsigned int mul;
! 578:
! 579: if ( !a || !b || !a->w || !b->w ) {
! 580: c->w = 0; return;
! 581: }
! 582: wa = a->w; wb = b->w;
! 583: w = wa + wb;
! 584: bzero((char *)c->b,w*sizeof(unsigned int));
! 585: for ( i = 0, ba = a->b, bb = b->b; i < b->w; i++, bb++ )
! 586: if ( mul = *bb )
! 587: mulup2_n1(ba,wa,mul,c->b+i);
! 588: c->w = w;
! 589: _adjup2_(c);
! 590: }
! 591:
! 592: void squareup2(n,nr)
! 593: UP2 n;
! 594: UP2 *nr;
! 595: {
! 596: int w,w2,i;
! 597: unsigned int s;
! 598: unsigned int *tb,*nb;
! 599: UP2 t;
! 600:
! 601: if ( !n )
! 602: *nr = 0;
! 603: else {
! 604: w = n->w; w2 = 2*w;
! 605: NEWUP2(t,w2); *nr = t;
! 606: tb = t->b;
! 607: nb = n->b;
! 608: for ( i = 0; i < w; i++ ) {
! 609: s = nb[i];
! 610: tb[2*i] = up2_sqtab[s&0xff]|(up2_sqtab[(s>>8)&0xff]<<16);
! 611: tb[2*i+1] = up2_sqtab[(s>>16)&0xff]|(up2_sqtab[s>>24]<<16);
! 612: }
! 613: t->w = tb[w2-1]?w2:w2-1;
! 614: }
! 615: }
! 616:
! 617: void _squareup2(n,nr)
! 618: UP2 n;
! 619: UP2 nr;
! 620: {
! 621: int w,w2,i;
! 622: unsigned int s;
! 623: unsigned int *tb,*nb;
! 624: UP2 t;
! 625:
! 626: if ( !n )
! 627: nr->w = 0;
! 628: else {
! 629: w = n->w; w2 = 2*w;
! 630: t = nr;
! 631: tb = t->b;
! 632: nb = n->b;
! 633: for ( i = 0; i < w; i++ ) {
! 634: s = nb[i];
! 635: tb[2*i] = up2_sqtab[s&0xff]|(up2_sqtab[(s>>8)&0xff]<<16);
! 636: tb[2*i+1] = up2_sqtab[(s>>16)&0xff]|(up2_sqtab[s>>24]<<16);
! 637: }
! 638: t->w = tb[w2-1]?w2:w2-1;
! 639: }
! 640: }
! 641:
! 642: void _adjup2(n)
! 643: UP2 n;
! 644: {
! 645: int i;
! 646: unsigned int *nb;
! 647:
! 648: nb = n->b;
! 649: for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );
! 650: i++;
! 651: n->w = i;
! 652: }
! 653:
! 654: void _adjup2_(n)
! 655: _UP2 n;
! 656: {
! 657: int i;
! 658: unsigned int *nb;
! 659:
! 660: nb = n->b;
! 661: for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );
! 662: i++;
! 663: n->w = i;
! 664: }
! 665:
! 666: void _addup2(a,b,c)
! 667: UP2 a,b,c;
! 668: {
! 669: int i,wa,wb,w;
! 670: UP2 t;
! 671: unsigned int *ab,*bb,*cb;
! 672:
! 673: if ( !a || !a->w ) {
! 674: _copyup2(b,c); return;
! 675: } else if ( !b || !b->w ) {
! 676: _copyup2(a,c); return;
! 677: }
! 678: wa = a->w; wb = b->w;
! 679: if ( wa < wb ) {
! 680: t = a; a = b; b = t;
! 681: w = wa; wa = wb; wb = w;
! 682: }
! 683: ab = a->b; bb = b->b; cb = c->b;
! 684: for ( i = 0; i < wb; i++ )
! 685: *cb++ = (*ab++)^(*bb++);
! 686: for ( ; i < wa; i++ )
! 687: *cb++ = *ab++;
! 688: c->w = wa;
! 689: _adjup2(c);
! 690: }
! 691:
! 692: /* a += b */
! 693:
! 694: void _addup2_destructive(a,b)
! 695: UP2 a,b;
! 696: {
! 697: int i,wa,wb;
! 698: unsigned int *ab,*bb;
! 699:
! 700: if ( !a->w ) {
! 701: _copyup2(b,a); return;
! 702: } else if ( !b->w )
! 703: return;
! 704: else {
! 705: wa = a->w; wb = b->w;
! 706: ab = a->b; bb = b->b;
! 707: if ( wa >= wb ) {
! 708: for ( i = 0; i < wb; i++ )
! 709: *ab++ ^= *bb++;
! 710: } else {
! 711: for ( i = 0; i < wa; i++ )
! 712: *ab++ ^= *bb++;
! 713: for ( ; i < wb; i++ )
! 714: *ab++ = *bb++;
! 715: a->w = wb;
! 716: }
! 717: _adjup2(a);
! 718: }
! 719: }
! 720:
! 721: void _addup2_(a,b,c)
! 722: _UP2 a,b,c;
! 723: {
! 724: int i,wa,wb,w;
! 725: _UP2 t;
! 726: unsigned int *ab,*bb,*cb;
! 727:
! 728: wa = a->w; wb = b->w;
! 729: if ( wa < wb ) {
! 730: t = a; a = b; b = t;
! 731: w = wa; wa = wb; wb = w;
! 732: }
! 733: ab = a->b; bb = b->b; cb = c->b;
! 734: for ( i = 0; i < wb; i++ )
! 735: *cb++ = (*ab++)^(*bb++);
! 736: for ( ; i < wa; i++ )
! 737: *cb++ = *ab++;
! 738: c->w = wa;
! 739: }
! 740:
! 741: void _addtoup2_(a,b)
! 742: _UP2 a,b;
! 743: {
! 744: int i,wa;
! 745: unsigned int *ab,*bb;
! 746:
! 747: wa = a->w;
! 748: ab = a->b; bb = b->b;
! 749: for ( i = 0; i < wa; i++ )
! 750: *bb++ ^= *ab++;
! 751: }
! 752:
! 753: /* 8bit x 8bit; also works if deg(a*b) < 32 */
! 754:
! 755: unsigned int mulup2_bb(a,b)
! 756: unsigned int a,b;
! 757: {
! 758: unsigned int t;
! 759:
! 760: if ( !a || !b )
! 761: return 0;
! 762: else {
! 763: t = 0;
! 764: while ( b ) {
! 765: if ( b & 1 )
! 766: t ^= a;
! 767: a <<= 1;
! 768: b >>= 1;
! 769: }
! 770: return t;
! 771: }
! 772: }
! 773:
! 774: void init_up2_tab()
! 775: {
! 776: unsigned int i,j;
! 777:
! 778: for ( i = 0; i < 256; i++ )
! 779: for ( j = 0; j <= i; j++ ) {
! 780: up2_bbtab[i<<8|j] = mulup2_bb(i,j);
! 781: up2_bbtab[j<<8|i] = up2_bbtab[i<<8|j];
! 782: }
! 783: for ( i = 0; i < 256; i++ )
! 784: up2_sqtab[i] = mulup2_bb(i,i);
! 785: }
! 786:
! 787: /*
! 788: compute q s.t. a*x^(BSH-1) = b*q+r
! 789: deg(b)=BSH-1, deg(a)<=BSH-1 => deg(q)<=BSH-1, deg(r)<=BSH-2
! 790: */
! 791:
! 792: INLINE unsigned int quoup2_11(a,b)
! 793: unsigned int a,b;
! 794: {
! 795: unsigned int q,i;
! 796:
! 797: q = 0;
! 798: for ( i = ((unsigned int)1)<<(BSH-1); i; i >>=1, b >>= 1 )
! 799: if ( i & a ) {
! 800: a ^= b;
! 801: q |= i;
! 802: }
! 803: return q;
! 804: }
! 805:
! 806: void divup2_1(a1,a2,e1,e2,qp,rp)
! 807: unsigned int a1,a2;
! 808: int e1,e2;
! 809: unsigned int *qp,*rp;
! 810: {
! 811: int i;
! 812: unsigned t,q;
! 813:
! 814: for ( i=e1-e2, t=1<<e1, a2<<=i, q = 0; i>=0; i--, t>>=1, a2>>=1 ) {
! 815: q<<=1;
! 816: if ( a1 & t ) {
! 817: q |= 1;
! 818: a1 ^= a2;
! 819: }
! 820: }
! 821: *qp = q; *rp = a1;
! 822: }
! 823:
! 824: void qrup2(a,b,q,r)
! 825: UP2 a,b;
! 826: UP2 *q,*r;
! 827: {
! 828: unsigned int msa,msb,t,q0;
! 829: int s,i,wq,wb;
! 830: UP2 as,bs,_q;
! 831:
! 832: if ( !b )
! 833: error("qrup2 : division by 0");
! 834: else if ( !a ) {
! 835: *q = 0; *r = 0; return;
! 836: } else if ( degup2(a) < degup2(b) ) {
! 837: *q = 0; *r = a; return;
! 838: }
! 839: msb = b->b[b->w-1];
! 840: for ( t = msb, s = 0; t; t>>=1, s++ );
! 841: s = BSH-s;
! 842:
! 843: W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);
! 844: _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);
! 845: as->b[as->w] = 0;
! 846:
! 847: wb = bs->w;
! 848: wq = as->w-wb;
! 849: NEWUP2(_q,wq+1); *q = _q;
! 850:
! 851: msb = bs->b[bs->w-1];
! 852: for ( i = wq; i >= 0; i-- ) {
! 853: msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));
! 854: if ( msa ) {
! 855: q0 = quoup2_11(msa,msb);
! 856: mulup2_n1(bs->b,wb,q0,as->b+i);
! 857: } else
! 858: q0 = 0;
! 859: _q->b[i] = q0;
! 860: }
! 861: for ( i = wq; i >= 0 && !_q->b[i]; i-- );
! 862: _q->w = i+1;
! 863:
! 864: for ( i = wb-1; i >= 0 && !as->b[i]; i-- );
! 865: if ( i < 0 )
! 866: *r = 0;
! 867: else {
! 868: as->w = i+1;
! 869: NEWUP2(*r,as->w);
! 870: _bshiftup2(as,s,*r);
! 871: }
! 872: }
! 873:
! 874: /* q->w >= a->w-b->w+2, r->w >= b->w */
! 875:
! 876: void _qrup2(a,b,q,r)
! 877: UP2 a,b;
! 878: UP2 q,r;
! 879: {
! 880: unsigned int msa,msb,t,q0;
! 881: int s,i,wq,wb;
! 882:
! 883: if ( !b || !b->w )
! 884: error("_qrup2 : division by 0");
! 885: else if ( !a || !a->w ) {
! 886: q->w = 0; r->w = 0; return;
! 887: } else if ( degup2(a) < degup2(b) ) {
! 888: q->w = 0; _copyup2(a,r); return;
! 889: }
! 890: msb = b->b[b->w-1];
! 891: for ( t = msb, s = BSH; t; t>>=1, s-- );
! 892:
! 893: _copyup2(a,r);
! 894: wb = b->w;
! 895: wq = r->w-wb;
! 896: r->b[r->w] = 0;
! 897:
! 898: msb = (b->b[b->w-1]<<s)|(b->w==1||!s?0:b->b[b->w-2]>>(BSH-s));
! 899: for ( i = wq; i >= 0; i-- ) {
! 900: msa = (s==BSH-1?0:r->b[wb+i]<<(s+1))|(wb+i-1<0?0:r->b[wb+i-1]>>(BSH-1-s));
! 901: if ( msa ) {
! 902: q0 = quoup2_11(msa,msb);
! 903: mulup2_n1(b->b,wb,q0,r->b+i);
! 904: } else
! 905: q0 = 0;
! 906: q->b[i] = q0;
! 907: }
! 908: for ( i = wq; i >= 0 && !q->b[i]; i-- );
! 909: q->w = i+1;
! 910:
! 911: for ( i = wb-1; i >= 0 && !r->b[i]; i-- );
! 912: if ( i < 0 )
! 913: r->w = 0;
! 914: else
! 915: r->w = i+1;
! 916: }
! 917:
! 918: void remup2(a,b,c)
! 919: UP2 a,b;
! 920: UP2 *c;
! 921: {
! 922: unsigned int msa,msb,t,q;
! 923: int s,i,wq,wb;
! 924: UP2 as,bs,r;
! 925:
! 926: if ( !b || !b->w )
! 927: error("remup2 : division by 0");
! 928: else if ( !a || !a->w ) {
! 929: *c = 0; return;
! 930: } else if ( degup2(a) < degup2(b) ) {
! 931: *c = a; return;
! 932: }
! 933: msb = b->b[b->w-1];
! 934: for ( t = msb, s = 0; t; t>>=1, s++ );
! 935: s = BSH-s;
! 936:
! 937: W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);
! 938: _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);
! 939: as->b[as->w] = 0;
! 940:
! 941: wb = bs->w;
! 942: wq = as->w-wb;
! 943: msb = bs->b[bs->w-1];
! 944: for ( i = wq; i >= 0; i-- ) {
! 945: msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));
! 946: if ( msa ) {
! 947: q = quoup2_11(msa,msb);
! 948: mulup2_n1(bs->b,wb,q,as->b+i);
! 949: }
! 950: }
! 951: for ( i = wb-1; i >= 0 && !as->b[i]; i-- );
! 952: if ( i < 0 )
! 953: *c = 0;
! 954: else {
! 955: as->w = i+1;
! 956: NEWUP2(r,as->w); *c = r;
! 957: _bshiftup2(as,s,r);
! 958: }
! 959: }
! 960:
! 961: void _remup2(a,b,c)
! 962: UP2 a,b,c;
! 963: {
! 964: unsigned int msa,msb,t,q;
! 965: int s,i,wq,wb;
! 966: UP2 as,bs;
! 967:
! 968: if ( !b || !b->w )
! 969: error("_remup2 : division by 0");
! 970: else if ( !a || !a->w ) {
! 971: c->w = 0; return;
! 972: } else if ( degup2(a) < degup2(b) ) {
! 973: _copyup2(a,c); return;
! 974: }
! 975: msb = b->b[b->w-1];
! 976: for ( t = msb, s = 0; t; t>>=1, s++ );
! 977: s = BSH-s;
! 978:
! 979: W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);
! 980: _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);
! 981: as->b[as->w] = 0;
! 982:
! 983: wb = bs->w;
! 984: wq = as->w-wb;
! 985: msb = bs->b[bs->w-1];
! 986: for ( i = wq; i >= 0; i-- ) {
! 987: msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));
! 988: if ( msa ) {
! 989: q = quoup2_11(msa,msb);
! 990: mulup2_n1(bs->b,wb,q,as->b+i);
! 991: }
! 992: }
! 993: for ( i = wb-1; i >= 0 && !as->b[i]; i-- );
! 994: if ( i < 0 )
! 995: c->w = 0;
! 996: else {
! 997: as->w = i+1;
! 998: _bshiftup2(as,s,c);
! 999: }
! 1000: }
! 1001:
! 1002: /* b = b->w|b->b[0]|b->b[1]|... -> b = x^b->[0]+x^b->[1]+... (b->w terms) */
! 1003:
! 1004: void remup2_sparse(a,b,c)
! 1005: UP2 a,b;
! 1006: UP2 *c;
! 1007: {
! 1008: int i,j,k,wa,wb,d,ds,db,dr,r;
! 1009: unsigned int ha,hb;
! 1010: unsigned int *tb;
! 1011:
! 1012: UP2 t,s;
! 1013:
! 1014: if ( !b )
! 1015: error("remup2_sparse : division by 0");
! 1016: else if ( !a ) {
! 1017: *c = 0; return;
! 1018: } else if ( degup2(a) < (db = degup2_sparse(b)) ) {
! 1019: *c = a; return;
! 1020: } else if ( b->w == (int)(b->b[0])+1 ) {
! 1021: NEWUP2(*c,a->w);
! 1022: _copyup2(a,*c);
! 1023: remup2_type1_destructive(*c,b->b[0]);
! 1024: if ( (*c)->w <= 0 )
! 1025: *c = 0;
! 1026: } else {
! 1027: W_NEWUP2(t,a->w); _copyup2(a,t);
! 1028: wa = a->w; wb = b->w; tb = t->b;
! 1029: for ( i = wa-1; (ds = BSH*i-db) >= 0; ) {
! 1030: if ( !(ha = tb[i]) )
! 1031: i--;
! 1032: else {
! 1033: tb[i] = 0;
! 1034: for ( j = 1; j < wb; j++ ) {
! 1035: d = ds+b->b[j];
! 1036: k = d/BSH; r = d%BSH;
! 1037: if ( !r )
! 1038: tb[k] ^= ha;
! 1039: else {
! 1040: tb[k] ^= (ha<<r);
! 1041: if ( k+1 <= i )
! 1042: tb[k+1] ^= (ha>>(BSH-r));
! 1043: }
! 1044: }
! 1045: if ( !tb[i] )
! 1046: i--;
! 1047: }
! 1048: }
! 1049: dr = db%BSH;
! 1050: hb = 1<<dr;
! 1051: i = db/BSH;
! 1052: while ( (ha=tb[i]) >= hb ) {
! 1053: ha >>= dr;
! 1054: t->b[i] &= (hb-1);
! 1055: for ( j = 1; j < wb; j++ ) {
! 1056: d = b->b[j];
! 1057: k = d/BSH; r = d%BSH;
! 1058: if ( !r )
! 1059: tb[k] ^= ha;
! 1060: else {
! 1061: tb[k] ^= (ha<<r);
! 1062: if ( k+1 <= i )
! 1063: tb[k+1] ^= (ha>>(BSH-r));
! 1064: }
! 1065: }
! 1066: }
! 1067: t->w = i+1;
! 1068: _adjup2(t);
! 1069: if ( t->w == 0 )
! 1070: *c = 0;
! 1071: else {
! 1072: NEWUP2(s,t->w); *c = s;
! 1073: _copyup2(t,s);
! 1074: }
! 1075: }
! 1076: }
! 1077:
! 1078: void remup2_sparse_destructive(a,b)
! 1079: UP2 a,b;
! 1080: {
! 1081: int i,j,k,wb,d,ds,db,dr,r;
! 1082: unsigned int ha,hb;
! 1083: unsigned int *ab,*bb;
! 1084:
! 1085: if ( !b )
! 1086: error("remup2_sparse_destructive : division by 0");
! 1087: else if ( !a || !a->w || ( degup2(a) < (db = degup2_sparse(b)) ) )
! 1088: return;
! 1089: else if ( b->w == 3 )
! 1090: remup2_3_destructive(a,b);
! 1091: else if ( b->w == 5 )
! 1092: remup2_5_destructive(a,b);
! 1093: else if ( b->w == (int)(b->b[0])+1 )
! 1094: remup2_type1_destructive(a,b->b[0]);
! 1095: else {
! 1096: wb = b->w; ab = a->b; bb = b->b;
! 1097: for ( i = a->w-1; (ds = (i<<5)-db) >= 0; ) {
! 1098: if ( !(ha = ab[i]) )
! 1099: i--;
! 1100: else {
! 1101: ab[i] = 0;
! 1102: for ( j = 1; j < wb; j++ ) {
! 1103: d = ds+bb[j]; k = d>>5; r = d&31;
! 1104: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1105: }
! 1106: if ( !ab[i] )
! 1107: i--;
! 1108: }
! 1109: }
! 1110: i = db>>5; dr = db&31;
! 1111: for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
! 1112: ha >>= dr;
! 1113: ab[i] &= hb;
! 1114: for ( j = 1; j < wb; j++ ) {
! 1115: d = bb[j]; k = d>>5; r = d&31;
! 1116: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1117: }
! 1118: }
! 1119: a->w = i+1;
! 1120: _adjup2(a);
! 1121: }
! 1122: }
! 1123:
! 1124: /* b = x^d+x^(d-1)+...+1 */
! 1125:
! 1126: void remup2_type1_destructive(a,d)
! 1127: UP2 a;
! 1128: int d;
! 1129: {
! 1130: int i,k,ds,db,dr;
! 1131: unsigned int ha,hb,r;
! 1132: unsigned int *ab;
! 1133:
! 1134: ab = a->b; db = d+1;
! 1135: for ( i = a->w-1; (ds = (i<<5)-db) >= 0; ) {
! 1136: if ( !(ha = ab[i]) ) i--;
! 1137: else {
! 1138: ab[i] = 0;
! 1139: k = ds>>5; r = ds&31;
! 1140: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1141: if ( !ab[i] ) i--;
! 1142: }
! 1143: }
! 1144: i = db>>5; dr = db&31;
! 1145: for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
! 1146: ha >>= dr; ab[i] &= hb; ab[0] ^= ha;
! 1147: }
! 1148: i = d>>5; hb = 1<<(d&31);
! 1149: if ( ab[i]&hb ) {
! 1150: ab[i] = (ab[i]^hb)^(hb-1);
! 1151: for ( k = i-1; k >= 0; k-- ) ab[k] ^= 0xffffffff;
! 1152: }
! 1153: a->w = i+1;
! 1154: _adjup2(a);
! 1155: }
! 1156:
! 1157: /* b = x^b->b[0]+x^b->b[1]+1 */
! 1158:
! 1159: void remup2_3_destructive(a,b)
! 1160: UP2 a,b;
! 1161: {
! 1162: int i,k,d,ds,db,db1,dr;
! 1163: unsigned int ha,hb,r;
! 1164: unsigned int *ab,*bb;
! 1165:
! 1166: ab = a->b; bb = b->b;
! 1167: db = bb[0]; db1 = bb[1];
! 1168: for ( i = a->w-1; (ds = (i<<5)-db) >= 0; ) {
! 1169: if ( !(ha = ab[i]) )
! 1170: i--;
! 1171: else {
! 1172: ab[i] = 0;
! 1173: d = ds+db1; k = d>>5; r = d&31;
! 1174: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1175: d = ds; k = d>>5; r = d&31;
! 1176: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1177: if ( !ab[i] )
! 1178: i--;
! 1179: }
! 1180: }
! 1181: i = db>>5; dr = db&31;
! 1182: k = db1>>5; r = db1&31;
! 1183: for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
! 1184: ha >>= dr;
! 1185: ab[i] &= hb;
! 1186: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1187: ab[0] ^= ha;
! 1188: }
! 1189: a->w = i+1;
! 1190: _adjup2(a);
! 1191: }
! 1192:
! 1193: /* b = x^b->b[0]+x^b->b[1]+x^b->b[2]+x^b->b[3]+1 */
! 1194:
! 1195: void remup2_5_destructive(a,b)
! 1196: UP2 a,b;
! 1197: {
! 1198: int i,d,ds,db,db1,db2,db3,dr;
! 1199: int k,k1,k2,k3;
! 1200: unsigned int ha,hb,r,r1,r2,r3;
! 1201: unsigned int *ab,*bb;
! 1202:
! 1203: ab = a->b; bb = b->b;
! 1204: db = bb[0]; db1 = bb[1]; db2 = bb[2]; db3 = bb[3];
! 1205: for ( i = a->w-1; (ds = (i<<5)-db) >= 0; ) {
! 1206: if ( !(ha = ab[i]) )
! 1207: i--;
! 1208: else {
! 1209: ab[i] = 0;
! 1210: d = ds+db1; k = d>>5; r = d&31;
! 1211: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1212: d = ds+db2; k = d>>5; r = d&31;
! 1213: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1214: d = ds+db3; k = d>>5; r = d&31;
! 1215: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1216: d = ds; k = d>>5; r = d&31;
! 1217: REMUP2_ONESTEP(ab,k,r,i,ha)
! 1218: if ( !ab[i] )
! 1219: i--;
! 1220: }
! 1221: }
! 1222: i = db>>5; dr = db&31;
! 1223: k1 = db1>>5; r1 = db1&31;
! 1224: k2 = db2>>5; r2 = db2&31;
! 1225: k3 = db3>>5; r3 = db3&31;
! 1226: for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
! 1227: ha >>= dr;
! 1228: ab[i] &= hb;
! 1229: REMUP2_ONESTEP(ab,k1,r1,i,ha)
! 1230: REMUP2_ONESTEP(ab,k2,r2,i,ha)
! 1231: REMUP2_ONESTEP(ab,k3,r3,i,ha)
! 1232: ab[0] ^= ha;
! 1233: }
! 1234: a->w = i+1;
! 1235: _adjup2(a);
! 1236: }
! 1237:
! 1238: void _invup2_1(f1,f2,a1,b1)
! 1239: unsigned int f1,f2,*a1,*b1;
! 1240: {
! 1241: unsigned int p1,p2,p3,q1,q2,q3,g1,g2,q,r;
! 1242: int d1,d2;
! 1243:
! 1244: p1 = 1; p2 = 0;
! 1245: q1 = 0; q2 = 1;
! 1246: g1 = f1; g2 = f2;
! 1247: d1 = degup2_1(g1); d2 = degup2_1(g2);
! 1248: while ( g2 ) {
! 1249: divup2_1(g1,g2,d1,d2,&q,&r);
! 1250: g1 = g2; g2 = r;
! 1251: GCD_COEF(q,p1,p2,p3,q1,q2,q3)
! 1252: p1 = p2; p2 = p3;
! 1253: q1 = q2; q2 = q3;
! 1254: d1 = d2; d2 = degup2_1(g2);
! 1255: }
! 1256: *a1 = p1; *b1 = q1;
! 1257: }
! 1258:
! 1259: void _gcdup2_1(f1,f2,gcd)
! 1260: unsigned int f1,f2,*gcd;
! 1261: {
! 1262: unsigned int g1,g2,q,r;
! 1263: int d1,d2;
! 1264:
! 1265: if ( !f1 )
! 1266: *gcd = f2;
! 1267: else if ( !f2 )
! 1268: *gcd = f1;
! 1269: else {
! 1270: g1 = f1; g2 = f2;
! 1271: d1 = degup2_1(g1); d2 = degup2_1(g2);
! 1272: while ( 1 ) {
! 1273: divup2_1(g1,g2,d1,d2,&q,&r);
! 1274: if ( !r ) {
! 1275: *gcd = g2;
! 1276: return;
! 1277: }
! 1278: g1 = g2; g2 = r;
! 1279: d1 = d2; d2 = degup2_1(g2);
! 1280: }
! 1281: }
! 1282: }
! 1283:
! 1284: static struct oEGT eg_a,eg_b;
! 1285:
! 1286: void up2_init_eg() {
! 1287: init_eg(&eg_a); init_eg(&eg_b);
! 1288: }
! 1289:
! 1290: void up2_show_eg() {
! 1291: print_eg("A",&eg_a);
! 1292: print_eg("B",&eg_b);
! 1293: printf("\n");
! 1294: }
! 1295:
! 1296: void invup2(a,m,inv)
! 1297: UP2 a,m;
! 1298: UP2 *inv;
! 1299: {
! 1300: int w,e1,e2,d1,d2;
! 1301: UP2 g1,g2,g3,a1,a2,a3,q,r,w1,w2,t;
! 1302: unsigned int p1,p2,p3,q0,q1,q2,q3,h1,h2,h3;
! 1303:
! 1304: if ( degup2(a) >= degup2(m) ) {
! 1305: remup2(a,m,&t); a = t;
! 1306: }
! 1307: w = m->w+2;
! 1308: W_NEWUP2(g1,w); W_NEWUP2(g2,w); W_NEWUP2(g3,w);
! 1309: W_NEWUP2(q,w); W_NEWUP2(r,w);
! 1310: W_NEWUP2(a1,w); W_NEWUP2(a2,w); W_NEWUP2(a3,w);
! 1311: W_NEWUP2(w1,w); W_NEWUP2(w2,w);
! 1312:
! 1313: a1->w = 1; a1->b[0] = 1;
! 1314: a2->w = 0;
! 1315: _copyup2(a,g1); _copyup2(m,g2);
! 1316:
! 1317: while ( 1 ) {
! 1318: d1 = degup2(g1); d2 = degup2(g2);
! 1319: if ( d1 < d2 ) {
! 1320: t = g1; g1 = g2; g2 = t;
! 1321: t = a1; a1 = a2; a2 = t;
! 1322: } else if ( d1 < 32 ) {
! 1323: _invup2_1(g1->b[0],g2->b[0],&p1,&p2);
! 1324: _mulup2_1(a1,p1,w1);
! 1325: _mulup2_1(a2,p2,w2);
! 1326: addup2(w1,w2,inv);
! 1327: return;
! 1328: } else if ( d1-d2 >= 16 ) {
! 1329: _qrup2(g1,g2,q,g3);
! 1330: t = g1; g1 = g2;
! 1331: if ( !g3->w ) {
! 1332: NEWUP2(t,a2->w); *inv = t;
! 1333: _copyup2(a2,t);
! 1334: return;
! 1335: } else {
! 1336: g2 = g3; g3 = t;
! 1337: }
! 1338: _mulup2(a2,q,w1); _addup2(a1,w1,a3);
! 1339: t = a1; a1 = a2; a2 = a3; a3 = t;
! 1340: } else {
! 1341: _bshiftup2(g1,d1-31,w1); h1 = w1->b[0];
! 1342: _bshiftup2(g2,d1-31,w2); h2 = w2->b[0];
! 1343: p1 = 1; p2 = 0;
! 1344: q1 = 0; q2 = 1;
! 1345: /* get_eg(&eg0); */
! 1346: e1 = degup2_1(h1); /* e1 must be 31 */
! 1347: while ( 1 ) {
! 1348: e2 = degup2_1(h2);
! 1349: if ( e2 <= 15 )
! 1350: break;
! 1351: divup2_1(h1,h2,e1,e2,&q0,&h3);
! 1352: GCD_COEF(q0,p1,p2,p3,q1,q2,q3)
! 1353: p1 = p2; p2 = p3;
! 1354: q1 = q2; q2 = q3;
! 1355: h1 = h2; h2 = h3;
! 1356: e1 = e2;
! 1357: }
! 1358: /* get_eg(&eg1); */
! 1359: _mulup2_1(g1,p1,w1); _mulup2_1(g2,q1,w2); _addup2(w1,w2,g3);
! 1360: _mulup2_1(g1,p2,w1); _mulup2_1(g2,q2,w2); _addup2(w1,w2,g2);
! 1361: t = g1; g1 = g3; g3 = t;
! 1362: _mulup2_1(a1,p1,w1); _mulup2_1(a2,q1,w2); _addup2(w1,w2,a3);
! 1363: _mulup2_1(a1,p2,w1); _mulup2_1(a2,q2,w2); _addup2(w1,w2,a2);
! 1364: t = a1; a1 = a3; a3 = t;
! 1365: /* get_eg(&eg2); */
! 1366: /* add_eg(&eg_a,&eg0,&eg1); */
! 1367: /* add_eg(&eg_b,&eg1,&eg2); */
! 1368: }
! 1369: }
! 1370: }
! 1371:
! 1372: void gcdup2(a,m,gcd)
! 1373: UP2 a,m;
! 1374: UP2 *gcd;
! 1375: {
! 1376: int w,e1,e2,d1,d2;
! 1377: UP2 g1,g2,g3,q,r,w1,w2,t;
! 1378: unsigned int p1,p2,p3,q0,q1,q2,q3,h1,h2,h3;
! 1379:
! 1380: if ( !a ) {
! 1381: *gcd = m; return;
! 1382: } else if ( !m ) {
! 1383: *gcd = a; return;
! 1384: }
! 1385: if ( degup2(a) >= degup2(m) ) {
! 1386: remup2(a,m,&t); a = t;
! 1387: if ( !a ) {
! 1388: *gcd = m; return;
! 1389: }
! 1390: }
! 1391: w = m->w+2;
! 1392: W_NEWUP2(g1,w); W_NEWUP2(g2,w); W_NEWUP2(g3,w);
! 1393: W_NEWUP2(q,w); W_NEWUP2(r,w);
! 1394: W_NEWUP2(w1,w); W_NEWUP2(w2,w);
! 1395:
! 1396: _copyup2(a,g1); _copyup2(m,g2);
! 1397:
! 1398: while ( 1 ) {
! 1399: d1 = degup2(g1); d2 = degup2(g2);
! 1400: if ( d1 < d2 ) {
! 1401: t = g1; g1 = g2; g2 = t;
! 1402: } else if ( d2 < 0 ) {
! 1403: NEWUP2(t,g1->w); *gcd = t;
! 1404: _copyup2(g1,t);
! 1405: return;
! 1406: } else if ( d1 < 32 ) {
! 1407: NEWUP2(t,1); t->w = 1; *gcd = t;
! 1408: _gcdup2_1(g1->b[0],g2->b[0],&t->b[0]);
! 1409: return;
! 1410: } else if ( d1-d2 >= 16 ) {
! 1411: _qrup2(g1,g2,q,g3);
! 1412: t = g1; g1 = g2;
! 1413: if ( !g3->w ) {
! 1414: NEWUP2(t,g2->w); *gcd = t;
! 1415: _copyup2(g2,t);
! 1416: return;
! 1417: } else {
! 1418: g2 = g3; g3 = t;
! 1419: }
! 1420: } else {
! 1421: _bshiftup2(g1,d1-31,w1); h1 = w1->b[0];
! 1422: _bshiftup2(g2,d1-31,w2); h2 = w2->b[0];
! 1423: p1 = 1; p2 = 0;
! 1424: q1 = 0; q2 = 1;
! 1425: e1 = degup2_1(h1); /* e1 must be 31 */
! 1426: while ( 1 ) {
! 1427: e2 = degup2_1(h2);
! 1428: if ( e2 <= 15 )
! 1429: break;
! 1430: divup2_1(h1,h2,e1,e2,&q0,&h3);
! 1431: GCD_COEF(q0,p1,p2,p3,q1,q2,q3)
! 1432: p1 = p2; p2 = p3;
! 1433: q1 = q2; q2 = q3;
! 1434: h1 = h2; h2 = h3;
! 1435: e1 = e2;
! 1436: }
! 1437: _mulup2_h(g1,p1,w1); _mulup2_h(g2,q1,w2); _addup2(w1,w2,g3);
! 1438: _mulup2_h(g1,p2,w1); _mulup2_h(g2,q2,w2); _addup2(w1,w2,g2);
! 1439: t = g1; g1 = g3; g3 = t;
! 1440: }
! 1441: }
! 1442: }
! 1443:
! 1444: void chsgnup2(a,c)
! 1445: UP2 a,*c;
! 1446: {
! 1447: *c = a;
! 1448: }
! 1449:
! 1450: void pwrmodup2(a,b,m,c)
! 1451: UP2 a;
! 1452: Q b;
! 1453: UP2 m;
! 1454: UP2 *c;
! 1455: {
! 1456: N n;
! 1457: UP2 y,t,t1;
! 1458: int k;
! 1459:
! 1460: if ( !b )
! 1461: *c = (UP2)ONEUP2;
! 1462: else if ( !a )
! 1463: *c = 0;
! 1464: else {
! 1465: y = (UP2)ONEUP2;
! 1466: n = NM(b);
! 1467: if ( !m )
! 1468: for ( k = n_bits(n)-1; k >= 0; k-- ) {
! 1469: squareup2(y,&t);
! 1470: if ( n->b[k/BSH] & (1<<(k%BSH)) )
! 1471: mulup2(t,a,&y);
! 1472: else
! 1473: y = t;
! 1474: }
! 1475: else
! 1476: for ( k = n_bits(n)-1; k >= 0; k-- ) {
! 1477: squareup2(y,&t);
! 1478: remup2(t,m,&t1); t = t1;
! 1479: if ( n->b[k/BSH] & (1<<(k%BSH)) ) {
! 1480: mulup2(t,a,&y);
! 1481: remup2(y,m,&t1); y = t1;
! 1482: }
! 1483: else
! 1484: y = t;
! 1485: }
! 1486: *c = y;
! 1487: }
! 1488: }
! 1489:
! 1490: void pwrmodup2_sparse(a,b,m,c)
! 1491: UP2 a;
! 1492: Q b;
! 1493: UP2 m;
! 1494: UP2 *c;
! 1495: {
! 1496: N n;
! 1497: UP2 y,t,t1;
! 1498: int k;
! 1499:
! 1500: if ( !b )
! 1501: *c = (UP2)ONEUP2;
! 1502: else if ( !a )
! 1503: *c = 0;
! 1504: else {
! 1505: y = (UP2)ONEUP2;
! 1506: n = NM(b);
! 1507: if ( !m )
! 1508: for ( k = n_bits(n)-1; k >= 0; k-- ) {
! 1509: squareup2(y,&t);
! 1510: if ( n->b[k/BSH] & (1<<(k%BSH)) )
! 1511: mulup2(t,a,&y);
! 1512: else
! 1513: y = t;
! 1514: }
! 1515: else
! 1516: for ( k = n_bits(n)-1; k >= 0; k-- ) {
! 1517: squareup2(y,&t);
! 1518: remup2_sparse(t,m,&t1); t = t1;
! 1519: if ( n->b[k/BSH] & (1<<(k%BSH)) ) {
! 1520: mulup2(t,a,&y);
! 1521: remup2_sparse(y,m,&t1); y = t1;
! 1522: }
! 1523: else
! 1524: y = t;
! 1525: }
! 1526: *c = y;
! 1527: }
! 1528: }
! 1529:
! 1530: int compup2(n1,n2)
! 1531: UP2 n1,n2;
! 1532: {
! 1533: int i;
! 1534: unsigned int *m1,*m2;
! 1535:
! 1536: if ( !n1 )
! 1537: if ( !n2 )
! 1538: return 0;
! 1539: else
! 1540: return -1;
! 1541: else if ( !n2 )
! 1542: return 1;
! 1543: else if ( n1->w > n2->w )
! 1544: return 1;
! 1545: else if ( n1->w < n2->w )
! 1546: return -1;
! 1547: else {
! 1548: for ( i = n1->w-1, m1 = n1->b+i, m2 = n2->b+i;
! 1549: i >= 0; i--, m1--, m2-- )
! 1550: if ( *m1 > *m2 )
! 1551: return 1;
! 1552: else if ( *m1 < *m2 )
! 1553: return -1;
! 1554: return 0;
! 1555: }
! 1556: }
! 1557:
! 1558: void _copyup2(n,r)
! 1559: UP2 n,r;
! 1560: {
! 1561: r->w = n->w;
! 1562: bcopy(n->b,r->b,n->w*sizeof(unsigned int));
! 1563: }
! 1564:
! 1565: void _bshiftup2(n,b,r)
! 1566: UP2 n;
! 1567: int b;
! 1568: UP2 r;
! 1569: {
! 1570: int w,l,nl,i,j;
! 1571: unsigned int msw;
! 1572: unsigned int *p,*pr;
! 1573:
! 1574: if ( b == 0 ) {
! 1575: _copyup2(n,r); return;
! 1576: }
! 1577: if ( b > 0 ) { /* >> */
! 1578: w = b / BSH; l = n->w-w;
! 1579: if ( l <= 0 ) {
! 1580: r->w = 0; return;
! 1581: }
! 1582: b %= BSH; p = n->b+w;
! 1583: if ( !b ) {
! 1584: r->w = l;
! 1585: bcopy(p,r->b,l*sizeof(unsigned int));
! 1586: return;
! 1587: }
! 1588: msw = p[l-1];
! 1589: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 1590: if ( b >= i ) {
! 1591: l--;
! 1592: if ( !l ) {
! 1593: r->w = 0; return;
! 1594: }
! 1595: r->w = l; pr = r->b;
! 1596: for ( j = 0; j < l; j++, p++ )
! 1597: *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 1598: } else {
! 1599: r->w = l; pr = r->b;
! 1600: for ( j = 1; j < l; j++, p++ )
! 1601: *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 1602: *pr = *p>>b;
! 1603: }
! 1604: } else { /* << */
! 1605: b = -b;
! 1606: w = b / BSH; b %= BSH; l = n->w; p = n->b;
! 1607: if ( !b ) {
! 1608: nl = l+w; r->w = nl;
! 1609: bzero((char *)r->b,w*sizeof(unsigned int));
! 1610: bcopy(p,r->b+w,l*sizeof(unsigned int));
! 1611: return;
! 1612: }
! 1613: msw = p[l-1];
! 1614: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 1615: if ( b + i > BSH ) {
! 1616: nl = l+w+1;
! 1617: r->w = nl; pr = r->b+w;
! 1618: bzero((char *)r->b,w*sizeof(unsigned int));
! 1619: *pr++ = *p++<<b;
! 1620: for ( j = 1; j < l; j++, p++ )
! 1621: *pr++ = (*p<<b)|(*(p-1)>>(BSH-b));
! 1622: *pr = *(p-1)>>(BSH-b);
! 1623: } else {
! 1624: nl = l+w;
! 1625: r->w = nl; pr = r->b+w;
! 1626: bzero((char *)r->b,w*sizeof(unsigned int));
! 1627: *pr++ = *p++<<b;
! 1628: for ( j = 1; j < l; j++, p++ )
! 1629: *pr++ = (*p<<b)|(*(p-1)>>(BSH-b));
! 1630: }
! 1631: }
! 1632: }
! 1633:
! 1634: void _bshiftup2_destructive(n,b)
! 1635: UP2 n;
! 1636: int b;
! 1637: {
! 1638: int w,l,nl,i,j;
! 1639: unsigned int msw;
! 1640: unsigned int *p,*pr;
! 1641:
! 1642: if ( b == 0 || !n->w )
! 1643: return;
! 1644: if ( b > 0 ) { /* >> */
! 1645: w = b / BSH; l = n->w-w;
! 1646: if ( l <= 0 ) {
! 1647: n->w = 0; return;
! 1648: }
! 1649: b %= BSH; p = n->b+w;
! 1650: if ( !b ) {
! 1651: n->w = l;
! 1652: bcopy(p,n->b,l*sizeof(unsigned int));
! 1653: return;
! 1654: }
! 1655: msw = p[l-1];
! 1656: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 1657: if ( b >= i ) {
! 1658: l--;
! 1659: if ( !l ) {
! 1660: n->w = 0; return;
! 1661: }
! 1662: n->w = l; pr = n->b;
! 1663: for ( j = 0; j < l; j++, p++ )
! 1664: *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 1665: } else {
! 1666: n->w = l; pr = n->b;
! 1667: for ( j = 1; j < l; j++, p++ )
! 1668: *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
! 1669: *pr = *p>>b;
! 1670: }
! 1671: } else { /* << */
! 1672: b = -b;
! 1673: w = b / BSH; b %= BSH; l = n->w; p = n->b;
! 1674: if ( !b ) {
! 1675: nl = l+w; n->w = nl;
! 1676: bcopy(p,n->b+w,l*sizeof(unsigned int));
! 1677: bzero((char *)n->b,w*sizeof(unsigned int));
! 1678: return;
! 1679: }
! 1680: msw = p[l-1];
! 1681: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
! 1682: if ( b + i > BSH ) {
! 1683: nl = l+w+1;
! 1684: n->w = nl; p = n->b+l-1; pr = n->b+w+l;
! 1685: *pr-- = *p>>(BSH-b);
! 1686: for ( j = l-1; j >= 1; j--, p-- )
! 1687: *pr-- = (*p<<b)|(*(p-1)>>(BSH-b));
! 1688: *pr = *p<<b;
! 1689: bzero((char *)n->b,w*sizeof(unsigned int));
! 1690: } else {
! 1691: nl = l+w;
! 1692: n->w = nl; p = n->b+l-1; pr = n->b+w+l-1;
! 1693: for ( j = l-1; j >= 1; j--, p-- )
! 1694: *pr-- = (*p<<b)|(*(p-1)>>(BSH-b));
! 1695: *pr = *p<<b;
! 1696: bzero((char *)n->b,w*sizeof(unsigned int));
! 1697: }
! 1698: }
! 1699: }
! 1700:
! 1701: void diffup2(f,r)
! 1702: UP2 f;
! 1703: UP2 *r;
! 1704: {
! 1705: int d,i,w;
! 1706: UP2 t;
! 1707:
! 1708: d = degup2(f);
! 1709: w = f->w;
! 1710: NEWUP2(t,w);
! 1711: _bshiftup2(f,1,t);
! 1712: for ( i = 0; i < d; i++ )
! 1713: if ( i%2 )
! 1714: t->b[i/BSH] &= ~(1<<(i%BSH));
! 1715: for ( i = w-1; i >= 0 && !t->b[i]; i-- );
! 1716: if ( i < 0 )
! 1717: *r = 0;
! 1718: else {
! 1719: *r = t;
! 1720: t->w = i+1;
! 1721: }
! 1722: }
! 1723:
! 1724: int sqfrcheckup2(f)
! 1725: UP2 f;
! 1726: {
! 1727: UP2 df,g;
! 1728:
! 1729: diffup2(f,&df);
! 1730: gcdup2(f,df,&g);
! 1731: if ( degup2(g) >= 1 )
! 1732: return 0;
! 1733: else
! 1734: return 1;
! 1735: }
! 1736:
! 1737: int irredcheckup2(f)
! 1738: UP2 f;
! 1739: {
! 1740: int n,w,i,j,k,hcol;
! 1741: unsigned int hbit;
! 1742: unsigned int *u;
! 1743: unsigned int **mat;
! 1744: UP2 p,t;
! 1745:
! 1746: if ( !sqfrcheckup2(f) )
! 1747: return 0;
! 1748: n = degup2(f);
! 1749: w = n/BSH + (n%BSH?1:0);
! 1750: mat = (unsigned int **)almat(n,w);
! 1751: W_NEWUP2(p,w+1);
! 1752: W_NEWUP2(t,w+1);
! 1753: p->w = 1; p->b[0] = 1; /* p = 1 mod f */
! 1754: for ( i = 1; i < n; i++ ) {
! 1755: _bshiftup2(p,-2,t); _remup2(t,f,p);
! 1756: _copy_up2bits(p,mat,i);
! 1757:
! 1758: }
! 1759: /* _print_frobmat(mat,n,w); */
! 1760: for ( j = 1; j < n; j++ ) {
! 1761: hcol = j/BSH; hbit = 1<<(j%BSH);
! 1762: for ( i = j-1; i < n; i++ )
! 1763: if ( mat[i][hcol] & hbit )
! 1764: break;
! 1765: if ( i == n )
! 1766: return 0;
! 1767: if ( i != j-1 ) {
! 1768: u = mat[i]; mat[i] = mat[j-1]; mat[j-1] = u;
! 1769: }
! 1770: for ( i = j; i < n; i++ )
! 1771: if ( mat[i][hcol] & hbit )
! 1772: for ( k = hcol; k < w; k++ )
! 1773: mat[i][k] ^= mat[j-1][k];
! 1774: }
! 1775: return 1;
! 1776: }
! 1777:
! 1778: int irredcheck_dddup2(f)
! 1779: UP2 f;
! 1780: {
! 1781: UP2 x,u,t,s,gcd;
! 1782: int n,i;
! 1783:
! 1784: if ( !sqfrcheckup2(f) )
! 1785: return -1;
! 1786: n = degup2(f);
! 1787:
! 1788: W_NEWUP2(u,1); u->w = 1; u->b[0] = 2; /* u = x mod f */
! 1789: x = u;
! 1790: for ( i = 1; 2*i <= n; i++ ) {
! 1791: squareup2(u,&t); remup2(t,f,&u); addup2(u,x,&s);
! 1792: gcdup2(f,s,&gcd);
! 1793: if ( degup2(gcd) > 0 )
! 1794: return 0;
! 1795: }
! 1796: return 1;
! 1797: }
! 1798:
! 1799: void _copy_up2bits(p,mat,pos)
! 1800: UP2 p;
! 1801: unsigned int **mat;
! 1802: int pos;
! 1803: {
! 1804: int d,col,j,jcol,jsh;
! 1805: unsigned int bit;
! 1806:
! 1807: d = degup2(p);
! 1808: col = pos/BSH; bit = 1<<(pos%BSH);
! 1809:
! 1810: for ( j = 0; j <= d; j++ ) {
! 1811: jcol = j/BSH; jsh = j%BSH;
! 1812: if ( p->b[jcol]&(1<<jsh) )
! 1813: mat[j][col] |= bit;
! 1814: }
! 1815: mat[pos][col] ^= bit;
! 1816: }
! 1817:
! 1818: int compute_multiplication_matrix(p0,mp)
! 1819: P p0;
! 1820: GF2MAT *mp;
! 1821: {
! 1822: UP2 p;
! 1823: int n,w,i,j,k,l;
! 1824: unsigned int **a,**b,**c,**d,**t,**m;
! 1825: GF2MAT am,bm;
! 1826:
! 1827: ptoup2(p0,&p);
! 1828: n = degup2(p);
! 1829: w = p->w;
! 1830: if ( !compute_representation_conversion_matrix(p0,&am,&bm) )
! 1831: return 0;
! 1832: a = am->body; b = bm->body;
! 1833: /* _print_frobmat(a,n,w); printf("\n"); */
! 1834: /* _print_frobmat(b,n,w); printf("\n"); */
! 1835: c = (unsigned int **)almat(n,w);
! 1836: for ( i = 0; i < n-1; i++ ) {
! 1837: j = i+1;
! 1838: c[i][j/BSH] |= 1<<(j%BSH);
! 1839: }
! 1840: bcopy(p->b,c[n-1],p->w*sizeof(unsigned int));
! 1841:
! 1842: /* _print_frobmat(b,n,w); printf("\n"); */
! 1843: t = (unsigned int **)almat(n,w);
! 1844: mulgf2mat(n,a,c,t);
! 1845: d = (unsigned int **)almat(n,w);
! 1846: mulgf2mat(n,t,b,d);
! 1847: /* _print_frobmat(d,n,w); printf("\n"); */
! 1848:
! 1849: m = (unsigned int **)almat(n,w);
! 1850: for ( i = 0; i < n; i++ )
! 1851: for ( j = 0; j < n; j++ ) {
! 1852: k = (n-1-(j-i))%n; l = (n-1+i)%n;
! 1853: if ( d[k][l/BSH] & 1<<(l%BSH) )
! 1854: m[n-1-i][(n-1-j)/BSH] |= 1<<((n-1-j)%BSH);
! 1855: }
! 1856: /* _print_frobmat(m,n,w); printf("\n"); */
! 1857: TOGF2MAT(n,n,m,*mp);
! 1858: return 1;
! 1859: }
! 1860:
! 1861: #define GF2N_PBTOPB 0
! 1862: #define GF2N_NBTOPB 1
! 1863:
! 1864: void compute_change_of_basis_matrix_with_root(P,P,int,GF2N,GF2MAT *,GF2MAT *);
! 1865:
! 1866: /*
! 1867: * if 'to' = GF2N_NBTOPB then p0 must be a normal poly.
! 1868: * rep0 x m01 -> rep1, rep1 x m10 -> rep0
! 1869: */
! 1870:
! 1871: void compute_change_of_basis_matrix(p0,p1,to,m01,m10)
! 1872: P p0,p1;
! 1873: int to;
! 1874: GF2MAT *m01,*m10;
! 1875: {
! 1876: UP2 up0;
! 1877: int n,w;
! 1878: unsigned int **p01,**p10;
! 1879: GF2N root;
! 1880:
! 1881: setmod_gf2n(p1);
! 1882:
! 1883: ptoup2(p0,&up0);
! 1884: n = degup2(up0); w = up0->w;
! 1885: find_root_up2(up0,&root);
! 1886: compute_change_of_basis_matrix_with_root(p0,p1,to,root,m01,m10);
! 1887: }
! 1888:
! 1889: void compute_change_of_basis_matrix_with_root(p0,p1,to,root,m01,m10)
! 1890: P p0,p1;
! 1891: int to;
! 1892: GF2N root;
! 1893: GF2MAT *m01,*m10;
! 1894: {
! 1895: UP2 up0,t,u,s;
! 1896: int n,w,i;
! 1897: unsigned int **a,**b,**g,**h;
! 1898: unsigned int **p01,**p10;
! 1899: P tmp;
! 1900:
! 1901: setmod_gf2n(p1);
! 1902:
! 1903: ptoup2(p0,&up0);
! 1904: n = degup2(up0); w = up0->w;
! 1905: substp(CO,p0,p0->v,(P)root,&tmp);
! 1906: if ( tmp )
! 1907: error("error in find_root_up2");
! 1908: u = root->body;
! 1909:
! 1910: p01 = (unsigned int **)almat(n,w);
! 1911: p10 = (unsigned int **)almat(n,w);
! 1912: if ( to == GF2N_PBTOPB ) {
! 1913: t = ONEUP2;
! 1914: bcopy(t->b,p01[0],t->w*sizeof(unsigned int));
! 1915: for ( i = 1; i < n; i++ ) {
! 1916: mulup2(t,u,&s); gen_simpup2_destructive(s,current_mod_gf2n); t = s;
! 1917: bcopy(t->b,p01[i],t->w*sizeof(unsigned int));
! 1918: }
! 1919: } else {
! 1920: t = u;
! 1921: bcopy(t->b,p01[n-1],t->w*sizeof(unsigned int));
! 1922: for ( i = 1; i < n; i++ ) {
! 1923: squareup2(t,&s); gen_simpup2_destructive(s,current_mod_gf2n); t = s;
! 1924: bcopy(t->b,p01[n-1-i],t->w*sizeof(unsigned int));
! 1925: }
! 1926: }
! 1927: if ( !invgf2mat(n,p01,p10) )
! 1928: error("compute_change_of_basis_matrix : something is wrong");
! 1929: TOGF2MAT(n,n,p01,*m01);
! 1930: TOGF2MAT(n,n,p10,*m10);
! 1931: }
! 1932:
! 1933: /*
! 1934: *
! 1935: * computation of representation conversion matrix
! 1936: * repn x np -> repp, repp x pn -> repn
! 1937: *
! 1938: */
! 1939:
! 1940: int compute_representation_conversion_matrix(p0,np,pn)
! 1941: P p0;
! 1942: GF2MAT *np,*pn;
! 1943: {
! 1944: UP2 x,s;
! 1945: int n,w,i;
! 1946: unsigned int **ntop,**pton;
! 1947:
! 1948: setmod_gf2n(p0);
! 1949:
! 1950: n = UDEG(p0); w = n/BSH + (n%BSH?1:0);
! 1951:
! 1952: /* x = alpha */
! 1953: NEWUP2(x,1); x->w = 1; x->b[0]=2;
! 1954: gen_simpup2_destructive(x,current_mod_gf2n);
! 1955:
! 1956: ntop = (unsigned int **)almat(n,w);
! 1957: pton = (unsigned int **)almat(n,w);
! 1958:
! 1959: bcopy(x->b,ntop[n-1],x->w*sizeof(unsigned int));
! 1960: for ( i = 1; i < n; i++ ) {
! 1961: squareup2(x,&s); gen_simpup2_destructive(s,current_mod_gf2n); x = s;
! 1962: bcopy(x->b,ntop[n-1-i],x->w*sizeof(unsigned int));
! 1963: }
! 1964: if ( !invgf2mat(n,ntop,pton) )
! 1965: return 0;
! 1966: TOGF2MAT(n,n,ntop,*np);
! 1967: TOGF2MAT(n,n,pton,*pn);
! 1968: return 1;
! 1969: }
! 1970:
! 1971: void mul_nb(mat,a,b,c)
! 1972: GF2MAT mat;
! 1973: unsigned int *a,*b,*c;
! 1974: {
! 1975: int n,w,i;
! 1976: unsigned int *wa,*wb,*t;
! 1977: unsigned int **m;
! 1978:
! 1979: n = mat->row;
! 1980: m = mat->body;
! 1981: w = n/BSH + (n%BSH?1:0);
! 1982: wa = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
! 1983: wb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
! 1984: t = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
! 1985: bcopy(a,wa,w*sizeof(unsigned int));
! 1986: bcopy(b,wb,w*sizeof(unsigned int));
! 1987: bzero((char *)c,w*sizeof(unsigned int));
! 1988: for ( i = n-1; i >= 0; i-- ) {
! 1989: mulgf2vectmat(n,wa,m,t);
! 1990: if ( mulgf2vectvect(n,t,wb) )
! 1991: c[i/BSH] |= 1<<(i%BSH);
! 1992: leftshift(wa,n);
! 1993: leftshift(wb,n);
! 1994: }
! 1995: }
! 1996:
! 1997:
! 1998: /*
! 1999: a=(c0,c1,...,c{n-1})->(c1,c2,...,c{n-1},c0)
! 2000: */
! 2001:
! 2002: void leftshift(a,n)
! 2003: unsigned int *a;
! 2004: int n;
! 2005: {
! 2006: int r,w,i;
! 2007: unsigned int msb;
! 2008:
! 2009: r = n%BSH;
! 2010: w = n/BSH + (r?1:0);
! 2011: msb = a[(n-1)/BSH]&(1<<((n-1)%BSH));
! 2012: for ( i = w-1; i > 0; i-- )
! 2013: a[i] = (a[i]<<1)|(a[i-1]>>(BSH-1));
! 2014: a[0] = (a[0]<<1)|(msb?1:0);
! 2015: if ( r )
! 2016: a[w-1] &= (1<<r)-1;
! 2017: }
! 2018:
! 2019: void mat_to_gf2mat(a,b)
! 2020: MAT a;
! 2021: unsigned int ***b;
! 2022: {
! 2023: int n,w,i,j;
! 2024: unsigned int **m;
! 2025:
! 2026: n = a->row;
! 2027: w = n/BSH + (n%BSH?1:0);
! 2028: *b = m = (unsigned int **)almat(n,w);
! 2029: for ( i = 0; i < n; i++ )
! 2030: for ( j = 0; j < n; j++ )
! 2031: if ( a->body[i][j] )
! 2032: m[i][j/BSH] |= 1<<(j%BSH);
! 2033: }
! 2034:
! 2035: void gf2mat_to_mat(a,n,b)
! 2036: unsigned int **a;
! 2037: MAT *b;
! 2038: {
! 2039: int w,i,j;
! 2040: MAT m;
! 2041:
! 2042: MKMAT(m,n,n); *b = m;
! 2043: w = n/BSH + (n%BSH?1:0);
! 2044: for ( i = 0; i < n; i++ )
! 2045: for ( j = 0; j < n; j++ )
! 2046: if ( a[i][j/BSH] & 1<<(j%BSH) )
! 2047: m->body[i][j] = (pointer)ONE;
! 2048: }
! 2049:
! 2050: void mulgf2mat(n,a,b,c)
! 2051: int n;
! 2052: unsigned int **a,**b,**c;
! 2053: {
! 2054: int i,j,k,w;
! 2055:
! 2056: w = n/BSH + (n%BSH?1:0);
! 2057: for ( i = 0; i < n; i++ ) {
! 2058: bzero((char *)c[i],w*sizeof(unsigned int));
! 2059: for ( j = 0; j < n; j++ )
! 2060: if ( a[i][j/BSH] & 1<<(j%BSH) )
! 2061: for ( k = 0; k < w; k++ )
! 2062: c[i][k] ^= b[j][k];
! 2063: }
! 2064: }
! 2065:
! 2066: /* c = a*b; where a, c are row vectors */
! 2067:
! 2068: void mulgf2vectmat(n,a,b,c)
! 2069: int n;
! 2070: unsigned int *a;
! 2071: unsigned int **b;
! 2072: unsigned int *c;
! 2073: {
! 2074: int j,k,w;
! 2075:
! 2076: w = n/BSH + (n%BSH?1:0);
! 2077: bzero((char *)c,w*sizeof(unsigned int));
! 2078: for ( j = 0; j < n; j++ )
! 2079: if ( a[j/BSH] & 1<<(j%BSH) )
! 2080: for ( k = 0; k < w; k++ )
! 2081: c[k] ^= b[j][k];
! 2082: }
! 2083:
! 2084: int mulgf2vectvect(n,a,b)
! 2085: int n;
! 2086: unsigned int *a,*b;
! 2087: {
! 2088: unsigned int t,r;
! 2089: int i,w;
! 2090:
! 2091: w = n/BSH + (n%BSH?1:0);
! 2092: for ( i = 0, r = 0; i < w; i++ )
! 2093: for ( t = a[i]&b[i]; t; t >>=1 )
! 2094: r ^= (t&1);
! 2095: return r;
! 2096: }
! 2097:
! 2098: int invgf2mat(n,a,b)
! 2099: int n;
! 2100: unsigned int **a,**b;
! 2101: {
! 2102: int i,j,k,hcol,hbit,w;
! 2103: unsigned int *u;
! 2104: unsigned int **s;
! 2105:
! 2106: w = n/BSH + (n%BSH?1:0);
! 2107: s = (unsigned int **)almat(n,w);
! 2108: for ( i = 0; i < n; i++ )
! 2109: bcopy(a[i],s[i],w*sizeof(unsigned int));
! 2110:
! 2111: for ( i = 0; i < n; i++ )
! 2112: b[i][i/BSH] |= 1<<(i%BSH);
! 2113:
! 2114: for ( j = 0; j < n; j++ ) {
! 2115: hcol = j/BSH; hbit = 1<<(j%BSH);
! 2116: for ( i = j; i < n; i++ )
! 2117: if ( s[i][hcol] & hbit )
! 2118: break;
! 2119: if ( i == n )
! 2120: return 0;
! 2121: if ( i != j ) {
! 2122: u = s[i]; s[i] = s[j]; s[j] = u;
! 2123: u = b[i]; b[i] = b[j]; b[j] = u;
! 2124: }
! 2125: for ( i = 0; i < n; i++ ) {
! 2126: if ( i == j )
! 2127: continue;
! 2128: if ( s[i][hcol] & hbit ) {
! 2129: for ( k = hcol; k < w; k++ )
! 2130: s[i][k] ^= s[j][k];
! 2131: for ( k = 0; k < w; k++ )
! 2132: b[i][k] ^= b[j][k];
! 2133: }
! 2134: }
! 2135: }
! 2136: return 1;
! 2137: }
! 2138:
! 2139: INLINE void _mulup2_11(a1,a2,ar)
! 2140: unsigned int a1,a2;
! 2141: unsigned int *ar;
! 2142: {
! 2143: GF2M_MUL_1(ar[1],ar[0],a1,a2);
! 2144: }
! 2145:
! 2146: void _mulup2_22(a1,a2,ar)
! 2147: unsigned int *a1,*a2,*ar;
! 2148: {
! 2149: unsigned int m[2];
! 2150:
! 2151: _mulup2_11(*a1,*a2,ar);
! 2152: _mulup2_11(*(a1+1),*(a2+1),ar+2);
! 2153: _mulup2_11(a1[0]^a1[1],a2[0]^a2[1],m);
! 2154: m[0] ^= ar[0]^ar[2]; m[1] ^= ar[1]^ar[3];
! 2155:
! 2156: ar[1] ^= m[0]; ar[2] ^= m[1];
! 2157: }
! 2158:
! 2159: #if 0
! 2160: void _mulup2_33(a1,a2,ar)
! 2161: unsigned int *a1,*a2,*ar;
! 2162: {
! 2163: unsigned int m[4];
! 2164: unsigned int c1[2],c2[2];
! 2165:
! 2166: c1[0] = a1[2]^a1[0]; c1[1] = a1[1];
! 2167: c2[0] = a2[2]^a2[0]; c2[1] = a2[1];
! 2168:
! 2169: _mulup2_22(a1,a2,ar);
! 2170: _mulup2_11(*(a1+2),*(a2+2),ar+4);
! 2171: _mulup2_22(c1,c2,m);
! 2172:
! 2173: m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];
! 2174: m[2] ^= ar[2]; m[3] ^= ar[3];
! 2175:
! 2176: ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
! 2177: }
! 2178: #else
! 2179: /* (ar[5]...ar[0]) = (a1[2]a1[1]a1[0])*(a2[2]a2[1]a2[0]) */
! 2180:
! 2181: void _mulup2_33(a1,a2,ar)
! 2182: unsigned int *a1,*a2,*ar;
! 2183: {
! 2184: unsigned int m[4];
! 2185: unsigned int c[2];
! 2186: unsigned int t,s;
! 2187:
! 2188: /* (ar[1],ar[0]) = a1[0]*a2[0] */
! 2189: _mulup2_11(a1[0],a2[0],ar);
! 2190:
! 2191: /* (ar[3],ar[2]) = hi(m) = a1[1]*a1[1] */
! 2192: _mulup2_11(a1[1],a2[1],ar+2); m[2] = ar[2]; m[3] = ar[3];
! 2193:
! 2194: /* (ar[5],ar[4]) = a1[2]*a2[2] */
! 2195: _mulup2_11(a1[2],a2[2],ar+4);
! 2196:
! 2197: /*
! 2198: * (ar[3],ar[2],ar[1],ar[0]) = (a1[1],a1[0])*(a2[1],a2[0])
! 2199: * a1[1]*a2[1] is already in (ar[3],ar[2])
! 2200: */
! 2201: /* c = (a1[1]+a1[0])*(a2[1]+a2[0]) */
! 2202: t = a1[1]^a1[0]; s = a2[1]^a2[0]; _mulup2_11(t,s,c);
! 2203: /* c += (ar[1],ar[0])+(ar[3],ar[2]) */
! 2204: c[0] ^= ar[2]^ar[0]; c[1] ^= ar[3]^ar[1];
! 2205: /* (ar[2],ar[1]) += c */
! 2206: ar[1] ^= c[0]; ar[2] ^= c[1];
! 2207:
! 2208: /*
! 2209: * m = (a1[1],a1[2]+a1[0])*(a2[1],a2[2]+a2[0])
! 2210: * a1[1]*a2[1] is already in hi(m)
! 2211: *
! 2212: */
! 2213: /* low(m) = (a1[0]+a1[2])*(a2[0]+a2[2]) */
! 2214: t = a1[2]^a1[0]; s = a2[2]^a2[0]; _mulup2_11(t,s,m);
! 2215: /* c = (a1[1]+(a1[0]+a1[2]))*(a2[1]+(a2[0]+a2[2])) */
! 2216: t ^= a1[1]; s ^= a2[1]; _mulup2_11(t,s,c);
! 2217: /* c += low(m)+hi(m) */
! 2218: c[0] ^= m[2]^m[0]; c[1] ^= m[3]^m[1];
! 2219: /* mid(m) += c */
! 2220: m[1] ^= c[0]; m[2] ^= c[1];
! 2221:
! 2222: /* m += (ar[3],ar[2],ar[1],ar[0])+(ar[5],ar[4]) */
! 2223: m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];
! 2224: m[2] ^= ar[2]; m[3] ^= ar[3];
! 2225:
! 2226: /* (ar[5],ar[4],ar[3],ar[2]) += m */
! 2227: ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
! 2228: }
! 2229: #endif
! 2230:
! 2231: void _mulup2_44(a1,a2,ar)
! 2232: unsigned int *a1,*a2,*ar;
! 2233: {
! 2234: unsigned int m[4];
! 2235: unsigned int c1[2],c2[2];
! 2236:
! 2237: c1[0] = a1[2]^a1[0]; c1[1] = a1[3]^a1[1];
! 2238: c2[0] = a2[2]^a2[0]; c2[1] = a2[3]^a2[1];
! 2239:
! 2240: _mulup2_22(a1,a2,ar);
! 2241: _mulup2_22(a1+2,a2+2,ar+4);
! 2242: _mulup2_22(c1,c2,m);
! 2243:
! 2244: m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];
! 2245: m[2] ^= ar[2]^ar[6]; m[3] ^= ar[3]^ar[7];
! 2246:
! 2247: ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
! 2248: }
! 2249:
! 2250: void _mulup2_55(a1,a2,ar)
! 2251: unsigned int *a1,*a2,*ar;
! 2252: {
! 2253: unsigned int m[6];
! 2254: unsigned int c1[3],c2[3];
! 2255:
! 2256: c1[0]=a1[3]^a1[0]; c1[1]=a1[4]^a1[1]; c1[2]=a1[2];
! 2257: c2[0]=a2[3]^a2[0]; c2[1]=a2[4]^a2[1]; c2[2]=a2[2];
! 2258:
! 2259: _mulup2_33(a1,a2,ar);
! 2260: _mulup2_22(a1+3,a2+3,ar+6);
! 2261: _mulup2_33(c1,c2,m);
! 2262:
! 2263: m[0] ^= ar[0]^ar[6]; m[1] ^= ar[1]^ar[7]; m[2] ^= ar[2]^ar[8];
! 2264: m[3] ^= ar[3]^ar[9]; m[4] ^= ar[4]; m[5] ^= ar[5];
! 2265:
! 2266: ar[3] ^= m[0]; ar[4] ^= m[1]; ar[5] ^= m[2];
! 2267: ar[6] ^= m[3]; ar[7] ^= m[4]; ar[8] ^= m[5];
! 2268: }
! 2269:
! 2270: void _mulup2_66(a1,a2,ar)
! 2271: unsigned int *a1,*a2,*ar;
! 2272: {
! 2273: unsigned int m[6];
! 2274: unsigned int c1[3],c2[3];
! 2275:
! 2276: c1[0]=a1[3]^a1[0]; c1[1]=a1[4]^a1[1]; c1[2]=a1[5]^a1[2];
! 2277: c2[0]=a2[3]^a2[0]; c2[1]=a2[4]^a2[1]; c2[2]=a2[5]^a2[2];
! 2278:
! 2279: _mulup2_33(a1,a2,ar);
! 2280: _mulup2_33(a1+3,a2+3,ar+6);
! 2281: _mulup2_33(c1,c2,m);
! 2282:
! 2283: m[0] ^= ar[0]^ar[6]; m[1] ^= ar[1]^ar[7]; m[2] ^= ar[2]^ar[8];
! 2284: m[3] ^= ar[3]^ar[9]; m[4] ^= ar[4]^ar[10]; m[5] ^= ar[5]^ar[11];
! 2285:
! 2286: ar[3] ^= m[0]; ar[4] ^= m[1]; ar[5] ^= m[2];
! 2287: ar[6] ^= m[3]; ar[7] ^= m[4]; ar[8] ^= m[5];
! 2288: }
! 2289:
! 2290: #if 0
! 2291: void printup2_(f,l)
! 2292: unsigned int *f;
! 2293: int l;
! 2294: {
! 2295: int i;
! 2296:
! 2297: for ( i = 32*l-1; i >= 0; i-- ) {
! 2298: if ( f[i>>5]&(1<<(i&31)) )
! 2299: fprintf(stderr,"+x^%d",i);
! 2300: }
! 2301: fprintf(stderr,"\n");
! 2302: }
! 2303: #endif
! 2304:
! 2305: void type1_bin_invup2(a,n,inv)
! 2306: UP2 a;
! 2307: int n;
! 2308: UP2 *inv;
! 2309: {
! 2310: int lf,lg,i,j,k,lg2,df,dg,l,w;
! 2311: unsigned int r;
! 2312: UP2 wf,wg,wb,wc,ret,t;
! 2313: unsigned int *p;
! 2314:
! 2315: lf = a->w;
! 2316: lg = (n>>5)+1;
! 2317:
! 2318: W_NEWUP2(wf,lf+1); _copyup2(a,wf);
! 2319: W_NEWUP2(wg,lg+1); wg->w = lg;
! 2320: for ( i = lg-1, p = wg->b; i>=0; i-- )
! 2321: p[i] = 0xffffffff;
! 2322: if ( r = (n+1)&31 )
! 2323: p[lg-1] &= (1<<r)-1;
! 2324: lg2 = lg*2;
! 2325: W_NEWUP2(wb,lg2+2); wb->w = 1; wb->b[0] = 1;
! 2326: W_NEWUP2(wc,lg2+2); wc->w = 0;
! 2327: k = 0;
! 2328: df = degup2(wf); dg = degup2(wg);
! 2329: while ( 1 ) {
! 2330: lf = wf->w;
! 2331: p = wf->b;
! 2332: for ( j = 0; j < lf && !p[j]; j++ );
! 2333: for ( l = j<<5, w = p[j]; !(w&1); w >>=1, l++ );
! 2334: _bshiftup2_destructive(wf,l);
! 2335: _bshiftup2_destructive(wc,-l);
! 2336: k += l;
! 2337: df -= l;
! 2338: if ( wf->w == 1 && wf->b[0] == 1 ) {
! 2339: k %= (n+1);
! 2340: if ( k != 1 ) {
! 2341: _bshiftup2_destructive(wb,k-(n+1));
! 2342: remup2_type1_destructive(wb,n);
! 2343: }
! 2344: NEWUP2(ret,wb->w);
! 2345: _copyup2(wb,ret);
! 2346: *inv = ret;
! 2347: return;
! 2348: }
! 2349: if ( df < dg ) {
! 2350: i = df; df = dg; dg = i;
! 2351: t = wf; wf = wg; wg = t;
! 2352: t = wb; wb = wc; wc = t;
! 2353: }
! 2354: /* df >= dg */
! 2355: _addup2_destructive(wf,wg);
! 2356: df = degup2(wf);
! 2357: _addup2_destructive(wb,wc);
! 2358: }
! 2359: }
! 2360:
! 2361: UP2 *compute_tab_gf2n(f)
! 2362: UP2 f;
! 2363: {
! 2364: GEN_UP2 mod;
! 2365: int m,n,w,i;
! 2366: UP2 t;
! 2367: UP2 *tab;
! 2368:
! 2369: mod = current_mod_gf2n;
! 2370: m = degup2(mod->dense);
! 2371: n = degup2(f);
! 2372: w = f->w;
! 2373: tab = (UP2 *)CALLOC(m,sizeof(UP2));
! 2374: NEWUP2(tab[0],w); tab[0]->w = 1; tab[0]->b[0] = 2;
! 2375: for ( i = 1; i < m; i++ ) {
! 2376: squareup2(tab[i-1],&t); remup2(t,f,&tab[i]);
! 2377: }
! 2378: return tab;
! 2379: }
! 2380:
! 2381: UP compute_trace_gf2n(tab,c,n)
! 2382: UP2 *tab;
! 2383: GF2N c;
! 2384: int n;
! 2385: {
! 2386: GEN_UP2 mod;
! 2387: int w,m,i,j;
! 2388: UP2 *trace;
! 2389: UP2 c1,c2,s,t;
! 2390: UP r;
! 2391: GF2N a;
! 2392:
! 2393: mod = current_mod_gf2n;
! 2394: w = mod->dense->w;
! 2395: m = degup2(mod->dense);
! 2396: trace = (UP2 *)ALLOCA(n*sizeof(UP2));
! 2397: for ( i = 0; i < n; i++ ) {
! 2398: W_NEWUP2(trace[i],w); trace[i]->w = 0;
! 2399: }
! 2400: W_NEWUP2(c1,2*w); _copyup2(c->body,c1);
! 2401: W_NEWUP2(c2,2*w);
! 2402:
! 2403: for ( i = 0; i < m; i++ ) {
! 2404: t = tab[i];
! 2405: /* trace += c^(2^i)*tab[i] */
! 2406: for ( j = degup2(t); j >= 0; j-- )
! 2407: if ( t->b[j/BSH] & (1<<(j%BSH)) )
! 2408: _addup2_destructive(trace[j],c1);
! 2409: _squareup2(c1,c2); gen_simpup2_destructive(c2,mod);
! 2410: t = c1; c1 = c2; c2 = t;
! 2411: }
! 2412: for ( i = n-1; i >= 0 && !trace[i]->w; i-- );
! 2413: if ( i < 0 )
! 2414: r = 0;
! 2415: else {
! 2416: r = UPALLOC(i); r->d = i;
! 2417: for ( j = 0; j <= i; j++ )
! 2418: if ( t = trace[j] ) {
! 2419: NEWUP2(s,t->w); _copyup2(t,s);
! 2420: MKGF2N(s,a); r->c[j] = (Num)a;
! 2421: }
! 2422: }
! 2423: return r;
! 2424: }
! 2425:
! 2426: void up2toup(f,r)
! 2427: UP2 f;
! 2428: UP *r;
! 2429: {
! 2430: int d,i;
! 2431: UP2 c;
! 2432: UP t;
! 2433: GF2N a;
! 2434:
! 2435: if ( !f || !f->w )
! 2436: *r = 0;
! 2437: else {
! 2438: d = degup2(f);
! 2439: t = UPALLOC(d); *r = t;
! 2440: t->d = d;
! 2441: for ( i = 0; i <= d; i++ ) {
! 2442: if ( f->b[i/BSH] & (1<<(i%BSH)) ) {
! 2443: NEWUP2(c,1);
! 2444: c->w = 1; c->b[0] = 1;
! 2445: MKGF2N(c,a);
! 2446: t->c[i] = (Num)a;
! 2447: }
! 2448: }
! 2449: }
! 2450: }
! 2451:
! 2452: void find_root_up2(f,r)
! 2453: UP2 f;
! 2454: GF2N *r;
! 2455: {
! 2456: int n;
! 2457: UP2 *tab;
! 2458: UP uf,trace,gcd,quo,rem;
! 2459: GF2N c;
! 2460: int i;
! 2461:
! 2462: /* computeation of tab : tab[i] = t^(2^i) mod f (i=0,...,n-1) */
! 2463: n = degup2(f);
! 2464: tab = compute_tab_gf2n(f);
! 2465: up2toup(f,&uf);
! 2466: while ( uf->d > 1 ) {
! 2467: randomgf2n(&c);
! 2468: if ( !c )
! 2469: continue;
! 2470: /* trace = (ct)+(ct)^2+...+(ct)^(2^(n-1)) */
! 2471: trace = compute_trace_gf2n(tab,c,n);
! 2472: gcdup(trace,uf,&gcd);
! 2473: if ( gcd->d > 0 && gcd->d < uf->d ) {
! 2474: /* fprintf(stderr,"deg(GCD)=%d\n",gcd->d); */
! 2475: if ( 2*gcd->d > uf->d ) {
! 2476: qrup(uf,gcd,&quo,&rem);
! 2477: uf = quo;
! 2478: } else {
! 2479: uf = gcd;
! 2480: }
! 2481: }
! 2482: }
! 2483: /* now deg(uf) = 1. The root of uf is ug->c[0]/uf->c[1] */
! 2484: divgf2n((GF2N)uf->c[0],(GF2N)uf->c[1],r);
! 2485: }