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