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