![[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|
Return to up2.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
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 ! fujimoto 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/up2.c,v 1.5 2015/08/06 10:01:52 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: {\
100: unsigned int _a,_b,_c,_d,_t,_s;\
101: unsigned int _ll,_cc;\
102: _a = (a);\
103: _b = (b);\
104: _c = _a&0xff00ff00; /* c = a4 0 a2 0 */\
105: _a ^= _c; /* a = 0 a3 0 a1 */\
106: _d = _b&0xff00ff00; /* d = 0 0 b2 0 */\
107: _b ^= _d; /* b = 0 0 0 b1 */\
108: _c |= (_d>>8); /* c = a4 00 a2 b2 */\
109: _b |= (_a<<8); /* b = a3 00 a1 b1 */\
110: _a = (_c>>16); /* a = a4 00 */\
111: _c &= 0xffff; /* c = a2 b2 */\
112: _d = (_b>>16); /* d = a3 00 */\
113: _b &= 0xffff; /* b = a1 b1 */\
114: _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\
115: _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\
116: _a ^= _c; _d^= _b;\
117: _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
118: _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
119: _cc ^= _ll;\
120: (low) = (unsigned int)(_ll^(_cc<<16));\
121: (high) = (unsigned int)(_cc>>16);\
122: }
123:
124: /*
125: * 32bit x 32bit -> 64bit
126: * This is an inline expansion of 4byte x 4byte Karatsuba multiplication
127: * Necessary indices for up2_bbtab are efficiently generated.
128: */
129:
130: #define GF2M_MUL_1(high,low,a,b)\
131: {\
132: unsigned int _a,_b,_c,_d,_t,_s;\
133: unsigned int _uu,_ll,_cc;\
134: _a = (a);\
135: _b = (b);\
136: _c = _a&0xff00ff00; /* _c = a4 0 a2 0 */\
137: _a ^= _c; /* a = 0 a3 0 a1 */\
138: _d = _b&0xff00ff00; /* _d = b4 0 b2 0 */\
139: _b ^= _d; /* b = 0 b3 0 b1 */\
140: _c |= (_d>>8); /* _c = a4 b4 a2 b2 */\
141: _b |= (_a<<8); /* b = a3 b3 a1 b1 */\
142: _a = (_c>>16); /* a = a4 b4 */\
143: _c &= 0xffff; /* _c = a2 b2 */\
144: _d = (_b>>16); /* _d = a3 b3 */\
145: _b &= 0xffff; /* b = a1 b1 */\
146: _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
147: _uu = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
148: _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\
149: _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\
150: _a ^= _c; _d^= _b;\
151: _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
152: _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
153: _cc ^= (_ll^_uu);\
154: (low) = (unsigned int)(_ll^(_cc<<16));\
155: (high) = (unsigned int)(_uu^(_cc>>16));\
156: }
157: #endif
158:
159: #define REMUP2_ONESTEP(a,q,s,n,w)\
160: if ( !s ) a[q] ^= w;\
161: else {\
162: a[q] ^= (w<<s);\
163: if ( q+1 <= n )\
164: a[q+1] ^= (w>>(32-s));\
165: }
166:
1.4 noro 167: void ptoup2(P n,UP2 *nr)
1.1 noro 168: {
169: DCP dc;
170: UP2 r,s;
171: int d,w,i;
172:
173: if ( !n )
174: *nr = 0;
175: else if ( OID(n) == O_N )
176: if ( 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: }
194: }
195:
1.4 noro 196: void ptoup2_sparse(P n,UP2 *nr)
1.1 noro 197: {
198: DCP dc;
199: UP2 s;
200: int d,i;
201:
202: if ( !n )
203: *nr = 0;
204: else if ( OID(n) == O_N )
205: if ( 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: }
221: }
222:
1.4 noro 223: void up2top(UP2 n,P *nr)
1.1 noro 224: {
225: int i,d;
226: DCP dc0,dc;
227:
228: if ( !n )
229: *nr = 0;
230: else if ( !(d = degup2(n)) )
231: *nr = (P)ONE;
232: else {
233: for ( i = d, dc0 = 0; i >= 0; i-- )
234: if ( n->b[i/BSH] & (1<<(i%BSH)) ) {
235: NEXTDC(dc0,dc); STOQ(i,DEG(dc)); COEF(dc) = (P)ONE;
236: }
237: if ( !up2_var )
238: up2_var = CO->v;
239: MKP(up2_var,dc0,*nr);
240: }
241: }
242:
1.4 noro 243: void up2tovect(UP2 n,VECT *nr)
1.1 noro 244: {
245: int i,d;
246: VECT v;
247:
248: if ( !n )
249: *nr = 0;
250: else {
251: d = degup2(n);
252: MKVECT(v,d+1); *nr = v;
253: for ( i = d; i >= 0; i-- )
254: if ( n->b[i/BSH] & (1<<(i%BSH)) )
255: v->body[i] = (pointer)ONE;
256: }
257: }
258:
1.4 noro 259: void up2ton(UP2 p,Q *n)
1.1 noro 260: {
261: N nm;
262: int w;
263:
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: }
273: }
274:
1.4 noro 275: void ntoup2(Q n,UP2 *p)
1.1 noro 276: {
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: }
288: }
289:
1.4 noro 290: void gen_simpup2(UP2 p,GEN_UP2 m,UP2 *r)
1.1 noro 291: {
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;
300: }
301:
1.4 noro 302: void gen_simpup2_destructive(UP2 p,GEN_UP2 m)
1.1 noro 303: {
304: UP2 t;
305:
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);
313: }
314:
1.4 noro 315: void gen_invup2(UP2 p,GEN_UP2 m,UP2 *r)
1.1 noro 316: {
317: if ( !m )
318: error("gen_invup2 : invalid modulus");
319: else
320: invup2(p,m->dense,r);
321: }
322:
1.4 noro 323: void gen_pwrmodup2(UP2 a,Q b,GEN_UP2 m,UP2 *c)
1.1 noro 324: {
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);
331: }
332:
1.4 noro 333: void simpup2(UP2 p,UP2 m,UP2 *r)
1.1 noro 334: {
335: if ( !lm_lazy && m )
336: remup2(p,m,r);
337: else
338: *r = p;
339: }
340:
1.4 noro 341: int degup2(UP2 a)
1.1 noro 342: {
343: unsigned int l,i,t;
344:
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: }
352: }
353:
1.4 noro 354: int degup2_sparse(UP2 a)
1.1 noro 355: {
356: if ( !a || !a->w )
357: return -1;
358: else
359: return a->b[0];
360: }
361:
1.4 noro 362: int degup2_1(unsigned int a)
1.1 noro 363: {
364: int i;
365:
366: for ( i = 0; a; a>>=1, i++);
367: return i-1;
368: }
369:
1.4 noro 370: void addup2(UP2 a,UP2 b,UP2 *c)
1.1 noro 371: {
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: }
399: }
400:
1.4 noro 401: void subup2(UP2 a,UP2 b,UP2 *c)
1.1 noro 402: {
403: addup2(a,b,c);
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: {
418: int i;
419: unsigned int _u,_l,t;
420:
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: }
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: {
430: int i;
431: unsigned int u,l;
432:
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: }
437: }
438:
1.4 noro 439: void _mulup2_1(UP2 a,unsigned int b,UP2 c)
1.1 noro 440: {
441: int w;
442:
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: }
453: }
454:
1.4 noro 455: void _mulup2_h(UP2 a,unsigned int b,UP2 c)
1.1 noro 456: {
457: int w;
458:
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: }
466: }
467:
1.4 noro 468: void mulup2(UP2 a,UP2 b,UP2 *c)
1.1 noro 469: {
470: UP2 t;
471: int wa,wb,w;
472:
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: }
497: }
498:
1.4 noro 499: void _kmulup2_(unsigned int *a,unsigned int *b,int w,unsigned int *c)
1.1 noro 500: {
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: }
510: }
511:
1.4 noro 512: void _mulup2_nn(unsigned int *a,unsigned int *b,int w,unsigned int *c)
1.1 noro 513: {
514: int wlow,whigh;
515: struct _oUP2 ablow,abhigh,alow,ahigh,blow,bhigh,aa,bb,mid,cmid;
516:
517: /* wlow >= whigh */
518: wlow = (w+1)/2; whigh = w-wlow;
519:
520: alow.w = wlow; alow.b = a;
521: ahigh.w = whigh; ahigh.b = a+wlow;
522:
523: blow.w = wlow; blow.b = b;
524: bhigh.w = whigh; bhigh.b = b+wlow;
525:
526: ablow.b = c;
527: abhigh.b = c+wlow*2;
528:
529: _kmulup2_(a,b,wlow,ablow.b); ablow.w = 2*wlow;
530: _kmulup2_(a+wlow,b+wlow,whigh,abhigh.b); abhigh.w = 2*whigh;
531:
532: W_NEW_UP2(aa,wlow);
533: _addup2_(&alow,&ahigh,&aa); aa.w = wlow;
534:
535: W_NEW_UP2(bb,wlow);
536: _addup2_(&blow,&bhigh,&bb); bb.w = wlow;
537:
538: W_NEW_UP2(mid,2*wlow);
539: _kmulup2_(aa.b,bb.b,wlow,mid.b); mid.w = 2*wlow;
540:
541: _addtoup2_(&ablow,&mid); _addtoup2_(&abhigh,&mid);
542:
543: cmid.w = 2*wlow; cmid.b = c+wlow;
544: _addtoup2_(&mid,&cmid);
545: }
546:
1.4 noro 547: void _mulup2(UP2 a,UP2 b,UP2 c)
1.1 noro 548: {
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);
565: }
566:
1.4 noro 567: void _mulup2_(_UP2 a,_UP2 b,_UP2 c)
1.1 noro 568: {
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);
585: }
586:
1.4 noro 587: void squareup2(UP2 n,UP2 *nr)
1.1 noro 588: {
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: }
608: }
609:
1.4 noro 610: void _squareup2(UP2 n,UP2 nr)
1.1 noro 611: {
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: }
631: }
632:
1.4 noro 633: void _adjup2(UP2 n)
1.1 noro 634: {
635: int i;
636: unsigned int *nb;
637:
638: nb = n->b;
639: for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );
640: i++;
641: n->w = i;
642: }
643:
1.4 noro 644: void _adjup2_(_UP2 n)
1.1 noro 645: {
646: int i;
647: unsigned int *nb;
648:
649: nb = n->b;
650: for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );
651: i++;
652: n->w = i;
653: }
654:
1.4 noro 655: void _addup2(UP2 a,UP2 b,UP2 c)
1.1 noro 656: {
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);
678: }
679:
680: /* a += b */
681:
1.4 noro 682: void _addup2_destructive(UP2 a,UP2 b)
1.1 noro 683: {
684: int i,wa,wb;
685: unsigned int *ab,*bb;
686:
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: }
706: }
707:
1.4 noro 708: void _addup2_(_UP2 a,_UP2 b,_UP2 c)
1.1 noro 709: {
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;
725: }
726:
1.4 noro 727: void _addtoup2_(_UP2 a,_UP2 b)
1.1 noro 728: {
729: int i,wa;
730: unsigned int *ab,*bb;
731:
732: wa = a->w;
733: ab = a->b; bb = b->b;
734: for ( i = 0; i < wa; i++ )
735: *bb++ ^= *ab++;
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: {
742: unsigned int t;
743:
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: }
756: }
757:
758: void init_up2_tab()
759: {
760: unsigned int i,j;
761:
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);
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: {
778: unsigned int q,i;
779:
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;
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: {
791: int i;
792: unsigned t,q;
793:
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;
802: }
803:
1.4 noro 804: void qrup2(UP2 a,UP2 b,UP2 *q,UP2 *r)
1.1 noro 805: {
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: }
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: {
856: unsigned int msa,msb,t,q0;
857: int s,i,wq,wb;
858:
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;
892: }
893:
1.4 noro 894: void remup2(UP2 a,UP2 b,UP2 *c)
1.1 noro 895: {
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: }
933: }
934:
1.4 noro 935: void _remup2(UP2 a,UP2 b,UP2 c)
1.1 noro 936: {
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: }
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: {
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: }
1047: }
1048:
1.4 noro 1049: void remup2_sparse_destructive(UP2 a,UP2 b)
1.1 noro 1050: {
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: }
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: {
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);
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: {
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);
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: {
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);
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: {
1206: unsigned int p1,p2,p3,q1,q2,q3,g1,g2,q,r;
1207: int d1,d2;
1208:
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;
1222: }
1223:
1.4 noro 1224: void _gcdup2_1(unsigned int f1,unsigned int f2,unsigned int *gcd)
1.1 noro 1225: {
1226: unsigned int g1,g2,q,r;
1227: int d1,d2;
1228:
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: }
1246: }
1247:
1248: static struct oEGT eg_a,eg_b;
1249:
1250: void up2_init_eg() {
1251: init_eg(&eg_a); init_eg(&eg_b);
1252: }
1253:
1254: void up2_show_eg() {
1255: print_eg("A",&eg_a);
1256: print_eg("B",&eg_b);
1257: printf("\n");
1258: }
1259:
1.4 noro 1260: void invup2(UP2 a,UP2 m,UP2 *inv)
1.1 noro 1261: {
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: }
1332: }
1333:
1.4 noro 1334: void gcdup2(UP2 a,UP2 m,UP2 *gcd)
1.1 noro 1335: {
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: }
1402: }
1403:
1.4 noro 1404: void chsgnup2(UP2 a,UP2 *c)
1.1 noro 1405: {
1406: *c = a;
1407: }
1408:
1.4 noro 1409: void pwrmodup2(UP2 a,Q b,UP2 m,UP2 *c)
1.1 noro 1410: {
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: }
1443: }
1444:
1.4 noro 1445: void pwrmodup2_sparse(UP2 a,Q b,UP2 m,UP2 *c)
1.1 noro 1446: {
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: }
1479: }
1480:
1.4 noro 1481: int compup2(UP2 n1,UP2 n2)
1.1 noro 1482: {
1483: int i;
1484: unsigned int *m1,*m2;
1485:
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: }
1506: }
1507:
1.4 noro 1508: void _copyup2(UP2 n,UP2 r)
1.1 noro 1509: {
1510: r->w = n->w;
1511: bcopy(n->b,r->b,n->w*sizeof(unsigned int));
1512: }
1513:
1.4 noro 1514: void _bshiftup2(UP2 n,int b,UP2 r)
1.1 noro 1515: {
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: }
1578: }
1579:
1.4 noro 1580: void _bshiftup2_destructive(UP2 n,int b)
1.1 noro 1581: {
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: }
1643: }
1644:
1.4 noro 1645: void diffup2(UP2 f,UP2 *r)
1.1 noro 1646: {
1647: int d,i,w;
1648: UP2 t;
1649:
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: }
1664: }
1665:
1.4 noro 1666: int sqfrcheckup2(UP2 f)
1.1 noro 1667: {
1668: UP2 df,g;
1669:
1670: diffup2(f,&df);
1671: gcdup2(f,df,&g);
1672: if ( degup2(g) >= 1 )
1673: return 0;
1674: else
1675: return 1;
1676: }
1677:
1.4 noro 1678: int irredcheckup2(UP2 f)
1.1 noro 1679: {
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;
1716: }
1717:
1.4 noro 1718: int irredcheck_dddup2(UP2 f)
1.1 noro 1719: {
1720: UP2 x,u,t,s,gcd;
1721: int n,i;
1722:
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;
1736: }
1737:
1.4 noro 1738: void _copy_up2bits(UP2 p,unsigned int **mat,int pos)
1.1 noro 1739: {
1740: int d,col,j,jcol,jsh;
1741: unsigned int bit;
1742:
1743: d = degup2(p);
1744: col = pos/BSH; bit = 1<<(pos%BSH);
1745:
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;
1752: }
1753:
1.4 noro 1754: int compute_multiplication_matrix(P p0,GF2MAT *mp)
1.1 noro 1755: {
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;
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: {
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);
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: {
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);
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: {
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;
1890: }
1891:
1.4 noro 1892: void mul_nb(GF2MAT mat,unsigned int *a,unsigned int *b,unsigned int *c)
1.1 noro 1893: {
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: }
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: {
1923: int r,w,i;
1924: unsigned int msb;
1925:
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;
1934: }
1935:
1.4 noro 1936: void mat_to_gf2mat(MAT a,unsigned int ***b)
1.1 noro 1937: {
1938: int n,w,i,j;
1939: unsigned int **m;
1940:
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);
1948: }
1949:
1.4 noro 1950: void gf2mat_to_mat(unsigned int **a,int n,MAT *b)
1.1 noro 1951: {
1952: int w,i,j;
1953: MAT m;
1954:
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;
1961: }
1962:
1.4 noro 1963: void mulgf2mat(int n,unsigned int **a,unsigned int **b,unsigned int **c)
1.1 noro 1964: {
1965: int i,j,k,w;
1966:
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: }
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: {
1981: int j,k,w;
1982:
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];
1989: }
1990:
1.4 noro 1991: int mulgf2vectvect(int n,unsigned int *a,unsigned int *b)
1.1 noro 1992: {
1993: unsigned int t,r;
1994: int i,w;
1995:
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;
2001: }
2002:
1.4 noro 2003: int invgf2mat(int n,unsigned int **a,unsigned int **b)
1.1 noro 2004: {
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;
2040: }
2041:
1.4 noro 2042: INLINE void _mulup2_11(unsigned int a1,unsigned int a2,unsigned int *ar)
1.1 noro 2043: {
2044: GF2M_MUL_1(ar[1],ar[0],a1,a2);
2045: }
2046:
1.4 noro 2047: void _mulup2_22(unsigned int *a1,unsigned int *a2,unsigned int *ar)
1.1 noro 2048: {
2049: unsigned int m[2];
2050:
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];
2055:
2056: ar[1] ^= m[0]; ar[2] ^= m[1];
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: {
2062: unsigned int m[4];
2063: unsigned int c1[2],c2[2];
2064:
2065: c1[0] = a1[2]^a1[0]; c1[1] = a1[1];
2066: c2[0] = a2[2]^a2[0]; c2[1] = a2[1];
2067:
2068: _mulup2_22(a1,a2,ar);
2069: _mulup2_11(*(a1+2),*(a2+2),ar+4);
2070: _mulup2_22(c1,c2,m);
2071:
2072: m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];
2073: m[2] ^= ar[2]; m[3] ^= ar[3];
2074:
2075: ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
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: {
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];
2123:
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];
2126: }
2127: #endif
2128:
1.4 noro 2129: void _mulup2_44(unsigned int *a1,unsigned int *a2,unsigned int *ar)
1.1 noro 2130: {
2131: unsigned int m[4];
2132: unsigned int c1[2],c2[2];
2133:
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];
2136:
2137: _mulup2_22(a1,a2,ar);
2138: _mulup2_22(a1+2,a2+2,ar+4);
2139: _mulup2_22(c1,c2,m);
2140:
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];
2143:
2144: ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
2145: }
2146:
1.4 noro 2147: void _mulup2_55(unsigned int *a1,unsigned int *a2,unsigned int *ar)
1.1 noro 2148: {
2149: unsigned int m[6];
2150: unsigned int c1[3],c2[3];
2151:
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];
2154:
2155: _mulup2_33(a1,a2,ar);
2156: _mulup2_22(a1+3,a2+3,ar+6);
2157: _mulup2_33(c1,c2,m);
2158:
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];
2161:
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];
2164: }
2165:
1.4 noro 2166: void _mulup2_66(unsigned int *a1,unsigned int *a2,unsigned int *ar)
1.1 noro 2167: {
2168: unsigned int m[6];
2169: unsigned int c1[3],c2[3];
2170:
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];
2173:
2174: _mulup2_33(a1,a2,ar);
2175: _mulup2_33(a1+3,a2+3,ar+6);
2176: _mulup2_33(c1,c2,m);
2177:
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];
2180:
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];
2183: }
2184:
2185: #if 0
1.4 noro 2186: void printup2_(unsigned int *f,int l)
1.1 noro 2187: {
2188: int i;
2189:
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");
2195: }
2196: #endif
2197:
1.4 noro 2198: void type1_bin_invup2(UP2 a,int n,UP2 *inv)
1.1 noro 2199: {
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: }
2249: }
2250:
1.4 noro 2251: UP2 *compute_tab_gf2n(UP2 f)
1.1 noro 2252: {
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;
2268: }
2269:
1.4 noro 2270: UP compute_trace_gf2n(UP2 *tab,GF2N c,int n)
1.1 noro 2271: {
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;
2310: }
2311:
1.4 noro 2312: void up2toup(UP2 f,UP *r)
1.1 noro 2313: {
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: }
2334: }
2335:
1.4 noro 2336: void find_root_up2(UP2 f,GF2N *r)
1.1 noro 2337: {
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);
2366: }