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