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