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

Annotation of OpenXM_contrib2/asir2000/engine/H.c, Revision 1.6

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.6     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/H.c,v 1.5 2001/07/04 07:19:20 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "inline.h"
                     52: #include "base.h"
                     53: #include <math.h>
                     54:
                     55: #if 1
                     56: #define Mulum mulum
                     57: #define Divum divum
                     58: #define Mulsum mulsum
                     59: #define Gcdum gcdum
                     60: #endif
                     61:
                     62: #define FCTR 0
                     63: #define SQFR 1
                     64: #define DDD 2
                     65: #define NEWDDD 3
                     66:
                     67: LUM LUMALLOC();
                     68:
1.6     ! noro       69: void berle(int index,int count,P f,ML *listp)
1.1       noro       70: {
                     71:        UM wf,wf1,wf2,wfs,gcd;
                     72:        ML flist;
                     73:        int fn,fn1,fm,m,n,fhd;
                     74:        register int i,j,inv,hd,*ptr,*ptr1;
                     75:
                     76:        n = UDEG(f);
                     77:        wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n);
                     78:        wfs = W_UMALLOC(n); gcd = W_UMALLOC(n);
                     79:        for ( j = 0, fn = n + 1; (j < count) && (fn > 1); ) {
                     80:                m = sprime[index++];
                     81:                if ( !rem(NM((Q)UCOEF(f)),m) )
                     82:                        continue;
                     83:                ptoum(m,f,wf); cpyum(wf,wf1);
                     84:                diffum(m,wf1,wf2); gcdum(m,wf1,wf2,gcd);
                     85:                if ( DEG(gcd) > 0 )
                     86:                        continue;
                     87:                hd = COEF(wf)[n]; inv = invm(hd,m);
                     88:                for ( i = n, ptr = COEF(wf); i >= 0; i-- )
                     89:                        ptr[i] = ( ptr[i] * inv ) % m;
                     90:                fn1 = berlecnt(m,wf);
                     91:                if ( fn1 < fn ) {
                     92:                        fn = fn1; fm = m; fhd = hd;
                     93:                        for ( i = n, ptr = COEF(wf), ptr1 = COEF(wfs); i >= 0; i-- )
                     94:                                ptr1[i] = ptr[i];
                     95:                }
                     96:                j++;
                     97:        }
                     98:        DEG(wfs) = n;
                     99:        *listp = flist = MLALLOC(fn); flist->n = fn; flist->mod = fm;
                    100: /*     berlemain(fm,wfs,(UM *)flist->c); */
                    101:        if ( fm == 2 )
                    102:                berlemain(fm,wfs,(UM *)flist->c);
                    103:        else
                    104:                newddd(fm,wfs,(UM *)flist->c);
                    105:        for ( i = DEG((UM)(flist->c[0])),
                    106:                ptr = COEF((UM)(flist->c[0])),
                    107:                hd = fhd, m = fm; i >= 0; i-- )
                    108:                ptr[i] = ( ptr[i] * hd ) % m;
                    109: }
                    110:
1.6     ! noro      111: int berlecnt(int mod,UM f)
1.1       noro      112: {
                    113:        register int i,j,**c;
                    114:        int d,dr,n;
                    115:        UM w,q;
                    116:        int **almat();
                    117:
                    118:        n = DEG(f); c = almat(n,n);
                    119:        w = W_UMALLOC(mod + n); q = W_UMALLOC(mod + n);
                    120:        for ( i = 1; ( d = ( mod * i ) ) < n; i++ )
                    121:                c[d][i - 1] = 1;
                    122:        DEG(w) = d; COEF(w)[d] = 1;
                    123:        for ( j = d - 1; j >= 0; j-- )
                    124:                COEF(w)[j] = 0;
                    125:        for ( ; ( i < n ) && ( ( dr = divum(mod,w,f,q) ) != -1 ); i++ ) {
                    126:                for ( j = dr; j >= 0; j-- )
                    127:                        COEF(w)[j + mod] = c[j][i - 1] = COEF(w)[j];
                    128:                for ( j = mod - 1; j >= 0; j-- )
                    129:                        COEF(w)[j] = 0;
                    130:                DEG(w) = dr + mod;
                    131:        }
                    132:        for ( i = 1; i < n; i++ )
                    133:                c[i][i - 1] = ( c[i][i - 1] + mod - 1 ) % mod;
                    134:        return berlecntmain(mod,n,n-1,c);
                    135: }
                    136:
                    137: /* XXX berlecntmain should not be used for large mod */
                    138:
1.6     ! noro      139: int berlecntmain(int mod,int n,int m,int **c)
1.1       noro      140: {
                    141:        register int *p1,*p2,i,j,k,l,a;
                    142:        int *tmp,inv;
                    143:        int cfs;
                    144:
                    145:        for ( cfs = 1, j = k = 0; j < m; j++ ) {
                    146:                for ( i = k; ( n > i ) && ( c[i][j] == 0 ); i++ );
                    147:                if ( i == n ) {
                    148:                        cfs++; continue;
                    149:                }
                    150:                if ( i != k ) {
                    151:                        tmp = c[i]; c[i] = c[k]; c[k] = tmp;
                    152:                }
                    153:                p1 = c[k]; inv = invm((p1[j] + mod) % mod,mod);
                    154:                for ( l = j; l < m; l++ )
                    155:                        p1[l] = ( p1[l] * inv ) % mod;
                    156:                for ( i = k + 1; i < n; c[i][j] = 0, i++ )
                    157:                        if ( i != k && ( a = -c[i][j] ) )
                    158:                                for ( l = j + 1, p2 = c[i]; l < m; l++ )
                    159:                                        p2[l] = (a*p1[l] + p2[l]) % mod;
                    160:                k++;
                    161:        }
                    162:        return ( cfs );
                    163: }
                    164:
1.6     ! noro      165: UM *berlemain(int mod,UM f,UM *fp)
1.1       noro      166: {
                    167:        UM wg,ws,wf,f0,gcd,q;
                    168:        int n;
                    169:        register int i;
                    170:
                    171:        n = DEG(f); wg = W_UMALLOC(n); mini(mod,f,wg);
                    172:        if ( DEG(wg) <= 0 ) {
                    173:                f0 = UMALLOC(n); cpyum(f,f0); *fp++ = f0;
                    174:                return ( fp );
                    175:        }
                    176:        f0 = W_UMALLOC(n); cpyum(f,f0);
                    177:        ws = W_UMALLOC(n); wf = W_UMALLOC(n);
                    178:        q = W_UMALLOC(n); gcd = W_UMALLOC(n);
                    179:        for ( i = 0; i < mod; i++ ) {
                    180:                cpyum(f0,wf); cpyum(wg,ws);
                    181:                COEF(ws)[0] = ( COEF(ws)[0] + mod - i ) % mod;
                    182:                gcdum(mod,wf,ws,gcd);
                    183:                if ( DEG(gcd) > 0 ) {
                    184:                        if ( DEG(gcd) < n ) {
                    185:                                divum(mod,f0,gcd,q); f0 = q; fp = berlemain(mod,gcd,fp);
                    186:                        }
                    187:                        break;
                    188:                }
                    189:        }
                    190:        fp = berlemain(mod,f0,fp);
                    191:        return ( fp );
                    192: }
                    193:
1.6     ! noro      194: void hensel(int index,int count,P f,ML *listp)
1.1       noro      195: {
                    196:        register int i,j;
                    197:        int q,n,bound;
                    198:        int *p;
                    199:        int **pp;
                    200:        ML blist,clist,bqlist,cqlist,rlist;
                    201:        UM *b;
                    202:        LUM fl,tl;
                    203:        LUM *l;
                    204:
                    205:        if ( UDEG(f) == 1 ) {
                    206:                *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
                    207:                return;
                    208:        }
                    209:        berle(index,count,f,&blist);
                    210:        if ( blist->n == 1 ) {
                    211:                *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
                    212:                return;
                    213:        }
                    214:        gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
                    215:        n = bqlist->n; q = bqlist->mod;
                    216:        bqlist->bound = cqlist->bound = bound = mignotte(q,f);
                    217:        if ( bound == 1 ) {
                    218:                *listp = rlist = MLALLOC(n);
                    219:                rlist->n = n; rlist->mod = q; rlist->bound = bound;
                    220:                for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
                    221:                        tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
                    222:                        for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
                    223:                                        pp[j][0] = p[j];
                    224:                }
                    225:        } else {
                    226:                W_LUMALLOC((int)UDEG(f),bound,fl);
                    227:                ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
                    228:        }
                    229: }
                    230:
1.6     ! noro      231: void hensel2(int index,int count,P f,ML *listp)
1.4       noro      232: {
                    233:        register int i,j;
                    234:        int mod,q,n,bound,dx;
                    235:        ML blist,clist,bqlist,cqlist,rlist;
                    236:        UM fm,qfm,gm,qgm,hm,qhm,qam,qbm,w;
                    237:        UM *b;
                    238:        LUM fl,tl;
1.6     ! noro      239:        int k;
1.4       noro      240:
                    241:        dx = UDEG(f);
                    242:        if ( dx == 1 ) {
                    243:                *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
                    244:                return;
                    245:        }
                    246:        berle(index,count,f,&blist);
                    247:        n = blist->n;
                    248:        mod = blist->mod;
                    249:
                    250:        if ( n == 1 ) {
                    251:                *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
                    252:                return;
                    253:        }
                    254:
                    255:        /* find k s.t. mod^k <= 2^27 < mod^(k+1); set q = mod^k */
                    256:        for ( k = 1, q = mod; q <= ((1<<27)/mod); q *= mod, k++ );
                    257:
                    258:        /* mignotte bound */
                    259:        bound = mignotte(q,f);
                    260:
                    261:        *listp = rlist = MLALLOC(n);
                    262:        rlist->n = n;
                    263:        rlist->mod = q;
                    264:        rlist->bound = bound;
                    265:
                    266:        if ( bound == 1 ) {
                    267:                gcdgen(f,blist,&clist);
                    268:                henprep(f,blist,clist,&bqlist,&cqlist);
                    269:
                    270:                for ( i = 0, b = (UM *)bqlist->c; i < n; i++ ) {
                    271:                        COEF(rlist)[i] = tl = LUMALLOC(DEG(b[i]),1);
                    272:                        for ( j = 0; j <= DEG(tl); j++ )
                    273:                                        COEF(tl)[j][0] = COEF(b[i])[j];
                    274:                        COEF(rlist)[i] = tl;
                    275:                }
                    276:        } else {
                    277:                /* fl = f mod q */
                    278:                fl = LUMALLOC(dx,bound);
                    279:                ptolum(q,bound,f,fl);
                    280:                /* fm = f mod mod */
                    281:                fm = W_UMALLOC(dx);
                    282:                ptoum(mod,f,fm);
                    283:                /* fm = f mod q */
                    284:                qfm = W_UMALLOC(dx);
                    285:                ptoum(q,f,qfm);
                    286:
                    287:                gm = W_UMALLOC(dx); qgm = W_UMALLOC(dx);
                    288:                hm = W_UMALLOC(dx); qhm = W_UMALLOC(dx);
                    289:                qam = W_UMALLOC(dx); qbm = W_UMALLOC(dx);
                    290:                w = W_UMALLOC(dx);
                    291:                for ( i = 0; i < n-1; i++ ) {
                    292:                        cpyum(COEF(blist)[i],gm);
                    293:                        cpyum(fm,w);
                    294:                        divum(mod,w,gm,hm);
                    295:
                    296:                        /* find am,bm s.t. qam*qgm+qbm*qhm=1 mod q, qgm=gm mod mod, qhm=hm mod mod */
                    297:                        henprep2(mod,q,k,qfm,gm,hm,qgm,qhm,qam,qbm);
                    298:
                    299:                        henmain2(fl,qgm,qhm,qam,qbm,q,bound,&tl);
                    300:                        rlist->c[i] = (pointer)tl;
                    301:                        cpyum(hm,fm);
                    302:                        cpyum(qhm,qfm);
                    303:                }
                    304:                rlist->c[i] = fl;
                    305:        }
                    306: }
                    307:
                    308: /*
                    309:        f = g0*h0 mod m -> f = gk*hk mod m^(bound), f is replaced by hk
                    310: */
                    311:
1.6     ! noro      312: void henmain2(LUM f,UM g0,UM h0,UM a0,UM b0,int m,int bound,LUM *gp)
1.4       noro      313: {
                    314:        int n,dg,dh,i,k,j,dg1,dh1;
1.6     ! noro      315:        UM wu,wr,ws,wt,q,wh1,wg1,wc,wd,we,wz;
1.4       noro      316:        LUM wb0,wb1,wb2,fk,gk,hk;
                    317:
                    318:        n = DEG(f); dg = DEG(g0); dh = DEG(h0);
                    319:
                    320:        W_LUMALLOC(n,bound,wb0);
                    321:        W_LUMALLOC(n,bound,wb1);
                    322:        W_LUMALLOC(n,bound,wb2);
                    323:
                    324:        wt = W_UMALLOC(2*n); ws = W_UMALLOC(2*n);
                    325:        wr = W_UMALLOC(2*n); wu = W_UMALLOC(2*n);
                    326:        q = W_UMALLOC(2*n);
                    327:        wh1 = W_UMALLOC(2*n); wg1 = W_UMALLOC(2*n);
                    328:
                    329:        /* gk = g0 */
                    330:        gk = LUMALLOC(n,bound);
                    331:        DEG(gk) = dg;
                    332:        for ( i = 0; i <= dg; i++ )
                    333:                COEF(gk)[i][0] = COEF(g0)[i];
                    334:
                    335:        /* hk = h0 */
                    336:        W_LUMALLOC(n,bound,hk);
                    337:        DEG(hk) = dh;
                    338:        for ( i = 0; i <= dh; i++ )
                    339:                COEF(hk)[i][0] = COEF(h0)[i];
                    340:
                    341:        /* fk = gk*hk */
                    342:        W_LUMALLOC(n,bound,fk);
                    343:        mullum(m,bound,gk,hk,fk);
                    344:
                    345:        wc = W_UMALLOC(2*n); wd = W_UMALLOC(2*n);
                    346:        we = W_UMALLOC(2*n); wz = W_UMALLOC(2*n);
                    347:
                    348: #if 0
                    349:        mulum(m,a0,g0,wc);
                    350:        mulum(m,b0,h0,wd);
                    351:        addum(m,wc,wd,wz);
                    352:        if ( DEG(wz) != 0 || COEF(wz)[0] != 1 )
                    353:                error("henmain2 : cannot happen(extgcd)");
                    354: #endif
                    355:
1.5       noro      356: #if 1
1.4       noro      357:        fprintf(stderr,"bound=%d\n",bound);
1.5       noro      358: #endif
1.4       noro      359:        for ( k = 1; k < bound; k++ ) {
1.5       noro      360: #if 1
1.4       noro      361:                fprintf(stderr,".");
1.5       noro      362: #endif
1.4       noro      363:                /* at this point, f = gk*hk mod y^k */
                    364:
                    365: #if 0
                    366:                for ( j = 0; j < k; j++ )
                    367:                        for ( i = 0; i <= n; i++ )
                    368:                                if ( COEF(f)[i][j] != COEF(f)[i][j] )
                    369:                                        error("henmain2 : cannot happen(f=fk)");
                    370: #endif
                    371:
                    372:                /* wt = (f-gk*hk)/y^k */
                    373:                for ( i = 0; i <= n; i++ )
                    374:                        COEF(ws)[i] = COEF(f)[i][k];
                    375:                degum(ws,n);
                    376:                for ( i = 0; i <= n; i++ )
                    377:                        COEF(wu)[i] = COEF(fk)[i][k];
                    378:                degum(wu,n);
                    379:                subum(m,ws,wu,wt);
                    380:
                    381:                /* compute wf1,wg1 s.t. wh1*g0+wg1*h0 = wt */
                    382:                mulum(m,a0,wt,wh1); DEG(wh1) = divum(m,wh1,h0,q);
                    383:                mulum(m,wh1,g0,wc); subum(m,wt,wc,wd); DEG(wd) = divum(m,wd,h0,wg1);
                    384:
                    385:                /* check */
                    386: #if 0
                    387:                if ( DEG(wd) >= 0 || DEG(wg1) > dg )
                    388:                        error("henmain2 : cannot happen(adj)");
                    389:
                    390:                mulum(m,wg1,h0,wc); mulum(m,wh1,g0,wd); addum(m,wc,wd,we);
                    391:                subum(m,we,wt,wz);
                    392:                if ( DEG(wz) >= 0 )
                    393:                        error("henmain2 : cannot happen(coef)");
                    394: #endif
                    395:
                    396:                /* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod m^bound */
                    397:
                    398:                /* wb0 = wh1*y^k */
                    399:                clearlum(n,bound,wb0);
                    400:                DEG(wb0) = dh1 = DEG(wh1);
                    401:                for ( i = 0; i <= dh1; i++ )
                    402:                        COEF(wb0)[i][k] = COEF(wh1)[i];
                    403:
                    404:                /* wb2 = gk*wb0 mod y^bound */
                    405:                clearlum(n,bound,wb2);
                    406:                mullum(m,bound,gk,wb0,wb2);
                    407:
                    408:                /* fk += wb2 */
                    409:                addtolum(m,bound,wb2,fk);
                    410:
                    411:                /* wb1 = wg1*y^k */
                    412:                clearlum(n,bound,wb1);
                    413:                DEG(wb1) = dg1 = DEG(wg1);
                    414:                for ( i = 0; i <= n; i++ )
                    415:                        COEF(wb1)[i][k] = COEF(wg1)[i];
                    416:
                    417:                /* wb2 = hk*wb1 mod y^bound */
                    418:                clearlum(n,bound,wb2);
                    419:                mullum(m,bound,hk,wb1,wb2);
                    420:
                    421:                /* fk += wb2 */
                    422:                addtolum(m,bound,wb2,fk);
                    423:
                    424:                /* fk += wg1*wh1*y^(2*k) mod y^bound) */
                    425:                if ( 2*k < bound ) {
                    426:                        clearlum(n,bound,wb2);
                    427:                        mullum(m,bound,wb0,wb1,wb2);
                    428:                        addtolum(m,bound,wb2,fk);
                    429:                }
                    430:
                    431:                /* gk += wg1*y^k, hk += wh1*y^k */
                    432:                for ( i = 0; i <= DEG(wg1); i++ )
                    433:                        COEF(gk)[i][k] = COEF(wg1)[i];
                    434:                for ( i = 0; i <= DEG(wh1); i++ )
                    435:                        COEF(hk)[i][k] = COEF(wh1)[i];
                    436:        }
1.5       noro      437: #if 1
1.4       noro      438:        fprintf(stderr,"\n");
1.5       noro      439: #endif
1.4       noro      440:        *gp = gk;
                    441:        clearlum(n,bound,f);
                    442:        DEG(f) = dh;
                    443:        for ( i = 0; i <= dh; i++ )
                    444:                for ( j = 0; j < bound; j++ )
                    445:                        COEF(f)[i][j] = COEF(hk)[i][j];
                    446: }
                    447:
1.6     ! noro      448: void clearlum(int n,int bound,LUM f)
1.4       noro      449: {
                    450:        int i;
                    451:
                    452:        for ( i = 0; i <= n; i++ )
                    453:                bzero(COEF(f)[i],bound*sizeof(int));
                    454: }
                    455:
                    456: /* g += f */
                    457:
1.6     ! noro      458: void addtolum(int m,int bound,LUM f,LUM g)
1.4       noro      459: {
                    460:        int n,i;
                    461:
                    462:        n = DEG(f);
                    463:        for ( i = 0; i <= n; i++ )
                    464:                addpadic(m,bound,COEF(f)[i],COEF(g)[i]);
                    465: }
                    466:
1.6     ! noro      467: void hsq(int index,int count,P f,int *nindex,DCP *dcp)
1.1       noro      468: {
                    469:        register int i,j,k;
                    470:        register int **pp,**fpp;
                    471:        register int *px,*py;
                    472:        int **wpp;
                    473:        int n,dr,tmp,m,b,e,np,dt;
                    474:        LUM fpa,wb0,wb1,lcpa,tpa,tlum;
                    475:        struct oDUM *dct;
                    476:        UM wt,wq0,wq,wr,wm,wm0,wa,ws,wb;
                    477:        LUM *llist,*ll;
                    478:        UM *dlist,*l,*c;
                    479:        ML list,fp,cfp;
                    480:        DCP dc;
                    481:
                    482:        sqfrum(index,count,f,nindex,&dct,&fp);
                    483:        np = fp->n; m = fp->mod;
                    484:        if ( ( np == 1 ) && ( dct[0].n == 1 ) ) {
                    485:                NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; *dcp = dc;
                    486:                return;
                    487:        }
                    488:        for ( i = 0, dt = 0; i < np; i++ )
                    489:                dt = MAX(DEG(dct[i].f),dt);
                    490:        b = mig(m,dt,f); fp->bound = b;
                    491:        if ( np == 1 ) {
                    492:                nthrootchk(f,dct,fp,dcp);
                    493:                return;
                    494:        }
                    495:        list = W_MLALLOC(np); list->n = np; list->mod = m; list->bound = 1;
                    496:        for ( i = 0, ll = (LUM *)list->c; i < np; i++ ) {
                    497:                W_LUMALLOC(DEG(dct[i].f),b,ll[i]);
                    498:                for ( j = 0, px = COEF(dct[i].f), pp = COEF(ll[i]);
                    499:                        j <= DEG(ll[i]); j++ )
                    500:                                pp[j][0] = px[j];
                    501:        }
                    502:        dtestsql(f,list,dct,&dc);
                    503:        if ( dc ) {
                    504:                *dcp = dc;
                    505:                return;
                    506:        }
                    507:        n = UDEG(f);
                    508:        W_LUMALLOC(n,b,fpa); W_LUMALLOC(0,b,lcpa);
                    509:        W_LUMALLOC(n,b,wb0); W_LUMALLOC(n,b,wb1);
                    510:        W_LUMALLOC(n,b,tpa);
                    511:        wt = W_UMALLOC(n); ws = W_UMALLOC(n);
                    512:        wr = W_UMALLOC(n);
                    513:        wq = W_UMALLOC(2*n); wq0 = W_UMALLOC(n);
                    514:        wm = W_UMALLOC(2*n); wm0 = W_UMALLOC(2*n);
                    515:        wa = W_UMALLOC(2*n);
                    516:        ptolum(m,b,f,fpa); DEG(lcpa) = 0;
                    517:        for ( i = 0, pp = COEF(lcpa), fpp = COEF(fpa); i < b; i++ )
                    518:                pp[0][i] = fpp[n][i];
                    519:        gcdgen(f,fp,&cfp);
                    520:        llist = (LUM *) ALLOCA(np*sizeof(LUM));
                    521:        dlist = (UM *) ALLOCA(np*sizeof(UM));
                    522:        l = (UM *)fp->c; c = (UM *)cfp->c;
                    523:        for ( i = 0; i < np; i++ ) {
                    524:                W_LUMALLOC(DEG(l[i]),b,llist[i]);
                    525:                for ( j = DEG(l[i]), pp = COEF(llist[i]), px = COEF(l[i]); j >= 0; j-- )
                    526:                        pp[j][0] = px[j];
                    527:                if ( ( e = dct[i].n ) != 1 ) {
                    528:                        wb = dct[i].f;
                    529:                        dlist[i] = W_UMALLOC(DEG(wb)*e); cpyum(l[i],dlist[i]);
                    530:                        divum(m,dlist[i],wb,wq); DEG(dlist[i])= DEG(wq);
                    531:                        for ( k = 0; k <= DEG(wq); k++ )
                    532:                                COEF(dlist[i])[k] = dmb(m,COEF(wq)[k],e,&tmp);
                    533:                }
                    534:        }
                    535:        for ( i = 1; i < b; i++ ) {
                    536:                mullum(m,i+1,lcpa,llist[0],wb0);
                    537:                for ( j = 1; j < np; j++ ) {
                    538:                        mullum(m,i+1,llist[j],wb0,wb1);
                    539:                        tlum = wb0; wb0 = wb1; wb1 = tlum;
                    540:                }
                    541:                for ( j = n, px = COEF(wt), pp = COEF(fpa), wpp = COEF(wb0);
                    542:                        j >= 0; j-- )
                    543:                        px[j] = ( pp[j][i] - wpp[j][i] + m ) % m;
                    544:                degum(wt,n);
                    545:                for ( j = n, px = COEF(wq0); j >= 0; j-- )
                    546:                        px[j] = 0;
                    547:                for ( j = 1; j < np; j++ ) {
                    548:                        mulum(m,wt,c[j],wm); dr = divum(m,wm,l[j],wq);
                    549:                        for ( k = DEG(wq), px = COEF(wq0), py = COEF(wq); k >= 0; k-- )
                    550:                                px[k] = ( px[k] + py[k] ) % m;
                    551:                        for ( k = dr, pp = COEF(llist[j]), px = COEF(wm); k >= 0; k-- )
                    552:                                pp[k][i] = px[k];
                    553:                }
                    554:                degum(wq0,n); mulum(m,wq0,l[0],wm);
                    555:                mulum(m,wt,c[0],wm0); addum(m,wm,wm0,wa);
                    556:                for ( j = DEG(wa), pp = COEF(llist[0]), px = COEF(wa); j >= 0; j-- )
                    557:                        pp[j][i] = px[j];
                    558:                for ( j = n, px = COEF(wq0); j >= 0; j-- )
                    559:                        px[j] = 0;
                    560:                for ( j = 0; j < np; j++ )
                    561:                        if ( dct[j].n == 1 )
                    562:                                for ( k = 0,
                    563:                                        pp = COEF(llist[j]),
                    564:                                        wpp = COEF(((LUM *)list->c)[j]);
                    565:                                        k <= DEG(llist[j]); k++ )
                    566:                                        wpp[k][i] = pp[k][i];
                    567:                        else {
                    568:                                pwrlum(m,i+1,((LUM *)list->c)[j],dct[j].n,tpa);
                    569:                                for ( k = 0,
                    570:                                        pp = COEF(llist[j]),
                    571:                                        wpp = COEF(tpa);
                    572:                                        k <= DEG(l[j]); k++ )
                    573:                                        COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
                    574:                                degum(wt,DEG(l[j])); dr = divum(m,wt,dlist[j],ws);
                    575:                                if ( dr >= 0 ) {
                    576:                                        *dcp = 0;
                    577:                                        return;
                    578:                                } else
                    579:                                        for ( k = 0,
                    580:                                                pp = COEF(((LUM *)list->c)[j]);
                    581:                                                k <= DEG(ws); k++ )
                    582:                                                pp[k][i] = COEF(ws)[k];
                    583:                        }
                    584:                list->bound = i+1; dtestsql(f,list,dct,&dc);
                    585:                if ( dc ) {
                    586:                        *dcp = dc;
                    587:                        return;
                    588:                }
                    589:        }
                    590:        *dcp = 0;
                    591: }
                    592:
1.6     ! noro      593: void gcdgen(P f,ML blist,ML *clistp)
1.1       noro      594: {
                    595:        register int i;
                    596:        int n,d,mod,np;
                    597:        UM wf,wm,wx,wy,wu,wv,wa,wb,wg,q,tum;
                    598:        UM *in,*out;
                    599:        ML clist;
                    600:
                    601:        n = UDEG(f); mod = blist->mod; np = blist->n;
                    602:        d = 2*n;
                    603:        q = W_UMALLOC(d); wf = W_UMALLOC(d);
                    604:        wm = W_UMALLOC(d); wx = W_UMALLOC(d);
                    605:        wy = W_UMALLOC(d); wu = W_UMALLOC(d);
                    606:        wv = W_UMALLOC(d); wg = W_UMALLOC(d);
                    607:        wa = W_UMALLOC(d); wb = W_UMALLOC(d);
                    608:        ptoum(mod,f,wf); DEG(wg) = 0; COEF(wg)[0] = 1;
                    609:        *clistp = clist = MLALLOC(np); clist->mod = mod; clist->n = np;
                    610:        for ( i = 0, in = (UM *)blist->c, out = (UM *)clist->c; i < np; i++ ) {
                    611:                divum(mod,wf,in[i],q); tum = wf; wf = q; q = tum;
                    612:                cpyum(wf,wx); cpyum(in[i],wy);
                    613:                eucum(mod,wx,wy,wa,wb); mulum(mod,wa,wg,wm);
                    614:                DEG(wm) = divum(mod,wm,in[i],q); out[i] = UMALLOC(DEG(wm));
                    615:                cpyum(wm,out[i]); mulum(mod,q,wf,wu);
                    616:                mulum(mod,wg,wb,wv); addum(mod,wu,wv,wg);
                    617:        }
                    618: }
                    619:
1.4       noro      620: /* find a,b s.t. qa*qg+qb*qh=1 mod q, qg=g mod mod, qh=h mod mod */
                    621: /* q = mod^k */
                    622:
1.6     ! noro      623: void henprep2(int mod,int q,int k,UM f,UM g,UM h,UM qg,UM qh,UM qa,UM qb)
1.4       noro      624: {
                    625:        int n;
                    626:        UM wg,wh,wa,wb;
                    627:        ML bl,cl,bql,cql;
                    628:        P ff;
                    629:
                    630:        n = DEG(f);
                    631:        wg = W_UMALLOC(2*n); wh = W_UMALLOC(2*n);
                    632:        wa = W_UMALLOC(2*n); wb = W_UMALLOC(2*n);
                    633:        cpyum(g,wg); cpyum(h,wh);
                    634:
                    635:        /* wa*g+wb*h = 1 mod mod */
                    636:        eucum(mod,wg,wh,wa,wb);
                    637:
                    638: #if 0
                    639:        /* check */
                    640:        wt = W_UMALLOC(2*n); ws = W_UMALLOC(2*n); wu = W_UMALLOC(2*n);
                    641:        mulum(mod,wa,g,wt);
                    642:        mulum(mod,wb,h,ws);
                    643:        addum(mod,wt,ws,wu);
                    644:        if ( DEG(wu) != 0 || COEF(wu)[0] != 1 )
                    645:                error("henprep 1");
                    646: #endif
                    647:
                    648:        bl = MLALLOC(2); bl->n = 2; bl->mod = mod; bl->c[0] = g; bl->c[1] = h;
                    649:        cl = MLALLOC(2); cl->n = 2; cl->mod = mod; cl->c[0] = wb; cl->c[1] = wa;
                    650:        umtop(CO->v,f,&ff); /* XXX */
                    651:        henprep(ff,bl,cl,&bql,&cql); /* XXX */
                    652:
                    653:        cpyum(bql->c[0],qg); cpyum(bql->c[1],qh);
                    654:        cpyum(cql->c[0],qb); cpyum(cql->c[1],qa);
                    655:
                    656: #if 0
                    657:        /* check */
                    658:        mulum(q,qa,qg,wt);
                    659:        mulum(q,qb,qh,ws);
                    660:        addum(q,wt,ws,wu);
                    661:        if ( DEG(wu) != 0 || COEF(wu)[0] != 1 )
                    662:                error("henprep 2");
                    663: #endif
                    664: }
                    665:
1.1       noro      666: /*
1.4       noro      667:        henprep(f,blist,clist,&bqlist,&cqlist);
1.1       noro      668:  */
                    669:
1.6     ! noro      670: void henprep(P f,ML blist,ML clist,ML *bqlistp,ML *cqlistp)
1.1       noro      671: {
                    672:        register int i,j,k,*px,*py,*pz;
                    673:        int n,pmax,dr,tmp,p,p1,mod,np,b,q;
                    674:        UM w,wm,wn,wa,wt,wq,wf,quot,tum,*in,*inc,*out,*outc;
                    675:        ML bqlist,cqlist;
                    676:
                    677:        n = UDEG(f); p = mod = blist->mod; np = blist->n;
                    678: /*     for ( b = 1, q = mod; q <= (unsigned int)(LBASE / (L)mod); q *= mod, b++ ); */
                    679:        for ( b = 1, q = mod; q <= ((1<<27) / mod); q *= mod, b++ );
                    680:        w = W_UMALLOC(n); ptoum(q,f,w);
                    681:        wm = W_UMALLOC(2*n); wn = W_UMALLOC(2*n);
                    682:        wa = W_UMALLOC(2*n); wt = W_UMALLOC(2*n);
                    683:        wq = W_UMALLOC(2*n); wf = W_UMALLOC(2*n);
                    684:        quot = W_UMALLOC(2*n);
                    685:        *bqlistp = bqlist = MLALLOC(np); *cqlistp = cqlist = MLALLOC(np);
                    686:        for ( i = 0; i < n+2; i++ )
                    687:                COEF(wq)[i] = 0;
                    688:        for ( i = 0,
                    689:                  in = (UM *)blist->c, inc = (UM *)clist->c,
                    690:                  out = (UM *)bqlist->c, outc = (UM *)cqlist->c;
                    691:                  i < np; i++ ) {
                    692:                out[i] = C_UMALLOC(n+1); cpyum(in[i],out[i]);
                    693:                outc[i] = C_UMALLOC(n+1); cpyum(inc[i],outc[i]);
                    694:        }
                    695:        for ( pmax = 1, i = b; i > 0; i-- )
                    696:                pmax *= mod;
                    697:        for ( i = 1; i < b; i++, p = p1 ) {
                    698:                cpyum(out[0],wm);
                    699:                for ( j = 1; j < np; j++ ) {
                    700:                        mulum(pmax,wm,out[j],wn);
                    701:                        tum = wm; wm = wn; wn = tum;
                    702:                }
                    703:                for ( j = n, px = COEF(w), py = COEF(wm), pz = COEF(wt); j >= 0; j-- ) {
                    704:                        tmp = ( ( px[j] - py[j] ) / p ) % mod;
                    705:                        pz[j] = ( tmp >= 0? tmp : tmp + mod );
                    706:                }
                    707:                degum(wt,n);
                    708:                for ( j = 1; j < np; j++ ) {
                    709:                        mulum(mod,wt,inc[j],wm); dr = divum(mod,wm,in[j],quot);
                    710:                        for ( k = DEG(quot); k >= 0; k-- )
                    711:                                COEF(wq)[k] = ( COEF(wq)[k] + COEF(quot)[k] ) % mod;
                    712:                        for ( k = dr, px = COEF(out[j]), py = COEF(wm); k >= 0; k-- )
                    713:                                px[k] += p * py[k];
                    714:                }
                    715:                degum(wq,n); mulum(mod,wq,in[0],wm);
                    716:                mulum(mod,wt,inc[0],wn); addum(mod,wm,wn,wa);
                    717:                for ( j = DEG(wa), px = COEF(out[0]), py = COEF(wa); j >= 0; j-- )
                    718:                        px[j] += p * py[j];
                    719:                for ( j = n, px = COEF(wq); j >= 0; j-- )
                    720:                        px[j] = 0;
                    721:                p1 = p * mod;
                    722:                for ( j = n, px = COEF(wt); j >= 1; j-- )
                    723:                        px[j] = 0;
                    724:                px[0] = 1;
                    725:                for ( j = 0; j < np; j++ ) {
                    726:                        cpyum(w,wf);
                    727:                        for ( k = DEG(wf), px = COEF(wf); k >= 0; k-- )
                    728:                                px[k] %= p1;
                    729:                        divum(p1,wf,out[j],quot); mulum(p1,outc[j],quot,wm);
                    730:                        for ( k = DEG(wm), px = COEF(wt), py = COEF(wm); k >= 0; k-- )
                    731:                                px[k] = ( px[k] - py[k] ) % p1;
                    732:                }
                    733:                degum(wt,n);
                    734:                for ( j = DEG(wt), px = COEF(wt); j >= 0; j-- )
                    735:                        px[j] = ((tmp=(px[j]/p)%mod)>= 0?tmp:tmp + mod);
                    736:                for ( j = 0; j < np; j++ ) {
                    737:                        mulum(mod,wt,outc[j],wm); dr = divum(mod,wm,in[j],quot);
                    738:                        for ( k = dr, px = COEF(outc[j]), py = COEF(wm); k >= 0; k-- )
                    739:                                px[k] += p * py[k];
                    740:                        degum(outc[j],MAX(DEG(outc[j]),dr));
                    741:                }
                    742:        }
                    743:        bqlist->n = cqlist->n = np;
                    744:        bqlist->mod = cqlist->mod = q;
                    745: }
                    746:
                    747: /*
                    748:        henmain(fl,bqlist,cqlist,listp)
                    749: */
                    750:
1.6     ! noro      751: void henmain(LUM f,ML bqlist,ML cqlist,ML *listp)
1.1       noro      752: {
                    753:        register int i,j,k;
                    754:        int *px,*py;
                    755:        int **pp,**pp1;
                    756:        int n,np,mod,bound,dr,tmp;
                    757:        UM wt,wq0,wq,wr,wm,wm0,wa,q;
                    758:        LUM wb0,wb1,tlum;
                    759:        UM *b,*c;
                    760:        LUM *l;
                    761:        ML list;
                    762:
                    763:        n = DEG(f); np = bqlist->n; mod = bqlist->mod; bound = bqlist->bound;
                    764:        *listp = list = MLALLOC(n);
                    765:        list->n = np; list->mod = mod; list->bound = bound;
                    766:        W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);
                    767:        wt = W_UMALLOC(n); wq0 = W_UMALLOC(n); wq = W_UMALLOC(n);
                    768:        wr = W_UMALLOC(n); wm = W_UMALLOC(2*n); wm0 = W_UMALLOC(2*n);
                    769:        wa = W_UMALLOC(2*n); q = W_UMALLOC(2*n);
                    770:        b = (UM *)bqlist->c; c = (UM *)cqlist->c; l = (LUM *)list->c;
                    771:        for ( i = 0; i < np; i++ ) {
                    772:                l[i] = LUMALLOC(DEG(b[i]),bound);
                    773:                for ( j = DEG(b[i]), pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- )
                    774:                        pp[j][0] = px[j];
                    775:        }
1.5       noro      776: #if 0
1.4       noro      777:        fprintf(stderr,"bound=%d\n",bound);
1.5       noro      778: #endif
1.1       noro      779:        for ( i = 1; i < bound; i++ ) {
1.5       noro      780: #if 0
1.4       noro      781:                fprintf(stderr,".");
1.5       noro      782: #endif
1.1       noro      783:                mullum(mod,i+1,l[0],l[1],wb0);
                    784:                for ( j = 2; j < np; j++ ) {
                    785:                        mullum(mod,i+1,l[j],wb0,wb1);
                    786:                        tlum = wb0; wb0 = wb1; wb1 = tlum;
                    787:                }
                    788:                for ( j = n, px = COEF(wt); j >= 0; j-- )
                    789:                        px[j] = 0;
                    790:                for ( j = n, pp = COEF(f), pp1 = COEF(wb0); j >= 0; j-- ) {
                    791:                        tmp = ( pp[j][i] - pp1[j][i] ) % mod;
                    792:                        COEF(wt)[j] = ( tmp < 0 ? tmp + mod : tmp );
                    793:                }
                    794:                degum(wt,n);
                    795:                for ( j = n, px = COEF(wq0); j >= 0; j-- )
                    796:                        px[j] = 0;
                    797:                for ( j = 1; j < np; j++ ) {
                    798:                        mulum(mod,wt,c[j],wm); dr = divum(mod,wm,b[j],q);
                    799:                        for ( k = DEG(q), px = COEF(wq0), py = COEF(q); k >= 0; k-- )
                    800:                                px[k] = ( px[k] + py[k] ) % mod;
                    801:                        for ( k = dr, pp = COEF(l[j]), px = COEF(wm); k >= 0; k-- )
                    802:                                pp[k][i] = px[k];
                    803:                }
                    804:                degum(wq0,n); mulum(mod,wq0,b[0],wm);
                    805:                mulum(mod,wt,c[0],wm0); addum(mod,wm,wm0,wa);
                    806:                for ( j = DEG(wa), pp = COEF(l[0]), px = COEF(wa); j >= 0; j-- )
                    807:                        pp[j][i] = px[j];
                    808:                for ( j = n, px = COEF(wq0); j >= 0; j-- )
                    809:                        px[j] = 0;
                    810:        }
1.5       noro      811: #if 0
1.4       noro      812:        fprintf(stderr,"\n");
1.5       noro      813: #endif
1.1       noro      814: }
                    815:
                    816: static double M;
                    817: static int E;
                    818:
1.6     ! noro      819: int mignotte(int q,P f)
1.1       noro      820: {
                    821:        int p;
                    822:        unsigned int *b;
                    823:        N c;
                    824:        DCP dc;
                    825:
                    826:        for ( dc = DC(f), M = 0, E = 0; dc; dc = NEXT(dc) ) {
                    827:                c = NM((Q)COEF(dc)); p = PL(c); b = BD(c);
                    828:                sqad(b[p-1],(p-1)*BSH);
                    829:        }
                    830:        if (E % 2) M *= 2; M = ceil(sqrt(M)); E /= 2;
                    831:        c = NM((Q)COEF(DC(f))); p = PL(c); M *= ((double)BD(c)[p-1]+1.0); E += (p-1) * BSH;
                    832:        return (int)ceil( (0.31*(E+UDEG(f)+1)+log10((double)M)) / log10((double)q) );
                    833: }
                    834:
1.6     ! noro      835: int mig(int q,int d,P f)
1.1       noro      836: {
                    837:        int p;
                    838:        unsigned int *b;
                    839:        N c;
                    840:        DCP dc;
                    841:
                    842:        for ( dc = DC(f), M = 0, E = 0; dc; dc = NEXT(dc) ) {
                    843:                c = NM((Q)COEF(dc)); p = PL(c); b = BD(c);
                    844:                sqad(b[p-1],(p-1)*BSH);
                    845:        }
                    846:        if (E % 2) M *= 2; M = ceil(sqrt(M)); E /= 2;
                    847:        c = NM((Q)COEF(DC(f))); p = PL(c);
                    848:        M *= (BD(c)[p-1]+1); E += (p-1) * BSH;
                    849:        return (int)ceil( (0.31*(E+d+1)+log10((double)M)) / log10((double)q) );
                    850: }
                    851:
1.6     ! noro      852: void sqad(unsigned int man,int exp)
1.1       noro      853: {
                    854:        int e,sqe;
                    855:        unsigned int t;
                    856:        double man1,d,sqm;
                    857:        int diff;
                    858:
                    859:        if ( man == BMASK ) {
                    860:                e = BSH; man1 = 1.0;
                    861:        } else {
                    862:                man += 1;
                    863:                for ( e = 0, t = man; t; e++, t >>= 1 );
                    864:                e--; d = (double)(1<<e);
                    865:                man1 = ((double)man)/d;
                    866:        }
                    867:        exp += e; sqm = man1 * man1; sqe = 2 * exp;
                    868:        if ( sqm >= 2.0 ) {
                    869:                sqm /= 2.0; sqe++;
                    870:        }
                    871:        diff = E - sqe;
                    872:        if ( diff > 18 )
                    873:                return;
                    874:        if ( diff < -18 ) {
                    875:                M = sqm; E = sqe;
                    876:                return;
                    877:        }
                    878:        if ( diff >= 0 )
                    879:                M += (sqm / (double)(1<<diff));
                    880:        else {
                    881:                M = ( ( M / (double)(1<<-diff)) + sqm ); E = sqe;
                    882:        }
                    883:        if ( M >= 2.0 ) {
                    884:                M /= 2.0; E++;
                    885:        }
                    886: }
                    887:
1.6     ! noro      888: void ptolum(int q,int bound,P f,LUM fl)
1.1       noro      889: {
                    890:        DCP dc;
                    891:        int i,j;
                    892:        int **pp;
                    893:        int d,br,s;
                    894:        unsigned int r;
                    895:        int *c;
                    896:        unsigned int *m,*w;
                    897:
                    898:        for ( dc = DC(f), pp = COEF(fl); dc; dc = NEXT(dc) ) {
                    899:                d = PL(NM((Q)COEF(dc))); m = BD(NM((Q)COEF(dc)));
                    900:                c = pp[QTOS(DEG(dc))]; w = (unsigned int *)W_ALLOC(d);
                    901:                for ( i = 0; i < d; i++ )
                    902:                        w[i] = m[i];
                    903:                for ( i = 0; d >= 1; ) {
                    904:                        for ( j = d - 1, r = 0; j >= 0; j-- ) {
                    905:                                DSAB(q,r,w[j],w[j],r)
                    906:                        }
                    907:                        c[i++] = (int)r;
                    908:                        if ( !w[d-1] )
                    909:                                d--;
                    910:                }
                    911:                if ( SGN((Q)COEF(dc)) < 0 )
                    912:                        for (i = 0, br = 0; i < bound; i++ )
                    913:                                if ( ( s = -(c[i] + br) ) < 0 ) {
                    914:                                        c[i] = s + q; br = 1;
                    915:                                } else {
                    916:                                        c[i] = 0; br = 0;
                    917:                                }
                    918:        }
                    919: }
                    920:
1.6     ! noro      921: void modfctrp(P p,int mod,int flag,DCP *dcp)
1.1       noro      922: {
                    923:        int cm,n,i,j,k;
                    924:        DCP dc,dc0;
                    925:        P zp;
                    926:        Q c,q;
                    927:        UM mp;
                    928:        UM *tl;
                    929:        struct oDUM *udc,*udc1;
                    930:
                    931:        if ( !p ) {
                    932:                *dcp = 0; return;
                    933:        }
                    934:        ptozp(p,1,&c,&zp);
                    935:        if ( DN(c) || !(cm = rem(NM(c),mod)) ) {
                    936:                *dcp = 0; return;
                    937:        }
                    938:        mp = W_UMALLOC(UDEG(p));
                    939:        ptoum(mod,zp,mp);
                    940:        if ( (n = DEG(mp)) < 0 ) {
                    941:                *dcp = 0; return;
                    942:        } else if ( n == 0 ) {
                    943:                cm = dmar(cm,COEF(mp)[0],0,mod); STOQ(cm,q);
                    944:                NEWDC(dc); COEF(dc) = (P)q; DEG(dc) = ONE;
                    945:                NEXT(dc) = 0; *dcp = dc;
                    946:                return;
                    947:        }
                    948:        if ( COEF(mp)[n] != 1 ) {
                    949:                cm = dmar(cm,COEF(mp)[n],0,mod);
                    950:                i = invm(COEF(mp)[n],mod);
                    951:                for ( j = 0; j <= n; j++ )
                    952:                        COEF(mp)[j] = dmar(COEF(mp)[j],i,0,mod);
                    953:        }
                    954:        W_CALLOC(n+1,struct oDUM,udc);
                    955:        gensqfrum(mod,mp,udc);
                    956:        switch ( flag ) {
                    957:                case FCTR:
                    958:                        tl = (UM *)ALLOCA((n+1)*sizeof(UM));
                    959:                        W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
                    960:                        for ( i = 0,j = 0; udc[i].f; i++ )
                    961:                                if ( DEG(udc[i].f) == 1 ) {
                    962:                                        udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
                    963:                                } else {
                    964:                                        bzero((char *)tl,(n+1)*sizeof(UM));
                    965:                                        berlemain(mod,udc[i].f,tl);
                    966:                                        for ( k = 0; tl[k]; k++, j++ ) {
                    967:                                                udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
                    968:                                        }
                    969:                                }
                    970:                        udc = udc1; break;
                    971:                case SQFR:
                    972:                        break;
                    973:                case DDD:
                    974:                        tl = (UM *)ALLOCA((n+1)*sizeof(UM));
                    975:                        W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
                    976:                        for ( i = 0,j = 0; udc[i].f; i++ )
                    977:                                if ( DEG(udc[i].f) == 1 ) {
                    978:                                        udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
                    979:                                } else {
                    980:                                        bzero((char *)tl,(n+1)*sizeof(UM));
                    981:                                        ddd(mod,udc[i].f,tl);
                    982:                                        for ( k = 0; tl[k]; k++, j++ ) {
                    983:                                                udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
                    984:                                        }
                    985:                                }
                    986:                        udc = udc1; break;
                    987:                case NEWDDD:
                    988:                        tl = (UM *)ALLOCA((n+1)*sizeof(UM));
                    989:                        W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
                    990:                        for ( i = 0,j = 0; udc[i].f; i++ )
                    991:                                if ( DEG(udc[i].f) == 1 ) {
                    992:                                        udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
                    993:                                } else {
                    994:                                        bzero((char *)tl,(n+1)*sizeof(UM));
                    995:                                        if ( mod == 2 )
                    996:                                                berlemain(mod,udc[i].f,tl);
                    997:                                        else
                    998:                                                newddd(mod,udc[i].f,tl);
                    999:                                        for ( k = 0; tl[k]; k++, j++ ) {
                   1000:                                                udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
                   1001:                                        }
                   1002:                                }
                   1003:                        udc = udc1; break;
                   1004:        }
                   1005:        NEWDC(dc0); STOQ(cm,q); COEF(dc0) = (P)q; DEG(dc0) = ONE; dc = dc0;
                   1006:        for ( n = 0; udc[n].f; n++ ) {
                   1007:                NEWDC(NEXT(dc)); dc = NEXT(dc);
                   1008:                STOQ(udc[n].n,DEG(dc)); umtop(VR(p),udc[n].f,&COEF(dc));
                   1009:        }
                   1010:        NEXT(dc) = 0; *dcp = dc0;
                   1011: }
                   1012:
1.6     ! noro     1013: void gensqfrum(int mod,UM p,struct oDUM *dc)
1.1       noro     1014: {
                   1015:        int n,i,j,d;
                   1016:        UM t,s,g,f,f1,b;
                   1017:
                   1018:        if ( (n = DEG(p)) == 1 ) {
                   1019:                dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
                   1020:                return;
                   1021:        }
                   1022:        t = W_UMALLOC(n); s = W_UMALLOC(n); g = W_UMALLOC(n);
                   1023:        f = W_UMALLOC(n); f1 = W_UMALLOC(n); b = W_UMALLOC(n);
                   1024:        diffum(mod,p,t); cpyum(p,s); Gcdum(mod,t,s,g);
                   1025:        if ( !DEG(g) ) {
                   1026:                dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
                   1027:                return;
                   1028:        }
                   1029:        cpyum(p,b); cpyum(p,t); Divum(mod,t,g,f);
                   1030:        for ( i = 0, d = 0; DEG(f); i++ ) {
                   1031:                while ( 1 ) {
                   1032:                        cpyum(b,t);
                   1033:                        if ( Divum(mod,t,f,s) >= 0 )
                   1034:                                break;
                   1035:                        else {
                   1036:                                cpyum(s,b); d++;
                   1037:                        }
                   1038:                }
                   1039:                cpyum(b,t); cpyum(f,s); Gcdum(mod,t,s,f1);
                   1040:                Divum(mod,f,f1,s); cpyum(f1,f);
                   1041:                dc[i].f = UMALLOC(DEG(s)); cpyum(s,dc[i].f); dc[i].n = d;
                   1042:        }
                   1043:        if ( DEG(b) > 0 ) {
                   1044:                d = 1;
                   1045:                while ( 1 ) {
                   1046:                        cpyum(b,t);
                   1047:                        for ( j = DEG(t); j >= 0; j-- )
                   1048:                                if ( COEF(t)[j] && (j % mod) )
                   1049:                                        break;
                   1050:                        if ( j >= 0 )
                   1051:                                break;
                   1052:                        else {
                   1053:                                DEG(s) = DEG(t)/mod;
                   1054:                                for ( j = 0; j <= DEG(t); j++ )
                   1055:                                        COEF(s)[j] = COEF(t)[j*mod];
                   1056:                                cpyum(s,b); d *= mod;
                   1057:                        }
                   1058:                }
                   1059:                gensqfrum(mod,b,dc+i);
                   1060:                for ( j = i; dc[j].f; j++ )
                   1061:                        dc[j].n *= d;
                   1062:        }
                   1063: }
                   1064:
                   1065: #if 0
1.6     ! noro     1066: void srchum(int mod,UM p1,UM p2,UM gr)
1.1       noro     1067: {
                   1068:        UM m,m1,m2,q,r,t,g1,g2;
                   1069:        int lc,d,d1,d2,i,j,k,l,l1,l2,l3,tmp,adj;
                   1070:        V v;
                   1071:
                   1072:        d = MAX(DEG(p1),DEG(p2));
                   1073:        g1 = W_UMALLOC(d); g2 = W_UMALLOC(d);
                   1074:        bzero((char *)g1,(d+2)*sizeof(int)); bzero((char *)g2,(d+2)*sizeof(int));
                   1075:        if ( d == DEG(p1) ) {
                   1076:                cpyum(p1,g1); cpyum(p2,g2);
                   1077:        } else {
                   1078:                cpyum(p1,g2); cpyum(p2,g1);
                   1079:        }
                   1080:        if ( ( d1 = DEG(g1) ) > ( d2 = DEG(g2) ) ) {
                   1081:                j = d1 - 1; adj = 1;
                   1082:        } else
                   1083:                j = d2;
                   1084:        lc = 1;
                   1085:        r = W_UMALLOC(d1+d2); q = W_UMALLOC(d1+d2);
                   1086:        m1 = W_UMALLOC(d1+d2); t = W_UMALLOC(d1+d2);
                   1087:        bzero((char *)r,(d1+d2+2)*sizeof(int)); bzero((char *)q,(d1+d2+2)*sizeof(int));
                   1088:        bzero((char *)m1,(d1+d2+2)*sizeof(int)); bzero((char *)t,(d1+d2+2)*sizeof(int));
                   1089:        m = W_UMALLOC(0); bzero((char *)m,2*sizeof(int));
                   1090:        adj = pwrm(mod,COEF(g2)[DEG(g2)],DEG(g1));
                   1091:        DEG(m) = 0; COEF(m)[0] = invm(COEF(g2)[DEG(g2)],mod);
                   1092:        Mulum(mod,g2,m,r); cpyum(r,g2);
                   1093:        while ( 1 ) {
                   1094:                if ( ( k = DEG(g2) ) < 0 ) {
                   1095:                        DEG(gr) = -1;
                   1096:                        return;
                   1097:                }
                   1098:                if ( k == j ) {
                   1099:                        if ( k == 0 ) {
                   1100:                                DEG(m) = 0; COEF(m)[0] = adj;
                   1101:                                Mulum(mod,g2,m,gr);
                   1102:                                return;
                   1103:                        } else {
                   1104:                                DEG(m) = 0;
                   1105:                                COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
                   1106:                                Mulum(mod,g1,m,r); DEG(r) = Divum(mod,r,g2,t);
                   1107:                                DEG(m) = 0; COEF(m)[0] = dmb(mod,lc,lc,&tmp);
                   1108:                                Divum(mod,r,m,q); cpyum(g2,g1); cpyum(q,g2);
                   1109:                                lc = COEF(g1)[DEG(g1)]; j = k - 1;
                   1110:                        }
                   1111:                } else {
                   1112:                        d = j - k;
                   1113:                        DEG(m) = 0; COEF(m)[0] = pwrm(mod,COEF(g2)[DEG(g2)],d);
                   1114:                        Mulum(mod,g2,m,m1); l = pwrm(mod,lc,d);
                   1115:                        DEG(m) = 0; COEF(m)[0] = l; Divum(mod,m1,m,t);
                   1116:                        if ( k == 0 ) {
                   1117:                                DEG(m) = 0; COEF(m)[0] = adj;
                   1118:                                Mulum(mod,t,m,gr);
                   1119:                                return;
                   1120:                        } else {
                   1121:                                DEG(m) = 0;
                   1122:                                COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
                   1123:                                Mulum(mod,g1,m,r); DEG(r) = Divum(mod,r,g2,q);
                   1124:                                l1 = dmb(mod,lc,lc,&tmp); l2 = dmb(mod,l,l1,&tmp);
                   1125:                                DEG(m) = 0; COEF(m)[0] = l2;
                   1126:                                Divum(mod,r,m,q); cpyum(t,g1); cpyum(q,g2);
                   1127:                                if ( d % 2 )
                   1128:                                        for ( i = DEG(g2); i >= 0; i-- )
                   1129:                                                COEF(g2)[i] = ( mod - COEF(g2)[i] ) % mod;
                   1130:                                lc = COEF(g1)[DEG(g1)]; j = k - 1;
                   1131:                        }
                   1132:                }
                   1133:        }
                   1134: }
                   1135:
1.6     ! noro     1136: UM *resberle(int mod,UM f,UM *fp)
1.1       noro     1137: {
                   1138:        UM w,wg,ws,wf,f0,gcd,q,res;
                   1139:        int n;
                   1140:        register int i;
                   1141:
                   1142:        n = DEG(f); wg = W_UMALLOC(n); mini(mod,f,wg);
                   1143:        if ( DEG(wg) <= 0 ) {
                   1144:                f0 = UMALLOC(n); cpyum(f,f0); *fp++ = f0;
                   1145:                return ( fp );
                   1146:        }
                   1147:        f0 = W_UMALLOC(n); cpyum(f,f0);
                   1148:        ws = W_UMALLOC(n); wf = W_UMALLOC(n);
                   1149:        q = W_UMALLOC(n); gcd = W_UMALLOC(n);
                   1150:        res = W_UMALLOC(2*n);
                   1151:        srchum(mod,f,wg,res);
                   1152:        for ( i = 0; i < mod; i++ ) {
                   1153:                if ( substum(mod,res,i) )
                   1154:                        continue;
                   1155:                cpyum(f0,wf); cpyum(wg,ws);
                   1156:                COEF(ws)[0] = ( COEF(ws)[0] + mod - i ) % mod;
                   1157:                Gcdum(mod,wf,ws,gcd);
                   1158:                if ( DEG(gcd) > 0 ) {
                   1159:                        if ( DEG(gcd) < n ) {
                   1160:                                Divum(mod,f0,gcd,q); f0 = q; fp = resberle(mod,gcd,fp);
                   1161:                        }
                   1162:                        break;
                   1163:                }
                   1164:        }
                   1165:        fp = resberle(mod,f0,fp);
                   1166:        return ( fp );
                   1167: }
                   1168:
1.6     ! noro     1169: int substum(int mod,UM p,int a)
1.1       noro     1170: {
                   1171:        int i,j,s;
                   1172:        int *c;
                   1173:
                   1174:        if ( DEG(p) < 0 )
                   1175:                return 0;
                   1176:        if ( DEG(p) == 0 )
                   1177:                return COEF(p)[0];
                   1178:        for ( i = DEG(p), c = COEF(p), s = c[i]; i >= 0; ) {
                   1179:                for ( j = i--; (i>=0) && !c[i]; i-- );
                   1180:                if ( i >= 0 )
                   1181:                        s = (s*pwrm(mod,a,j-i)%mod+c[i])%mod;
                   1182:                else
                   1183:                        s = s*pwrm(mod,a,j)%mod;
                   1184:        }
                   1185:        return s;
                   1186: }
                   1187: #endif
                   1188:
1.6     ! noro     1189: void ddd(int mod,UM f,UM *r)
1.1       noro     1190: {
                   1191:        register int i,j;
                   1192:        int d,n;
                   1193:        UM q,s,t,u,v,w,g,x,m;
                   1194:        UM *base;
                   1195:
                   1196:        n = DEG(f);
                   1197:        if ( n == 1 ) {
                   1198:                r[0] = UMALLOC(1); cpyum(f,r[0]); r[1] = 0; return;
                   1199:        }
                   1200:        base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
                   1201:        w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
                   1202:        base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
                   1203:        t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
                   1204:        pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
                   1205:        for ( i = 2; i < n; i++ ) {
                   1206:                mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
                   1207:                base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
                   1208:        }
                   1209:        v = W_UMALLOC(n); cpyum(f,v);
                   1210:        DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
                   1211:        x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
                   1212:        t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
                   1213:        for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
                   1214:                for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
                   1215:                        if ( COEF(w)[i] ) {
                   1216:                                Mulsum(mod,base[i],COEF(w)[i],s);
                   1217:                                addum(mod,s,t,u); cpyum(u,t);
                   1218:                        }
                   1219:                cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
                   1220:                if ( DEG(g) >= 1 ) {
                   1221:                        canzas(mod,g,d,base,r+j); j += DEG(g)/d;
                   1222:                        Divum(mod,v,g,q); cpyum(q,v);
                   1223:                        DEG(w) = Divum(mod,w,v,q);
                   1224:                        for ( i = 0; i < DEG(v); i++ )
                   1225:                                DEG(base[i]) = Divum(mod,base[i],v,q);
                   1226:                }
                   1227:        }
                   1228:        if ( DEG(v) ) {
                   1229:                r[j] = UMALLOC(DEG(v)); cpyum(v,r[j]); j++;
                   1230:        }
                   1231:        r[j] = 0;
                   1232: }
                   1233:
                   1234: #if 0
1.6     ! noro     1235: void canzas(int mod,UM f,int d,UM *base,UM *r)
1.1       noro     1236: {
                   1237:        UM t,s,u,w,g,o,q;
                   1238:        N n1,n2,n3,n4,n5;
                   1239:        UM *b;
                   1240:        int n,m,i;
                   1241:
                   1242:        if ( DEG(f) == d ) {
                   1243:                r[0] = UMALLOC(d); cpyum(f,r[0]);
                   1244:                return;
                   1245:        } else {
                   1246:                n = DEG(f); b = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)b,n*sizeof(UM));
                   1247:                for ( i = 0, m = 0; i < n; i++ )
                   1248:                        m = MAX(m,DEG(base[i]));
                   1249:                q = W_UMALLOC(m);
                   1250:                for ( i = 0; i < n; i++ ) {
                   1251:                        b[i] = W_UMALLOC(DEG(base[i])); cpyum(base[i],b[i]);
                   1252:                        DEG(b[i]) = Divum(mod,b[i],f,q);
                   1253:                }
                   1254:                t = W_UMALLOC(2*d);
                   1255:                s = W_UMALLOC(DEG(f)); u = W_UMALLOC(DEG(f));
                   1256:                w = W_UMALLOC(DEG(f)); g = W_UMALLOC(DEG(f));
                   1257:                o = W_UMALLOC(0); DEG(o) = 0; COEF(o)[0] = 1;
                   1258:                STON(mod,n1); pwrn(n1,d,&n2); subn(n2,ONEN,&n3);
                   1259:                STON(2,n4); divsn(n3,n4,&n5);
                   1260:                while ( 1 ) {
                   1261:                        randum(mod,2*d,t); spwrum(mod,f,b,t,n5,s);
                   1262:                        subum(mod,s,o,u); cpyum(f,w); Gcdum(mod,w,u,g);
                   1263:                        if ( (DEG(g) >= 1) && (DEG(g) < DEG(f)) ) {
                   1264:                                canzas(mod,g,d,b,r);
                   1265:                                cpyum(f,w); Divum(mod,w,g,s);
                   1266:                                canzas(mod,s,d,b,r+DEG(g)/d);
                   1267:                                return;
                   1268:                        }
                   1269:                }
                   1270:        }
                   1271: }
                   1272: #else
1.6     ! noro     1273: void canzas(int mod,UM f,int d,UM *base,UM *r)
1.1       noro     1274: {
                   1275:        UM t,s,u,w,g,o,q;
                   1276:        N n1,n2,n3,n4,n5;
                   1277:        UM *b;
                   1278:        int n,m,i;
                   1279:
                   1280:        if ( DEG(f) == d ) {
                   1281:                r[0] = UMALLOC(d); cpyum(f,r[0]);
                   1282:                return;
                   1283:        } else {
                   1284:                n = DEG(f); b = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)b,n*sizeof(UM));
                   1285:                for ( i = 0, m = 0; i < n; i++ )
                   1286:                        m = MAX(m,DEG(base[i]));
                   1287:                q = W_UMALLOC(m);
                   1288:                for ( i = 0; i < n; i++ ) {
                   1289:                        b[i] = W_UMALLOC(DEG(base[i])); cpyum(base[i],b[i]);
                   1290:                        DEG(b[i]) = Divum(mod,b[i],f,q);
                   1291:                }
                   1292:                t = W_UMALLOC(2*d);
                   1293:                s = W_UMALLOC(DEG(f)); u = W_UMALLOC(DEG(f));
                   1294:                w = W_UMALLOC(DEG(f)); g = W_UMALLOC(DEG(f));
                   1295:                o = W_UMALLOC(0); DEG(o) = 0; COEF(o)[0] = 1;
                   1296:                STON(mod,n1); pwrn(n1,d,&n2); subn(n2,ONEN,&n3);
                   1297:                STON(2,n4); divsn(n3,n4,&n5);
                   1298:                while ( 1 ) {
                   1299:                        randum(mod,2*d,t); spwrum0(mod,f,t,n5,s);
                   1300:                        subum(mod,s,o,u); cpyum(f,w); Gcdum(mod,w,u,g);
                   1301:                        if ( (DEG(g) >= 1) && (DEG(g) < DEG(f)) ) {
                   1302:                                canzas(mod,g,d,b,r);
                   1303:                                cpyum(f,w); Divum(mod,w,g,s);
                   1304:                                canzas(mod,s,d,b,r+DEG(g)/d);
                   1305:                                return;
                   1306:                        }
                   1307:                }
                   1308:        }
                   1309: }
                   1310: #endif
                   1311:
1.6     ! noro     1312: void randum(int mod,int d,UM p)
1.1       noro     1313: {
                   1314:        unsigned int n;
                   1315:        int i;
                   1316:
                   1317:        n = ((unsigned int)random()) % d; DEG(p) = n; COEF(p)[n] = 1;
                   1318:        for ( i = 0; i < (int)n; i++ )
                   1319:                COEF(p)[i] = ((unsigned int)random()) % mod;
                   1320: }
                   1321:
1.6     ! noro     1322: void pwrmodum(int mod,UM p,int e,UM f,UM pr)
1.1       noro     1323: {
                   1324:        UM wt,ws,q;
                   1325:
                   1326:        if ( e == 0 ) {
                   1327:                DEG(pr) = 0; COEF(pr)[0] = 1;
                   1328:        } else if ( DEG(p) < 0 )
                   1329:                DEG(pr) = -1;
                   1330:        else if ( e == 1 ) {
                   1331:                q = W_UMALLOC(DEG(p)); cpyum(p,pr);
                   1332:                DEG(pr) = divum(mod,pr,f,q);
                   1333:        } else if ( DEG(p) == 0 ) {
                   1334:                DEG(pr) = 0; COEF(pr)[0] = pwrm(mod,COEF(p)[0],e);
                   1335:        } else {
                   1336:                wt = W_UMALLOC(2*DEG(f)); ws = W_UMALLOC(2*DEG(f));
                   1337:                q = W_UMALLOC(2*DEG(f));
                   1338:                pwrmodum(mod,p,e/2,f,wt);
                   1339:                if ( !(e%2) )  {
                   1340:                        mulum(mod,wt,wt,pr); DEG(pr) = divum(mod,pr,f,q);
                   1341:                } else {
                   1342:                        mulum(mod,wt,wt,ws); DEG(ws) = divum(mod,ws,f,q);
                   1343:                        mulum(mod,ws,p,pr); DEG(pr) = divum(mod,pr,f,q);
                   1344:                }
                   1345:        }
                   1346: }
                   1347:
1.6     ! noro     1348: void spwrum(int mod,UM m,UM *base,UM f,N e,UM r)
1.1       noro     1349: {
                   1350:        int a,n,i;
                   1351:        N e1,an;
                   1352:        UM t,s,u,q,r1,r2;
                   1353:
                   1354:        if ( !e ) {
                   1355:                DEG(r) = 0; COEF(r)[0] = 1;
                   1356:        } else if ( UNIN(e) )
                   1357:                cpyum(f,r);
                   1358:        else if ( (PL(e) == 1) && (BD(e)[0] < (unsigned int)mod) )
                   1359:                spwrum0(mod,m,f,e,r);
                   1360:        else {
                   1361:                a = divin(e,mod,&e1); STON(a,an);
                   1362:                n = DEG(m); t = W_UMALLOC(n); s = W_UMALLOC(n);
                   1363:                u = W_UMALLOC(2*n); q = W_UMALLOC(2*n);
                   1364:                for ( DEG(t) = -1, i = 0; i <= DEG(f); i++ )
                   1365:                        if ( COEF(f)[i] ) {
                   1366:                                Mulsum(mod,base[i],COEF(f)[i],s);
                   1367:                                addum(mod,s,t,u); cpyum(u,t);
                   1368:                        }
                   1369:                r1 = W_UMALLOC(n); spwrum0(mod,m,f,an,r1);
                   1370:                r2 = W_UMALLOC(n); spwrum(mod,m,base,t,e1,r2);
                   1371:                Mulum(mod,r1,r2,u); DEG(u) = Divum(mod,u,m,q);
                   1372:                cpyum(u,r);
                   1373:        }
                   1374: }
                   1375:
1.6     ! noro     1376: void spwrum0(int mod,UM m,UM f,N e,UM r)
1.1       noro     1377: {
                   1378:        UM t,s,q;
                   1379:        N e1;
                   1380:        int a;
                   1381:
                   1382:        if ( !e ) {
                   1383:                DEG(r) = 0; COEF(r)[0] = 1;
                   1384:        } else if ( UNIN(e) )
                   1385:                cpyum(f,r);
                   1386:        else {
                   1387:                a = divin(e,2,&e1);
                   1388:                t = W_UMALLOC(2*DEG(m)); spwrum0(mod,m,f,e1,t);
                   1389:                s = W_UMALLOC(2*DEG(m)); q = W_UMALLOC(2*DEG(m));
                   1390:                Mulum(mod,t,t,s); DEG(s) = Divum(mod,s,m,q);
                   1391:                if ( a ) {
                   1392:                        Mulum(mod,s,f,t); DEG(t) = Divum(mod,t,m,q); cpyum(t,r);
                   1393:         } else
                   1394:                        cpyum(s,r);
                   1395:        }
                   1396: }
                   1397:
                   1398: #if 0
1.6     ! noro     1399: void Mulum(int mod,UM p1,UM p2,UM pr)
1.1       noro     1400: {
                   1401:        register int *pc1,*pcr;
                   1402:        register int mul,i,j,d1,d2;
                   1403:        int *c1,*c2,*cr;
                   1404:
                   1405:        if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
                   1406:                DEG(pr) = -1;
                   1407:                return;
                   1408:        }
                   1409:        c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
                   1410:        bzero((char *)cr,(d1+d2+1)*sizeof(int));
                   1411:        for ( i = 0; i <= d2; i++, cr++ )
                   1412:                if ( mul = *c2++ )
                   1413:                        for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ )
                   1414:                                *pcr = (*pc1 * mul + *pcr) % mod;
                   1415:        DEG(pr) = d1 + d2;
                   1416: }
                   1417:
1.6     ! noro     1418: void Mulsum(int mod,UM p,int n,UM pr)
1.1       noro     1419: {
                   1420:        register int *sp,*dp;
                   1421:        register int i;
                   1422:
                   1423:        for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
                   1424:                  i >= 0; i--, dp--, sp-- )
                   1425:                *dp = (*sp * n) % mod;
                   1426: }
                   1427:
1.6     ! noro     1428: int Divum(int mod,UM p1,UM p2,UM pq)
1.1       noro     1429: {
                   1430:        register int *pc1,*pct;
                   1431:        register int tmp,i,j,inv;
                   1432:        int *c1,*c2,*ct;
                   1433:        int d1,d2,dd,hd;
                   1434:
                   1435:        if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
                   1436:                DEG(pq) = -1;
                   1437:                return( d1 );
                   1438:        }
                   1439:        c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
                   1440:        if ( ( hd = c2[d2] ) != 1 ) {
                   1441:                inv = invm(hd,mod);
                   1442:                for ( pc1 = c2 + d2; pc1 >= c2; pc1-- )
                   1443:                        *pc1 = (*pc1 * inv) % mod;
                   1444:        } else
                   1445:                inv = 1;
                   1446:        for ( i = dd, ct = c1+d1; i >= 0; i-- )
                   1447:                if ( tmp = *ct-- ) {
                   1448:                        tmp = mod - tmp;
                   1449:                        for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- )
                   1450:                                *pct = (*pc1 * tmp + *pct) % mod;
                   1451:                }
                   1452:        if ( inv != 1 ) {
                   1453:                for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )
                   1454:                        *pc1 = (*pc1 * inv) % mod;
                   1455:                for ( pc1 = c2, pct = c2+d2, inv = hd; pc1 <= pct; pc1++ )
                   1456:                        *pc1 = (*pc1 * inv) % mod;
                   1457:        }
                   1458:        for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
                   1459:        for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
                   1460:                *pct-- = *pc1--;
                   1461:        return( i );
                   1462: }
                   1463:
1.6     ! noro     1464: void Gcdum(int mod,UM p1,UM p2,UM pr)
1.1       noro     1465: {
                   1466:        register int *sp,*dp;
                   1467:        register int i,inv;
                   1468:        UM t1,t2,q,tum;
                   1469:        int drem;
                   1470:
                   1471:        if ( DEG(p1) < 0 )
                   1472:                cpyum(p2,pr);
                   1473:        else if ( DEG(p2) < 0 )
                   1474:                cpyum(p1,pr);
                   1475:        else {
                   1476:                if ( DEG(p1) >= DEG(p2) ) {
                   1477:                        t1 = p1; t2 = p2;
                   1478:                } else {
                   1479:                        t1 = p2; t2 = p1;
                   1480:                }
                   1481:                q = W_UMALLOC(DEG(t1));
                   1482:                while ( ( drem = Divum(mod,t1,t2,q) ) >= 0 ) {
                   1483:                        tum = t1; t1 = t2; t2 = tum; DEG(t2) = drem;
                   1484:                }
                   1485:                inv = invm(COEF(t2)[DEG(t2)],mod);
                   1486:                Mulsum(mod,t2,inv,pr);
                   1487:        }
                   1488: }
                   1489: #endif
                   1490:
1.6     ! noro     1491: void mult_mod_tab(UM p,int mod,UM *tab,UM r,int d)
1.1       noro     1492: {
                   1493:        UM w,w1,c;
                   1494:        int n,i;
                   1495:        int *pc;
                   1496:
                   1497:        w = W_UMALLOC(d); w1 = W_UMALLOC(d);
                   1498:        c = W_UMALLOC(1); DEG(c) = 0;
                   1499:        n = DEG(p); DEG(r) = -1;
                   1500:        for ( i = 0, pc = COEF(p); i <= n; i++ )
                   1501:                if ( pc[i] ) {
                   1502:                        COEF(c)[0] = pc[i];
                   1503:                        mulum(mod,tab[i],c,w);
                   1504:                        addum(mod,r,w,w1);
                   1505:                        cpyum(w1,r);
                   1506:                }
                   1507: }
                   1508:
1.6     ! noro     1509: void make_qmat(UM p,int mod,UM *tab,int ***mp)
1.1       noro     1510: {
                   1511:        int n,i,j;
                   1512:        int *c;
                   1513:        UM q,r;
                   1514:        int **mat;
                   1515:
                   1516:        n = DEG(p);
                   1517:        *mp = mat = almat(n,n);
                   1518:        for ( j = 0; j < n; j++ ) {
                   1519:                r = W_UMALLOC(DEG(tab[j])); q = W_UMALLOC(DEG(tab[j]));
                   1520:                cpyum(tab[j],r); DEG(r) = divum(mod,r,p,q);
                   1521:                for ( i = 0, c = COEF(r); i <= DEG(r); i++ )
                   1522:                        mat[i][j] = c[i];
                   1523:        }
                   1524:        for ( i = 0; i < n; i++ )
                   1525:                mat[i][i] = (mat[i][i]+mod-1) % mod;
                   1526: }
                   1527:
1.6     ! noro     1528: void null_mod(int **mat,int mod,int n,int *ind)
1.1       noro     1529: {
                   1530:        int i,j,l,s,h,inv;
                   1531:        int *t,*u;
                   1532:
                   1533:        bzero((char *)ind,n*sizeof(int));
                   1534:        ind[0] = 0;
                   1535:        for ( i = j = 0; j < n; i++, j++ ) {
                   1536:                for ( ; j < n; j++ ) {
                   1537:                        for ( l = i; l < n; l++ )
                   1538:                                if ( mat[l][j] )
                   1539:                                        break;
                   1540:                        if ( l < n ) {
                   1541:                                t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
                   1542:                        } else
                   1543:                                ind[j] = 1;
                   1544:                }
                   1545:                if ( j == n )
                   1546:                        break;
                   1547:                inv = invm(mat[i][j],mod);
                   1548:                for ( s = j, t = mat[i]; s < n; s++ )
                   1549:                        t[s] =  dmar(t[s],inv,0,mod);
                   1550:                for ( l = 0; l < n; l++ ) {
                   1551:                        if ( l == i )
                   1552:                                continue;
                   1553:                        for ( s = j, u = mat[l], h = (mod-u[j])%mod; s < n; s++ )
                   1554:                                u[s] = dmar(h,t[s],u[s],mod);
                   1555:                }
                   1556:        }
                   1557: }
                   1558:
1.6     ! noro     1559: void null_to_sol(int **mat,int *ind,int mod,int n,UM *r)
1.1       noro     1560: {
                   1561:        int i,j,k,l;
                   1562:        int *c;
                   1563:        UM w;
                   1564:
                   1565:        for ( i = 0, l = 0; i < n; i++ ) {
                   1566:                if ( !ind[i] )
                   1567:                        continue;
                   1568:                w = UMALLOC(n);
                   1569:                for ( j = k = 0, c = COEF(w); j < n; j++ )
                   1570:                        if ( ind[j] )
                   1571:                                c[j] = 0;
                   1572:                        else
                   1573:                                c[j] = mat[k++][i];
                   1574:                c[i] = mod-1;
                   1575:                for ( j = n; j >= 0; j-- )
                   1576:                        if ( c[j] )
                   1577:                                break;
                   1578:                DEG(w) = j;
                   1579:                r[l++] = w;
                   1580:        }
                   1581: }
                   1582: /*
                   1583: make_qmat(p,mod,tab,mp)
                   1584: null_mod(mat,mod,n,ind)
                   1585: null_to_sol(mat,ind,mod,n,r)
                   1586: */
                   1587:
1.6     ! noro     1588: void newddd(int mod,UM f,UM *r)
1.1       noro     1589: {
                   1590:        register int i,j;
                   1591:        int d,n;
                   1592:        UM q,s,t,u,v,w,g,x,m;
                   1593:        UM *base;
                   1594:
                   1595:        n = DEG(f); base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
                   1596:        w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
                   1597:        base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
                   1598:        t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
                   1599:        pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
                   1600:        for ( i = 2; i < n; i++ ) {
                   1601: /*             fprintf(stderr,"i=%d\n",i); */
                   1602:                mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
                   1603:                base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
                   1604:        }
                   1605:        v = W_UMALLOC(n); cpyum(f,v);
                   1606:        DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
                   1607:        x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
                   1608:        t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
                   1609:        for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
                   1610: /*             fprintf(stderr,"d=%d\n",d); */
                   1611:                for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
                   1612:                        if ( COEF(w)[i] ) {
                   1613:                                Mulsum(mod,base[i],COEF(w)[i],s);
                   1614:                                addum(mod,s,t,u); cpyum(u,t);
                   1615:                        }
                   1616:                cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
                   1617:                if ( DEG(g) >= 1 ) {
                   1618:                        berlekamp(g,mod,d,base,r+j); j += DEG(g)/d;
                   1619:                        Divum(mod,v,g,q); cpyum(q,v);
                   1620:                        DEG(w) = Divum(mod,w,v,q);
                   1621:                        for ( i = 0; i < DEG(v); i++ )
                   1622:                                DEG(base[i]) = Divum(mod,base[i],v,q);
                   1623:                }
                   1624:        }
                   1625:        if ( DEG(v) ) {
                   1626:                r[j] = UMALLOC(DEG(v)); cpyum(v,r[j]); j++;
                   1627:        }
                   1628:        r[j] = 0;
                   1629: }
                   1630:
1.6     ! noro     1631: int nfctr_mod(UM f,int mod)
1.1       noro     1632: {
                   1633:        register int i,j;
                   1634:        int d,n;
                   1635:        UM q,s,t,u,v,w,g,x,m;
                   1636:        UM *base;
                   1637:
                   1638:        n = DEG(f); base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
                   1639:        w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
                   1640:        base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
                   1641:        t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
                   1642:        pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
                   1643:        for ( i = 2; i < n; i++ ) {
                   1644: /*             fprintf(stderr,"i=%d\n",i); */
                   1645:                mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
                   1646:                base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
                   1647:        }
                   1648:        v = W_UMALLOC(n); cpyum(f,v);
                   1649:        DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
                   1650:        x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
                   1651:        t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
                   1652:        for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
                   1653: /*             fprintf(stderr,"d=%d\n",d); */
                   1654:                for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
                   1655:                        if ( COEF(w)[i] ) {
                   1656:                                Mulsum(mod,base[i],COEF(w)[i],s);
                   1657:                                addum(mod,s,t,u); cpyum(u,t);
                   1658:                        }
                   1659:                cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
                   1660:                if ( DEG(g) >= 1 ) {
                   1661:                        j += DEG(g)/d;
                   1662:                        Divum(mod,v,g,q); cpyum(q,v);
                   1663:                        DEG(w) = Divum(mod,w,v,q);
                   1664:                        for ( i = 0; i < DEG(v); i++ )
                   1665:                                DEG(base[i]) = Divum(mod,base[i],v,q);
                   1666:                }
                   1667:        }
                   1668:        if ( DEG(v) ) j++;
                   1669:        return j;
                   1670: }
                   1671:
1.6     ! noro     1672: int irred_check(UM f,int mod)
1.1       noro     1673: {
                   1674:        register int i,j;
                   1675:        int d,n;
                   1676:        UM q,s,t,u,v,w,g,x,m,f1,b;
                   1677:        UM *base;
                   1678:
                   1679:        if ( (n = DEG(f)) == 1 )
                   1680:                return 1;
                   1681:        t = W_UMALLOC(n); s = W_UMALLOC(n); g = W_UMALLOC(n);
                   1682:        f1 = W_UMALLOC(n); b = W_UMALLOC(n);
                   1683:        diffum(mod,f,t); cpyum(f,s); Gcdum(mod,t,s,g);
                   1684:        if ( DEG(g) )
                   1685:                return 0;
                   1686:
                   1687:        base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
                   1688:        w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
                   1689:        base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
                   1690:        t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
                   1691:        pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
                   1692:        for ( i = 2; i < n; i++ ) {
                   1693: /*             fprintf(stderr,"i=%d\n",i); */
                   1694:                mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
                   1695:                base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
                   1696:        }
                   1697:        v = W_UMALLOC(n); cpyum(f,v);
                   1698:        DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
                   1699:        x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
                   1700:        t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
                   1701:        for ( j = 0, d = 1; 2*d <= n; d++ ) {
                   1702: /*             fprintf(stderr,"d=%d\n",d); */
                   1703:                for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
                   1704:                        if ( COEF(w)[i] ) {
                   1705:                                Mulsum(mod,base[i],COEF(w)[i],s);
                   1706:                                addum(mod,s,t,u); cpyum(u,t);
                   1707:                        }
                   1708:                cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
                   1709:                if ( DEG(g) >= 1 )
                   1710:                        return 0;
                   1711:        }
                   1712:        return 1;
                   1713: }
                   1714:
1.6     ! noro     1715: int berlekamp(UM p,int mod,int df,UM *tab,UM *r)
1.1       noro     1716: {
                   1717:        int n,i,j,k,nf,d,nr;
                   1718:        int **mat;
                   1719:        int *ind;
                   1720:        UM mp,w,q,gcd,w1,w2;
                   1721:        UM *u;
                   1722:        int *root;
                   1723:
                   1724:        n = DEG(p);
                   1725:        ind = ALLOCA(n*sizeof(int));
                   1726:        make_qmat(p,mod,tab,&mat);
                   1727:        null_mod(mat,mod,n,ind);
                   1728:        for ( i = 0, d = 0; i < n; i++ )
                   1729:                if ( ind[i] )
                   1730:                        d++;
                   1731:        if ( d == 1 ) {
                   1732:                r[0] = UMALLOC(n); cpyum(p,r[0]); return 1;
                   1733:        }
                   1734:        u = ALLOCA(d*sizeof(UM *));
                   1735:        r[0] = UMALLOC(n); cpyum(p,r[0]);
                   1736:        null_to_sol(mat,ind,mod,n,u);
                   1737:        root = ALLOCA(d*sizeof(int));
                   1738:        w = W_UMALLOC(n); mp = W_UMALLOC(d);
                   1739:        w1 = W_UMALLOC(n); w2 = W_UMALLOC(n);
                   1740:        for ( i = 1, nf = 1; i < d; i++ ) {
                   1741:                minipoly_mod(mod,u[i],p,mp);
                   1742:                nr = find_root(mod,mp,root);
                   1743:                for ( j = 0; j < nf; j++ ) {
                   1744:                        if ( DEG(r[j]) == df )
                   1745:                                continue;
                   1746:                        for ( k = 0; k < nr; k++ ) {
                   1747:                                cpyum(u[i],w1); cpyum(r[j],w2);
                   1748:                                COEF(w1)[0] = (mod-root[k]) % mod;
                   1749:                                gcdum(mod,w1,w2,w);
                   1750:                                if ( DEG(w) > 0 && DEG(w) < DEG(r[j]) ) {
                   1751:                                        gcd = UMALLOC(DEG(w));
                   1752:                                        q = UMALLOC(DEG(r[j])-DEG(w));
                   1753:                                        cpyum(w,gcd); divum(mod,r[j],w,q);
                   1754:                                        r[j] = q; r[nf++] = gcd;
                   1755:                                }
                   1756:                                if ( nf == d )
                   1757:                                        return d;
                   1758:                        }
                   1759:                }
                   1760:        }
1.6     ! noro     1761:        /* NOTREACHED */
        !          1762:        error("berlekamp : cannot happen");
        !          1763:        return -1;
1.1       noro     1764: }
                   1765:
1.6     ! noro     1766: void minipoly_mod(int mod,UM f,UM p,UM mp)
1.1       noro     1767: {
                   1768:        struct p_pair *list,*l,*l1,*lprev;
                   1769:        int n,d;
                   1770:        UM u,p0,p1,np0,np1,q,w;
                   1771:
                   1772:        list = (struct p_pair *)MALLOC(sizeof(struct p_pair));
                   1773:        list->p0 = u = W_UMALLOC(0); DEG(u) = 0; COEF(u)[0] = 1;
                   1774:        list->p1 = W_UMALLOC(0); cpyum(list->p0,list->p1);
                   1775:        list->next = 0;
                   1776:        n = DEG(p); w = UMALLOC(2*n);
                   1777:        p0 = UMALLOC(2*n); cpyum(list->p0,p0);
                   1778:        p1 = UMALLOC(2*n); cpyum(list->p1,p1);
                   1779:        q = W_UMALLOC(2*n);
                   1780:        while ( 1 ) {
                   1781:                COEF(p0)[DEG(p0)] = 0; DEG(p0)++; COEF(p0)[DEG(p0)] = 1;
                   1782:                mulum(mod,f,p1,w); DEG(w) = divum(mod,w,p,q); cpyum(w,p1);
                   1783:                np0 = UMALLOC(n); np1 = UMALLOC(n);
                   1784:                lnf_mod(mod,n,p0,p1,list,np0,np1);
                   1785:                if ( DEG(np1) < 0 ) {
                   1786:                        cpyum(np0,mp); return;
                   1787:                } else {
                   1788:                        l1 = (struct p_pair *)MALLOC(sizeof(struct p_pair));
                   1789:                        l1->p0 = np0; l1->p1 = np1;
                   1790:                        for ( l = list, lprev = 0, d = DEG(np1);
                   1791:                                l && (DEG(l->p1) > d); lprev = l, l = l->next );
                   1792:                        if ( lprev ) {
                   1793:                                lprev->next = l1; l1->next = l;
                   1794:                        } else {
                   1795:                                l1->next = list; list = l1;
                   1796:                        }
                   1797:                }
                   1798:        }
                   1799: }
                   1800:
1.6     ! noro     1801: void lnf_mod(int mod,int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1)
1.1       noro     1802: {
                   1803:        int inv,h,d1;
                   1804:        UM t0,t1,s0,s1;
                   1805:        struct p_pair *l;
                   1806:
                   1807:        cpyum(p0,np0); cpyum(p1,np1);
                   1808:        t0 = W_UMALLOC(n); t1 = W_UMALLOC(n);
                   1809:        s0 = W_UMALLOC(n); s1 = W_UMALLOC(n);
                   1810:        for ( l = list; l; l = l->next ) {
                   1811:                d1 = DEG(np1);
                   1812:                if ( d1 == DEG(l->p1) ) {
                   1813:                        inv = invm((mod-COEF(l->p1)[d1])%mod,mod);
                   1814:                        h = dmar(COEF(np1)[d1],inv,0,mod);
                   1815:                        mulsum(mod,l->p0,h,t0); addum(mod,np0,t0,s0); cpyum(s0,np0);
                   1816:                        mulsum(mod,l->p1,h,t1); addum(mod,np1,t1,s1); cpyum(s1,np1);
                   1817:                }
                   1818:        }
                   1819: }
                   1820:
1.6     ! noro     1821: int find_root(int mod,UM p,int *root)
1.1       noro     1822: {
                   1823:        UM *r;
                   1824:        int i,j;
                   1825:
                   1826:        r = ALLOCA((DEG(p)+1)*sizeof(UM));
                   1827:        ddd(mod,p,r);
                   1828:        for ( i = 0, j = 0; r[i]; i++ )
                   1829:                if ( DEG(r[i]) == 1 )
                   1830:                        root[j++] = (mod - COEF(r[i])[0]) % mod;
                   1831:        return j;
                   1832: }
                   1833:
1.6     ! noro     1834: void showum(UM p)
1.1       noro     1835: {
                   1836:        int i;
                   1837:        int *c;
                   1838:
                   1839:        for ( i = DEG(p), c = COEF(p); i >= 0; i-- )
                   1840:                if ( c[i] )
                   1841:                        printf("+%dx^%d",c[i],i);
                   1842:        printf("\n");
                   1843: }
                   1844:
1.6     ! noro     1845: void showumat(int **mat,int n)
1.1       noro     1846: {
                   1847:        int i,j;
                   1848:
                   1849:        for ( i = 0; i < n; i++ ) {
                   1850:                for ( j = 0; j < n; j++ )
                   1851:                        printf("%d ",mat[i][j]);
                   1852:                printf("\n");
                   1853:        }
                   1854: }

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