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