[BACK]Return to up2.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/up2.c, Revision 1.3

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3     ! noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.3     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/up2.c,v 1.2 2000/08/21 08:31:28 noro Exp $
1.2       noro       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: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>