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