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

Annotation of OpenXM_contrib2/asir2000/engine/nd.c, Revision 1.24

1.24    ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.23 2003/08/05 05:01:18 noro Exp $ */
1.2       noro        2:
1.1       noro        3: #include "ca.h"
                      4: #include "inline.h"
                      5:
                      6: #if defined(__GNUC__)
                      7: #define INLINE inline
                      8: #elif defined(VISUAL)
                      9: #define INLINE __inline
                     10: #else
                     11: #define INLINE
                     12: #endif
                     13:
                     14: #define REDTAB_LEN 32003
                     15:
                     16: typedef struct oPGeoBucket {
                     17:        int m;
                     18:        struct oND *body[32];
                     19: } *PGeoBucket;
                     20:
                     21: typedef struct oND {
                     22:        struct oNM *body;
                     23:        int nv;
                     24:        int sugar;
                     25: } *ND;
                     26:
1.3       noro       27: typedef struct oNDV {
                     28:        struct oNMV *body;
                     29:        int nv;
                     30:        int sugar;
                     31:        int len;
                     32: } *NDV;
                     33:
1.1       noro       34: typedef struct oNM {
                     35:        struct oNM *next;
1.14      noro       36:        union {
                     37:                int m;
                     38:                Q z;
                     39:        } c;
1.1       noro       40:        int td;
                     41:        unsigned int dl[1];
                     42: } *NM;
                     43:
1.3       noro       44: typedef struct oNMV {
1.14      noro       45:        union {
                     46:                int m;
                     47:                Q z;
                     48:        } c;
1.3       noro       49:        int td;
                     50:        unsigned int dl[1];
                     51: } *NMV;
                     52:
1.13      noro       53: typedef struct oRHist {
                     54:        struct oRHist *next;
                     55:        int index;
1.20      noro       56:        int td,sugar;
1.13      noro       57:        unsigned int dl[1];
                     58: } *RHist;
                     59:
1.1       noro       60: typedef struct oND_pairs {
                     61:        struct oND_pairs *next;
                     62:        int i1,i2;
                     63:        int td,sugar;
                     64:        unsigned int lcm[1];
                     65: } *ND_pairs;
                     66:
1.21      noro       67: double nd_scale=2;
1.1       noro       68: static unsigned int **nd_bound;
1.19      noro       69: int nd_nvar;
1.1       noro       70: int is_rlex;
                     71: int nd_epw,nd_bpe,nd_wpd;
                     72: unsigned int nd_mask[32];
                     73: unsigned int nd_mask0,nd_mask1;
1.20      noro       74:
1.1       noro       75: NM _nm_free_list;
                     76: ND _nd_free_list;
                     77: ND_pairs _ndp_free_list;
1.20      noro       78:
                     79: static NDV *nd_ps;
                     80: static NDV *nd_psq;
                     81: int *nd_psl;
                     82: RHist *nd_psh;
                     83: int nd_psn,nd_pslen;
                     84:
1.13      noro       85: RHist *nd_red;
1.1       noro       86: int nd_red_len;
                     87:
                     88: int nd_found,nd_create,nd_notfirst;
1.13      noro       89: int nm_adv;
1.20      noro       90: int nmv_adv;
                     91: int nmv_len;
                     92: NDV ndv_red;
1.1       noro       93:
1.20      noro       94: extern int Top,Reverse;
1.1       noro       95:
                     96: #define HTD(d) ((d)->body->td)
                     97: #define HDL(d) ((d)->body->dl)
1.14      noro       98: #define HCM(d) ((d)->body->c.m)
1.16      noro       99: #define HCQ(d) ((d)->body->c.z)
1.14      noro      100: #define CM(a) ((a)->c.m)
1.16      noro      101: #define CQ(a) ((a)->c.z)
1.14      noro      102: #define DL(a) ((a)->dl)
                    103: #define TD(a) ((a)->td)
                    104: #define SG(a) ((a)->sugar)
                    105: #define LEN(a) ((a)->len)
1.1       noro      106:
1.20      noro      107: #define NM_ADV(m) (m = (NM)(((char *)m)+nm_adv))
1.15      noro      108: #define NEWRHist(r) \
                    109: ((r)=(RHist)MALLOC(sizeof(struct oRHist)+(nd_wpd-1)*sizeof(unsigned int)))
1.1       noro      110: #define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list)
                    111: #define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list)
                    112: #define MKND(n,m,d) if(!_nd_free_list)_ND_alloc(); (d)=_nd_free_list; _nd_free_list = (ND)BDY(_nd_free_list); (d)->nv=(n); BDY(d)=(m)
                    113:
1.13      noro      114: #define NEXTRHist(r,c) \
                    115: if(!(r)){NEWRHist(r);(c)=(r);}else{NEWRHist(NEXT(c));(c)=NEXT(c);}
1.1       noro      116: #define NEXTNM(r,c) \
                    117: if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);}
                    118: #define NEXTNM2(r,c,s) \
                    119: if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
                    120: #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)
                    121: #define FREENDP(m) NEXT(m)=_ndp_free_list; _ndp_free_list=(m)
                    122: #define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m)
                    123:
                    124: #define NEXTND_pairs(r,c) \
                    125: if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
                    126:
1.20      noro      127: void nd_removecont(int mod,ND p);
1.21      noro      128: void nd_removecont2(ND p1,ND p2);
1.20      noro      129: void ndv_removecont(int mod,NDV p);
1.23      noro      130: void ndv_dehomogenize(NDV p);
1.16      noro      131: void ndv_mul_c_q(NDV p,Q mul);
                    132: void nd_mul_c_q(ND p,Q mul);
                    133:
1.20      noro      134: void GC_gcollect();
                    135: NODE append_one(NODE,int);
                    136:
1.21      noro      137: void removecont_array(Q *c,int n);
1.1       noro      138: ND_pairs crit_B( ND_pairs d, int s );
                    139: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp);
1.23      noro      140: void nd_gr_trace(LIST f,LIST v,int m,int homo,struct order_spec *ord,LIST *rp);
                    141: void nd_setup(int mod,NODE f);
                    142: void nd_setup_trace(int mod,NODE f);
1.16      noro      143: int nd_newps(int mod,ND a);
1.20      noro      144: int nd_newps_trace(int mod,ND nf,ND nfq);
1.1       noro      145: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest );
                    146: NODE update_base(NODE nd,int ndp);
                    147: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest );
                    148: int crit_2( int dp1, int dp2 );
                    149: ND_pairs crit_F( ND_pairs d1 );
                    150: ND_pairs crit_M( ND_pairs d1 );
                    151: ND_pairs nd_newpairs( NODE g, int t );
                    152: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t);
1.23      noro      153: NODE nd_gb(int m);
                    154: NODE nd_gb_trace(int m);
1.1       noro      155: void nd_free_private_storage();
                    156: void _NM_alloc();
                    157: void _ND_alloc();
                    158: int ndl_td(unsigned int *d);
1.23      noro      159: int ndl_dehomogenize(unsigned int *p);
1.19      noro      160: ND nd_add(int mod,ND p1,ND p2);
1.17      noro      161: ND nd_add_q(ND p1,ND p2);
                    162: ND nd_mul_nm(int mod,ND p,NM m0);
                    163: ND nd_mul_ind_nm(int mod,int index,NM m0);
1.16      noro      164: int nd_sp(int mod,ND_pairs p,ND *nf);
1.6       noro      165: int nd_find_reducer(ND g);
1.23      noro      166: int nd_find_reducer_direct(ND g,NDV *ps,int len);
1.16      noro      167: int nd_nf(int mod,ND g,int full,ND *nf);
1.1       noro      168: ND nd_reduce(ND p1,ND p2);
                    169: ND nd_reduce_special(ND p1,ND p2);
1.23      noro      170: NODE nd_reduceall(int m,NODE f);
1.1       noro      171: void nd_free(ND p);
1.23      noro      172: void ndv_free(NDV p);
1.1       noro      173: void ndl_print(unsigned int *dl);
                    174: void nd_print(ND p);
1.16      noro      175: void nd_print_q(ND p);
                    176: void ndv_print(NDV p);
                    177: void ndv_print_q(NDV p);
1.1       noro      178: void ndp_print(ND_pairs d);
                    179: int nd_length(ND p);
1.19      noro      180: void nd_mul_c(int mod,ND p,int mul);
1.1       noro      181: void nd_free_redlist();
                    182: void nd_append_red(unsigned int *d,int td,int i);
                    183: unsigned int *nd_compute_bound(ND p);
1.5       noro      184: unsigned int *dp_compute_bound(DP p);
1.20      noro      185: ND_pairs nd_reconstruct(int mod,int trace,ND_pairs ndp);
1.23      noro      186: void nd_reconstruct_direct(int mod,NDV *ps,int len);
1.1       noro      187: void nd_setup_parameters();
1.11      noro      188: void nd_realloc(ND p,int obpe);
1.6       noro      189: ND nd_copy(ND p);
1.1       noro      190: void ndl_dup(int obpe,unsigned int *d,unsigned int *r);
1.4       noro      191:
                    192: #define NMV_ADV(m) (m = (NMV)(((char *)m)+nmv_adv))
                    193: #define NEWNDV(d) ((d)=(NDV)MALLOC(sizeof(struct oNDV)))
1.14      noro      194: #define MKNDV(n,m,l,d) NEWNDV(d); NV(d)=(n); BDY(d)=(m); LEN(d) = l;
1.19      noro      195: void ndv_mul_c(int mod,NDV p,int mul);
                    196: ND ndv_add(int mod,ND p1,NDV p2);
1.16      noro      197: ND ndv_add_q(ND p1,NDV p2);
                    198: NDV ndtondv(int mod,ND p);
1.23      noro      199: ND ndvtond(int mod,NDV p);
1.16      noro      200: void ndv_mul_nm(int mod,NDV pv,NM m,NDV r);
                    201: ND ndv_mul_nm_create(int mod,NDV p,NM m0);
1.11      noro      202: void ndv_realloc(NDV p,int obpe,int oadv);
1.16      noro      203: NDV dptondv(int,DP);
                    204: DP ndvtodp(int,NDV);
1.17      noro      205: ND dptond(int,DP);
                    206: DP ndtodp(int,ND);
1.1       noro      207:
                    208: void nd_free_private_storage()
                    209: {
                    210:        _nd_free_list = 0;
                    211:        _nm_free_list = 0;
1.5       noro      212:        _ndp_free_list = 0;
1.13      noro      213:        bzero(nd_red,sizeof(REDTAB_LEN*sizeof(RHist)));
1.1       noro      214:        GC_gcollect();
                    215: }
                    216:
                    217: void _NM_alloc()
                    218: {
                    219:        NM p;
                    220:        int i;
                    221:
1.11      noro      222:        for ( i = 0; i < 1024; i++ ) {
1.1       noro      223:                p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
                    224:                p->next = _nm_free_list; _nm_free_list = p;
                    225:        }
                    226: }
                    227:
                    228: void _ND_alloc()
                    229: {
                    230:        ND p;
                    231:        int i;
                    232:
                    233:        for ( i = 0; i < 1024; i++ ) {
                    234:                p = (ND)GC_malloc(sizeof(struct oND));
                    235:                p->body = (NM)_nd_free_list; _nd_free_list = p;
                    236:        }
                    237: }
                    238:
                    239: void _NDP_alloc()
                    240: {
                    241:        ND_pairs p;
                    242:        int i;
                    243:
1.11      noro      244:        for ( i = 0; i < 1024; i++ ) {
1.1       noro      245:                p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs)
                    246:                        +(nd_wpd-1)*sizeof(unsigned int));
                    247:                p->next = _ndp_free_list; _ndp_free_list = p;
                    248:        }
                    249: }
                    250:
                    251: INLINE nd_length(ND p)
                    252: {
                    253:        NM m;
                    254:        int i;
                    255:
                    256:        if ( !p )
                    257:                return 0;
                    258:        else {
                    259:                for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
                    260:                return i;
                    261:        }
                    262: }
                    263:
                    264: INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2)
                    265: {
                    266:        unsigned int u1,u2;
                    267:        int i,j;
                    268:
                    269:        switch ( nd_bpe ) {
                    270:                case 4:
                    271:                        for ( i = 0; i < nd_wpd; i++ ) {
                    272:                                u1 = d1[i]; u2 = d2[i];
                    273:                                if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0;
                    274:                                if ( (u1&0xf000000) < (u2&0xf000000) ) return 0;
                    275:                                if ( (u1&0xf00000) < (u2&0xf00000) ) return 0;
                    276:                                if ( (u1&0xf0000) < (u2&0xf0000) ) return 0;
                    277:                                if ( (u1&0xf000) < (u2&0xf000) ) return 0;
                    278:                                if ( (u1&0xf00) < (u2&0xf00) ) return 0;
                    279:                                if ( (u1&0xf0) < (u2&0xf0) ) return 0;
                    280:                                if ( (u1&0xf) < (u2&0xf) ) return 0;
                    281:                        }
                    282:                        return 1;
                    283:                        break;
                    284:                case 6:
                    285:                        for ( i = 0; i < nd_wpd; i++ ) {
                    286:                                u1 = d1[i]; u2 = d2[i];
                    287:                                if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0;
                    288:                                if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0;
                    289:                                if ( (u1&0x3f000) < (u2&0x3f000) ) return 0;
                    290:                                if ( (u1&0xfc0) < (u2&0xfc0) ) return 0;
                    291:                                if ( (u1&0x3f) < (u2&0x3f) ) return 0;
                    292:                        }
                    293:                        return 1;
                    294:                        break;
                    295:                case 8:
                    296:                        for ( i = 0; i < nd_wpd; i++ ) {
                    297:                                u1 = d1[i]; u2 = d2[i];
                    298:                                if ( (u1&0xff000000) < (u2&0xff000000) ) return 0;
                    299:                                if ( (u1&0xff0000) < (u2&0xff0000) ) return 0;
                    300:                                if ( (u1&0xff00) < (u2&0xff00) ) return 0;
                    301:                                if ( (u1&0xff) < (u2&0xff) ) return 0;
                    302:                        }
                    303:                        return 1;
                    304:                        break;
                    305:                case 16:
                    306:                        for ( i = 0; i < nd_wpd; i++ ) {
                    307:                                u1 = d1[i]; u2 = d2[i];
                    308:                                if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0;
                    309:                                if ( (u1&0xffff) < (u2&0xffff) ) return 0;
                    310:                        }
                    311:                        return 1;
                    312:                        break;
                    313:                case 32:
                    314:                        for ( i = 0; i < nd_wpd; i++ )
                    315:                                if ( d1[i] < d2[i] ) return 0;
                    316:                        return 1;
                    317:                        break;
                    318:                default:
                    319:                        for ( i = 0; i < nd_wpd; i++ ) {
                    320:                                u1 = d1[i]; u2 = d2[i];
                    321:                                for ( j = 0; j < nd_epw; j++ )
                    322:                                        if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0;
                    323:                        }
                    324:                        return 1;
                    325:        }
                    326: }
                    327:
1.23      noro      328: /* returns the exponent of homo variable */
                    329:
                    330: int ndl_dehomogenize(unsigned int *d)
                    331: {
                    332:        unsigned int mask;
                    333:        unsigned int h;
                    334:        int i;
                    335:
                    336:        if ( is_rlex ) {
                    337:                if ( nd_bpe == 32 ) {
                    338:                        h = d[0];
                    339:                        for ( i = 1; i < nd_wpd; i++ )
                    340:                                d[i-1] = d[i];
                    341:                        d[i-1] = 0;
                    342:                        return h;
                    343:                } else {
                    344:                        mask = (1<<(nd_epw*nd_bpe))-1;
                    345:                        h = (d[0]>>((nd_epw-1)*nd_bpe))&nd_mask0;
                    346:                        for ( i = 0; i < nd_wpd; i++ )
                    347:                                d[i] = ((d[i]<<nd_bpe)&mask)
                    348:                                        |(i+1<nd_wpd?((d[i+1]>>((nd_epw-1)*nd_bpe))&nd_mask0):0);
                    349:                        return h;
                    350:                }
                    351:        } else {
                    352:                /* do nothing */
                    353:                return d[nd_wpd-1];
                    354:        }
                    355: }
                    356:
1.1       noro      357: void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
                    358: {
                    359:        unsigned int t1,t2,u,u1,u2;
                    360:        int i,j;
                    361:
                    362:        switch ( nd_bpe ) {
                    363:                case 4:
                    364:                        for ( i = 0; i < nd_wpd; i++ ) {
                    365:                                u1 = d1[i]; u2 = d2[i];
                    366:                                t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2;
                    367:                                t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2;
                    368:                                t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2;
                    369:                                t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2;
                    370:                                t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2;
                    371:                                t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2;
                    372:                                t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2;
                    373:                                t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2;
                    374:                                d[i] = u;
                    375:                        }
                    376:                        break;
                    377:                case 6:
                    378:                        for ( i = 0; i < nd_wpd; i++ ) {
                    379:                                u1 = d1[i]; u2 = d2[i];
                    380:                                t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2;
                    381:                                t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2;
                    382:                                t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2;
                    383:                                t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2;
                    384:                                t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2;
                    385:                                d[i] = u;
                    386:                        }
                    387:                        break;
                    388:                case 8:
                    389:                        for ( i = 0; i < nd_wpd; i++ ) {
                    390:                                u1 = d1[i]; u2 = d2[i];
                    391:                                t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2;
                    392:                                t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2;
                    393:                                t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2;
                    394:                                t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2;
                    395:                                d[i] = u;
                    396:                        }
                    397:                        break;
                    398:                case 16:
                    399:                        for ( i = 0; i < nd_wpd; i++ ) {
                    400:                                u1 = d1[i]; u2 = d2[i];
                    401:                                t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2;
                    402:                                t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2;
                    403:                                d[i] = u;
                    404:                        }
                    405:                        break;
                    406:                case 32:
                    407:                        for ( i = 0; i < nd_wpd; i++ ) {
                    408:                                u1 = d1[i]; u2 = d2[i];
                    409:                                d[i] = u1>u2?u1:u2;
                    410:                        }
                    411:                        break;
                    412:                default:
                    413:                        for ( i = 0; i < nd_wpd; i++ ) {
                    414:                                u1 = d1[i]; u2 = d2[i];
                    415:                                for ( j = 0, u = 0; j < nd_epw; j++ ) {
                    416:                                        t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2;
                    417:                                }
                    418:                                d[i] = u;
                    419:                        }
                    420:                        break;
                    421:        }
                    422: }
                    423:
                    424: int ndl_td(unsigned int *d)
                    425: {
                    426:        unsigned int t,u;
                    427:        int i,j;
                    428:
                    429:        for ( t = 0, i = 0; i < nd_wpd; i++ ) {
                    430:                u = d[i];
                    431:                for ( j = 0; j < nd_epw; j++, u>>=nd_bpe )
                    432:                        t += (u&nd_mask0);
                    433:        }
                    434:        return t;
                    435: }
                    436:
                    437: INLINE int ndl_compare(unsigned int *d1,unsigned int *d2)
                    438: {
                    439:        int i;
                    440:
                    441:        for ( i = 0; i < nd_wpd; i++, d1++, d2++ )
                    442:                if ( *d1 > *d2 )
                    443:                        return is_rlex ? -1 : 1;
                    444:                else if ( *d1 < *d2 )
                    445:                        return is_rlex ? 1 : -1;
                    446:        return 0;
                    447: }
                    448:
                    449: INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)
                    450: {
                    451:        int i;
                    452:
                    453:        for ( i = 0; i < nd_wpd; i++ )
                    454:                if ( d1[i] != d2[i] )
                    455:                        return 0;
                    456:        return 1;
                    457: }
                    458:
1.6       noro      459: INLINE void ndl_copy(unsigned int *d1,unsigned int *d2)
                    460: {
                    461:        int i;
                    462:
                    463:        switch ( nd_wpd ) {
                    464:                case 1:
                    465:                        d2[0] = d1[0];
                    466:                        break;
                    467:                case 2:
                    468:                        d2[0] = d1[0];
                    469:                        d2[1] = d1[1];
                    470:                        break;
                    471:                default:
                    472:                        for ( i = 0; i < nd_wpd; i++ )
                    473:                                d2[i] = d1[i];
                    474:                        break;
                    475:        }
                    476: }
                    477:
1.1       noro      478: INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d)
                    479: {
                    480:        int i;
                    481:
1.6       noro      482:        switch ( nd_wpd ) {
                    483:                case 1:
                    484:                        d[0] = d1[0]+d2[0];
                    485:                        break;
                    486:                case 2:
                    487:                        d[0] = d1[0]+d2[0];
                    488:                        d[1] = d1[1]+d2[1];
                    489:                        break;
                    490:                default:
                    491:                        for ( i = 0; i < nd_wpd; i++ )
                    492:                                d[i] = d1[i]+d2[i];
                    493:                        break;
                    494:        }
                    495: }
                    496:
                    497: INLINE void ndl_add2(unsigned int *d1,unsigned int *d2)
                    498: {
                    499:        int i;
                    500:
                    501:        switch ( nd_wpd ) {
                    502:                case 1:
                    503:                        d2[0] += d1[0];
                    504:                        break;
                    505:                case 2:
                    506:                        d2[0] += d1[0];
                    507:                        d2[1] += d1[1];
                    508:                        break;
                    509:                default:
                    510:                        for ( i = 0; i < nd_wpd; i++ )
                    511:                                d2[i] += d1[i];
                    512:                        break;
1.1       noro      513:        }
                    514: }
                    515:
                    516: void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
                    517: {
                    518:        int i;
                    519:
                    520:        for ( i = 0; i < nd_wpd; i++ )
                    521:                d[i] = d1[i]-d2[i];
                    522: }
                    523:
                    524: int ndl_disjoint(unsigned int *d1,unsigned int *d2)
                    525: {
                    526:        unsigned int t1,t2,u,u1,u2;
                    527:        int i,j;
                    528:
                    529:        switch ( nd_bpe ) {
                    530:                case 4:
                    531:                        for ( i = 0; i < nd_wpd; i++ ) {
                    532:                                u1 = d1[i]; u2 = d2[i];
                    533:                                t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0;
                    534:                                t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0;
                    535:                                t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0;
                    536:                                t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0;
                    537:                                t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0;
                    538:                                t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0;
                    539:                                t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0;
                    540:                                t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0;
                    541:                        }
                    542:                        return 1;
                    543:                        break;
                    544:                case 6:
                    545:                        for ( i = 0; i < nd_wpd; i++ ) {
                    546:                                u1 = d1[i]; u2 = d2[i];
                    547:                                t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0;
                    548:                                t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0;
                    549:                                t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0;
                    550:                                t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0;
                    551:                                t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0;
                    552:                        }
                    553:                        return 1;
                    554:                        break;
                    555:                case 8:
                    556:                        for ( i = 0; i < nd_wpd; i++ ) {
                    557:                                u1 = d1[i]; u2 = d2[i];
                    558:                                t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0;
                    559:                                t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0;
                    560:                                t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0;
                    561:                                t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0;
                    562:                        }
                    563:                        return 1;
                    564:                        break;
                    565:                case 16:
                    566:                        for ( i = 0; i < nd_wpd; i++ ) {
                    567:                                u1 = d1[i]; u2 = d2[i];
                    568:                                t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0;
                    569:                                t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0;
                    570:                        }
                    571:                        return 1;
                    572:                        break;
                    573:                case 32:
                    574:                        for ( i = 0; i < nd_wpd; i++ )
                    575:                                if ( d1[i] && d2[i] ) return 0;
                    576:                        return 1;
                    577:                        break;
                    578:                default:
                    579:                        for ( i = 0; i < nd_wpd; i++ ) {
                    580:                                u1 = d1[i]; u2 = d2[i];
                    581:                                for ( j = 0; j < nd_epw; j++ ) {
                    582:                                        if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0;
                    583:                                        u1 >>= nd_bpe; u2 >>= nd_bpe;
                    584:                                }
                    585:                        }
                    586:                        return 1;
                    587:                        break;
                    588:        }
                    589: }
                    590:
1.5       noro      591: int ndl_check_bound2(int index,unsigned int *d2)
1.1       noro      592: {
1.5       noro      593:        unsigned int u2;
                    594:        unsigned int *d1;
                    595:        int i,j,ind,k;
1.1       noro      596:
1.5       noro      597:        d1 = nd_bound[index];
                    598:        ind = 0;
                    599:        switch ( nd_bpe ) {
                    600:                case 4:
                    601:                        for ( i = 0; i < nd_wpd; i++ ) {
                    602:                                u2 = d2[i];
                    603:                                if ( d1[ind++]+((u2>>28)&0xf) >= 0x10 ) return 1;
                    604:                                if ( d1[ind++]+((u2>>24)&0xf) >= 0x10 ) return 1;
                    605:                                if ( d1[ind++]+((u2>>20)&0xf) >= 0x10 ) return 1;
                    606:                                if ( d1[ind++]+((u2>>16)&0xf) >= 0x10 ) return 1;
                    607:                                if ( d1[ind++]+((u2>>12)&0xf) >= 0x10 ) return 1;
                    608:                                if ( d1[ind++]+((u2>>8)&0xf) >= 0x10 ) return 1;
                    609:                                if ( d1[ind++]+((u2>>4)&0xf) >= 0x10 ) return 1;
                    610:                                if ( d1[ind++]+(u2&0xf) >= 0x10 ) return 1;
                    611:                        }
                    612:                        return 0;
                    613:                        break;
                    614:                case 6:
                    615:                        for ( i = 0; i < nd_wpd; i++ ) {
                    616:                                u2 = d2[i];
                    617:                                if ( d1[ind++]+((u2>>24)&0x3f) >= 0x40 ) return 1;
                    618:                                if ( d1[ind++]+((u2>>18)&0x3f) >= 0x40 ) return 1;
                    619:                                if ( d1[ind++]+((u2>>12)&0x3f) >= 0x40 ) return 1;
                    620:                                if ( d1[ind++]+((u2>>6)&0x3f) >= 0x40 ) return 1;
                    621:                                if ( d1[ind++]+(u2&0x3f) >= 0x40 ) return 1;
                    622:                        }
                    623:                        return 0;
                    624:                        break;
                    625:                case 8:
                    626:                        for ( i = 0; i < nd_wpd; i++ ) {
                    627:                                u2 = d2[i];
                    628:                                if ( d1[ind++]+((u2>>24)&0xff) >= 0x100 ) return 1;
                    629:                                if ( d1[ind++]+((u2>>16)&0xff) >= 0x100 ) return 1;
                    630:                                if ( d1[ind++]+((u2>>8)&0xff) >= 0x100 ) return 1;
                    631:                                if ( d1[ind++]+(u2&0xff) >= 0x100 ) return 1;
                    632:                        }
                    633:                        return 0;
                    634:                        break;
                    635:                case 16:
                    636:                        for ( i = 0; i < nd_wpd; i++ ) {
                    637:                                u2 = d2[i];
                    638:                                if ( d1[ind++]+((u2>>16)&0xffff) > 0x10000 ) return 1;
                    639:                                if ( d1[ind++]+(u2&0xffff) > 0x10000 ) return 1;
                    640:                        }
                    641:                        return 0;
                    642:                        break;
                    643:                case 32:
                    644:                        for ( i = 0; i < nd_wpd; i++ )
                    645:                                if ( d1[i]+d2[i]<d1[i] ) return 1;
                    646:                        return 0;
                    647:                        break;
                    648:                default:
                    649:                        for ( i = 0; i < nd_wpd; i++ ) {
                    650:                                u2 = d2[i];
                    651:                                k = (nd_epw-1)*nd_bpe;
                    652:                                for ( j = 0; j < nd_epw; j++, k -= nd_bpe )
                    653:                                        if ( d1[ind++]+((u2>>k)&nd_mask0) > nd_mask0 ) return 1;
                    654:                        }
                    655:                        return 0;
                    656:                        break;
                    657:        }
1.1       noro      658: }
                    659:
1.23      noro      660: int ndl_check_bound2_direct(unsigned int *d1,unsigned int *d2)
                    661: {
                    662:        unsigned int u1,u2;
                    663:        int i,j,k;
                    664:
                    665:        switch ( nd_bpe ) {
                    666:                case 4:
                    667:                        for ( i = 0; i < nd_wpd; i++ ) {
                    668:                                u1 = d1[i]; u2 = d2[i];
                    669:                                if ( ((u1>>28)&0xf)+((u2>>28)&0xf) >= 0x10 ) return 1;
                    670:                                if ( ((u1>>24)&0xf)+((u2>>24)&0xf) >= 0x10 ) return 1;
                    671:                                if ( ((u1>>20)&0xf)+((u2>>20)&0xf) >= 0x10 ) return 1;
                    672:                                if ( ((u1>>16)&0xf)+((u2>>16)&0xf) >= 0x10 ) return 1;
                    673:                                if ( ((u1>>12)&0xf)+((u2>>12)&0xf) >= 0x10 ) return 1;
                    674:                                if ( ((u1>>8)&0xf)+((u2>>8)&0xf) >= 0x10 ) return 1;
                    675:                                if ( ((u1>>4)&0xf)+((u2>>4)&0xf) >= 0x10 ) return 1;
                    676:                                if ( (u1&0xf)+(u2&0xf) >= 0x10 ) return 1;
                    677:                        }
                    678:                        return 0;
                    679:                        break;
                    680:                case 6:
                    681:                        for ( i = 0; i < nd_wpd; i++ ) {
                    682:                                u1 = d1[i]; u2 = d2[i];
                    683:                                if ( ((u1>>24)&0x3f)+((u2>>24)&0x3f) >= 0x40 ) return 1;
                    684:                                if ( ((u1>>18)&0x3f)+((u2>>18)&0x3f) >= 0x40 ) return 1;
                    685:                                if ( ((u1>>12)&0x3f)+((u2>>12)&0x3f) >= 0x40 ) return 1;
                    686:                                if ( ((u1>>6)&0x3f)+((u2>>6)&0x3f) >= 0x40 ) return 1;
                    687:                                if ( (u1&0x3f)+(u2&0x3f) >= 0x40 ) return 1;
                    688:                        }
                    689:                        return 0;
                    690:                        break;
                    691:                case 8:
                    692:                        for ( i = 0; i < nd_wpd; i++ ) {
                    693:                                u1 = d1[i]; u2 = d2[i];
                    694:                                if ( ((u1>>24)&0xff)+((u2>>24)&0xff) >= 0x100 ) return 1;
                    695:                                if ( ((u1>>16)&0xff)+((u2>>16)&0xff) >= 0x100 ) return 1;
                    696:                                if ( ((u1>>8)&0xff)+((u2>>8)&0xff) >= 0x100 ) return 1;
                    697:                                if ( (u1&0xff)+(u2&0xff) >= 0x100 ) return 1;
                    698:                        }
                    699:                        return 0;
                    700:                        break;
                    701:                case 16:
                    702:                        for ( i = 0; i < nd_wpd; i++ ) {
                    703:                                u1 = d1[i]; u2 = d2[i];
                    704:                                        if ( ((u1>>16)&0xffff)+((u2>>16)&0xffff) > 0x10000 ) return 1;
                    705:                                        if ( (u2&0xffff)+(u2&0xffff) > 0x10000 ) return 1;
                    706:                        }
                    707:                        return 0;
                    708:                        break;
                    709:                case 32:
                    710:                        for ( i = 0; i < nd_wpd; i++ )
                    711:                                if ( d1[i]+d2[i]<d1[i] ) return 1;
                    712:                        return 0;
                    713:                        break;
                    714:                default:
                    715:                        for ( i = 0; i < nd_wpd; i++ ) {
                    716:                                u1 = d1[i]; u2 = d2[i];
                    717:                                k = (nd_epw-1)*nd_bpe;
                    718:                                for ( j = 0; j < nd_epw; j++, k -= nd_bpe )
                    719:                                        if ( ((u1>>k)&nd_mask0)+((u2>>k)&nd_mask0) > nd_mask0 )
                    720:                                                return 1;
                    721:                        }
                    722:                        return 0;
                    723:                        break;
                    724:        }
                    725: }
                    726:
1.6       noro      727: INLINE int ndl_hash_value(int td,unsigned int *d)
1.1       noro      728: {
                    729:        int i;
                    730:        int r;
                    731:
                    732:        r = td;
                    733:        for ( i = 0; i < nd_wpd; i++ )
                    734:                r = ((r<<16)+d[i])%REDTAB_LEN;
                    735:        return r;
                    736: }
                    737:
1.9       noro      738: INLINE int nd_find_reducer(ND g)
1.1       noro      739: {
1.13      noro      740:        RHist r;
1.6       noro      741:        int d,k,i;
1.1       noro      742:
                    743:        d = ndl_hash_value(HTD(g),HDL(g));
1.13      noro      744:        for ( r = nd_red[d], k = 0; r; r = NEXT(r), k++ ) {
1.14      noro      745:                if ( HTD(g) == TD(r) && ndl_equal(HDL(g),DL(r)) ) {
1.1       noro      746:                        if ( k > 0 ) nd_notfirst++;
                    747:                        nd_found++;
1.13      noro      748:                        return r->index;
1.1       noro      749:                }
                    750:        }
                    751:
1.13      noro      752:        if ( Reverse )
                    753:                for ( i = nd_psn-1; i >= 0; i-- ) {
                    754:                        r = nd_psh[i];
1.14      noro      755:                        if ( HTD(g) >= TD(r) && ndl_reducible(HDL(g),DL(r)) ) {
1.13      noro      756:                                nd_create++;
                    757:                                nd_append_red(HDL(g),HTD(g),i);
                    758:                                return i;
                    759:                        }
                    760:                }
                    761:        else
                    762:                for ( i = 0; i < nd_psn; i++ ) {
                    763:                        r = nd_psh[i];
1.14      noro      764:                        if ( HTD(g) >= TD(r) && ndl_reducible(HDL(g),DL(r)) ) {
1.13      noro      765:                                nd_create++;
                    766:                                nd_append_red(HDL(g),HTD(g),i);
                    767:                                return i;
                    768:                        }
1.1       noro      769:                }
1.6       noro      770:        return -1;
1.1       noro      771: }
                    772:
1.23      noro      773: INLINE int nd_find_reducer_direct(ND g,NDV *ps,int len)
                    774: {
                    775:        int i;
                    776:        NDV r;
                    777:
                    778:        if ( Reverse )
                    779:                for ( i = len-1; i >= 0; i-- ) {
                    780:                        r = ps[i];
                    781:                        if ( HTD(g) >= HTD(r) && ndl_reducible(HDL(g),HDL(r)) ) {
                    782:                                nd_append_red(HDL(g),HTD(g),i);
                    783:                                return i;
                    784:                        }
                    785:                }
                    786:        else
                    787:                for ( i = 0; i < len; i++ ) {
                    788:                        r = ps[i];
                    789:                        if ( HTD(g) >= HTD(r) && ndl_reducible(HDL(g),HDL(r)) ) {
                    790:                                nd_append_red(HDL(g),HTD(g),i);
                    791:                                return i;
                    792:                        }
                    793:                }
                    794:        return -1;
                    795: }
                    796:
1.19      noro      797: ND nd_add(int mod,ND p1,ND p2)
1.1       noro      798: {
                    799:        int n,c;
                    800:        int t;
                    801:        ND r;
                    802:        NM m1,m2,mr0,mr,s;
                    803:
                    804:        if ( !p1 )
                    805:                return p2;
                    806:        else if ( !p2 )
                    807:                return p1;
                    808:        else {
                    809:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
1.14      noro      810:                        if ( TD(m1) > TD(m2) )
1.1       noro      811:                                c = 1;
1.14      noro      812:                        else if ( TD(m1) < TD(m2) )
1.1       noro      813:                                c = -1;
                    814:                        else
1.14      noro      815:                                c = ndl_compare(DL(m1),DL(m2));
1.1       noro      816:                        switch ( c ) {
                    817:                                case 0:
1.19      noro      818:                                        t = ((CM(m1))+(CM(m2))) - mod;
1.1       noro      819:                                        if ( t < 0 )
1.19      noro      820:                                                t += mod;
1.1       noro      821:                                        s = m1; m1 = NEXT(m1);
                    822:                                        if ( t ) {
1.14      noro      823:                                                NEXTNM2(mr0,mr,s); CM(mr) = (t);
1.1       noro      824:                                        } else {
                    825:                                                FREENM(s);
                    826:                                        }
                    827:                                        s = m2; m2 = NEXT(m2); FREENM(s);
                    828:                                        break;
                    829:                                case 1:
                    830:                                        s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
                    831:                                        break;
                    832:                                case -1:
                    833:                                        s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
                    834:                                        break;
                    835:                        }
                    836:                }
                    837:                if ( !mr0 )
                    838:                        if ( m1 )
                    839:                                mr0 = m1;
                    840:                        else if ( m2 )
                    841:                                mr0 = m2;
                    842:                        else
                    843:                                return 0;
                    844:                else if ( m1 )
                    845:                        NEXT(mr) = m1;
                    846:                else if ( m2 )
                    847:                        NEXT(mr) = m2;
                    848:                else
                    849:                        NEXT(mr) = 0;
                    850:                BDY(p1) = mr0;
1.14      noro      851:                SG(p1) = MAX(SG(p1),SG(p2));
1.1       noro      852:                FREEND(p2);
                    853:                return p1;
                    854:        }
                    855: }
                    856:
1.17      noro      857: ND nd_add_q(ND p1,ND p2)
                    858: {
                    859:        int n,c;
                    860:        ND r;
                    861:        NM m1,m2,mr0,mr,s;
                    862:        Q t;
                    863:
                    864:        if ( !p1 )
                    865:                return p2;
                    866:        else if ( !p2 )
                    867:                return p1;
                    868:        else {
                    869:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                    870:                        if ( TD(m1) > TD(m2) )
                    871:                                c = 1;
                    872:                        else if ( TD(m1) < TD(m2) )
                    873:                                c = -1;
                    874:                        else
                    875:                                c = ndl_compare(DL(m1),DL(m2));
                    876:                        switch ( c ) {
                    877:                                case 0:
                    878:                                        addq(CQ(m1),CQ(m2),&t);
                    879:                                        s = m1; m1 = NEXT(m1);
                    880:                                        if ( t ) {
                    881:                                                NEXTNM2(mr0,mr,s); CQ(mr) = (t);
                    882:                                        } else {
                    883:                                                FREENM(s);
                    884:                                        }
                    885:                                        s = m2; m2 = NEXT(m2); FREENM(s);
                    886:                                        break;
                    887:                                case 1:
                    888:                                        s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
                    889:                                        break;
                    890:                                case -1:
                    891:                                        s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
                    892:                                        break;
                    893:                        }
                    894:                }
                    895:                if ( !mr0 )
                    896:                        if ( m1 )
                    897:                                mr0 = m1;
                    898:                        else if ( m2 )
                    899:                                mr0 = m2;
                    900:                        else
                    901:                                return 0;
                    902:                else if ( m1 )
                    903:                        NEXT(mr) = m1;
                    904:                else if ( m2 )
                    905:                        NEXT(mr) = m2;
                    906:                else
                    907:                        NEXT(mr) = 0;
                    908:                BDY(p1) = mr0;
                    909:                SG(p1) = MAX(SG(p1),SG(p2));
                    910:                FREEND(p2);
                    911:                return p1;
                    912:        }
                    913: }
                    914:
1.1       noro      915: #if 1
                    916: /* ret=1 : success, ret=0 : overflow */
1.16      noro      917: int nd_nf(int mod,ND g,int full,ND *rp)
1.1       noro      918: {
1.11      noro      919:        ND d;
1.1       noro      920:        NM m,mrd,tail;
1.7       noro      921:        NM mul;
1.10      noro      922:        int n,sugar,psugar,sugar0,stat,index;
1.6       noro      923:        int c,c1,c2;
1.17      noro      924:        RHist h;
1.11      noro      925:        NDV p,red;
1.16      noro      926:        Q cg,cred,gcd;
1.21      noro      927:        double hmag;
1.1       noro      928:
                    929:        if ( !g ) {
                    930:                *rp = 0;
                    931:                return 1;
                    932:        }
1.21      noro      933:        if ( !mod )
                    934:                hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
                    935:
1.14      noro      936:        sugar0 = sugar = SG(g);
1.1       noro      937:        n = NV(g);
1.7       noro      938:        mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
1.1       noro      939:        for ( d = 0; g; ) {
1.6       noro      940:                index = nd_find_reducer(g);
                    941:                if ( index >= 0 ) {
1.17      noro      942:                        h = nd_psh[index];
                    943:                        ndl_sub(HDL(g),DL(h),DL(mul));
                    944:                        TD(mul) = HTD(g)-TD(h);
1.10      noro      945: #if 0
1.14      noro      946:                        if ( d && (SG(p)+TD(mul)) > sugar ) {
1.10      noro      947:                                goto afo;
                    948:                        }
                    949: #endif
1.14      noro      950:                        if ( ndl_check_bound2(index,DL(mul)) ) {
1.6       noro      951:                                nd_free(g); nd_free(d);
                    952:                                return 0;
                    953:                        }
1.16      noro      954:                        if ( mod ) {
1.17      noro      955:                                p = nd_ps[index];
1.19      noro      956:                                c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
                    957:                                DMAR(c1,c2,0,mod,c); CM(mul) = c;
1.16      noro      958:                        } else {
1.20      noro      959:                                p = nd_psq[index];
1.17      noro      960:                                igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);
1.16      noro      961:                                chsgnq(cg,&CQ(mul));
1.20      noro      962:                                nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);
1.16      noro      963:                        }
1.20      noro      964:                        ndv_mul_nm(mod,p,mul,ndv_red);
                    965:                        g = ndv_add(mod,g,ndv_red);
1.14      noro      966:                        sugar = MAX(sugar,SG(ndv_red));
1.22      noro      967:                        if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {
1.21      noro      968:                                nd_removecont2(d,g);
                    969:                                hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
                    970:                        }
1.1       noro      971:                } else if ( !full ) {
                    972:                        *rp = g;
                    973:                        return 1;
                    974:                } else {
1.10      noro      975: afo:
1.1       noro      976:                        m = BDY(g);
                    977:                        if ( NEXT(m) ) {
                    978:                                BDY(g) = NEXT(m); NEXT(m) = 0;
                    979:                        } else {
                    980:                                FREEND(g); g = 0;
                    981:                        }
                    982:                        if ( d ) {
                    983:                                NEXT(tail)=m;
                    984:                                tail=m;
                    985:                        } else {
                    986:                                MKND(n,m,d);
                    987:                                tail = BDY(d);
                    988:                        }
                    989:                }
                    990:        }
                    991:        if ( d )
1.14      noro      992:                SG(d) = sugar;
1.1       noro      993:        *rp = d;
                    994:        return 1;
                    995: }
1.23      noro      996:
                    997: int nd_nf_direct(int mod,ND g,NDV *ps,int len,int full,ND *rp)
                    998: {
                    999:        ND d;
                   1000:        NM m,mrd,tail;
                   1001:        NM mul;
                   1002:        int n,sugar,psugar,sugar0,stat,index;
                   1003:        int c,c1,c2;
                   1004:        RHist h;
                   1005:        NDV p,red;
                   1006:        Q cg,cred,gcd;
                   1007:        double hmag;
                   1008:
                   1009:        if ( !g ) {
                   1010:                *rp = 0;
                   1011:                return 1;
                   1012:        }
                   1013: #if 0
                   1014:        if ( !mod )
                   1015:                hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
                   1016: #else
                   1017:        /* XXX */
                   1018:        hmag = 0;
                   1019: #endif
                   1020:
                   1021:        sugar0 = sugar = SG(g);
                   1022:        n = NV(g);
                   1023:        mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
                   1024:        for ( d = 0; g; ) {
                   1025:                index = nd_find_reducer_direct(g,ps,len);
                   1026:                if ( index >= 0 ) {
                   1027:                        p = ps[index];
                   1028:                        ndl_sub(HDL(g),HDL(p),DL(mul));
                   1029:                        TD(mul) = HTD(g)-HTD(p);
                   1030:                        if ( ndl_check_bound2_direct(HDL(p),DL(mul)) ) {
                   1031:                                nd_free(g); nd_free(d);
                   1032:                                return 0;
                   1033:                        }
                   1034:                        if ( mod ) {
                   1035:                                c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
                   1036:                                DMAR(c1,c2,0,mod,c); CM(mul) = c;
                   1037:                        } else {
                   1038:                                igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);
                   1039:                                chsgnq(cg,&CQ(mul));
                   1040:                                nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);
                   1041:                        }
                   1042:                        ndv_mul_nm(mod,p,mul,ndv_red);
                   1043:                        g = ndv_add(mod,g,ndv_red);
                   1044:                        sugar = MAX(sugar,SG(ndv_red));
                   1045:                        if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {
                   1046:                                nd_removecont2(d,g);
                   1047:                                hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
                   1048:                        }
                   1049:                } else if ( !full ) {
                   1050:                        *rp = g;
                   1051:                        return 1;
                   1052:                } else {
                   1053:                        m = BDY(g);
                   1054:                        if ( NEXT(m) ) {
                   1055:                                BDY(g) = NEXT(m); NEXT(m) = 0;
                   1056:                        } else {
                   1057:                                FREEND(g); g = 0;
                   1058:                        }
                   1059:                        if ( d ) {
                   1060:                                NEXT(tail)=m;
                   1061:                                tail=m;
                   1062:                        } else {
                   1063:                                MKND(n,m,d);
                   1064:                                tail = BDY(d);
                   1065:                        }
                   1066:                }
                   1067:        }
                   1068:        if ( d )
                   1069:                SG(d) = sugar;
                   1070:        *rp = d;
                   1071:        return 1;
                   1072: }
1.1       noro     1073: #else
                   1074:
                   1075: ND nd_remove_head(ND p)
                   1076: {
                   1077:        NM m;
                   1078:
                   1079:        m = BDY(p);
                   1080:        if ( !NEXT(m) ) {
                   1081:                FREEND(p);
                   1082:                p = 0;
                   1083:        } else
                   1084:                BDY(p) = NEXT(m);
                   1085:        FREENM(m);
                   1086:        return p;
                   1087: }
                   1088:
                   1089: PGeoBucket create_pbucket()
                   1090: {
                   1091:     PGeoBucket g;
                   1092:
                   1093:        g = CALLOC(1,sizeof(struct oPGeoBucket));
                   1094:        g->m = -1;
                   1095:        return g;
                   1096: }
                   1097:
1.19      noro     1098: void add_pbucket(int mod,PGeoBucket g,ND d)
1.1       noro     1099: {
                   1100:        int l,k,m;
                   1101:
                   1102:        l = nd_length(d);
                   1103:        for ( k = 0, m = 1; l > m; k++, m <<= 2 );
                   1104:        /* 4^(k-1) < l <= 4^k */
1.19      noro     1105:        d = nd_add(mod,g->body[k],d);
1.1       noro     1106:        for ( ; d && nd_length(d) > 1<<(2*k); k++ ) {
                   1107:                g->body[k] = 0;
1.19      noro     1108:                d = nd_add(int mod,g->body[k+1],d);
1.1       noro     1109:        }
                   1110:        g->body[k] = d;
                   1111:        g->m = MAX(g->m,k);
                   1112: }
                   1113:
1.19      noro     1114: int head_pbucket(int mod,PGeoBucket g)
1.1       noro     1115: {
                   1116:        int j,i,c,k,nv,sum;
                   1117:        unsigned int *di,*dj;
                   1118:        ND gi,gj;
                   1119:
                   1120:        k = g->m;
                   1121:        while ( 1 ) {
                   1122:                j = -1;
                   1123:                for ( i = 0; i <= k; i++ ) {
                   1124:                        if ( !(gi = g->body[i]) )
                   1125:                                continue;
                   1126:                        if ( j < 0 ) {
                   1127:                                j = i;
                   1128:                                gj = g->body[j];
                   1129:                                dj = HDL(gj);
1.14      noro     1130:                                sum = HCM(gj);
1.1       noro     1131:                        } else {
                   1132:                                di = HDL(gi);
                   1133:                                nv = NV(gi);
                   1134:                                if ( HTD(gi) > HTD(gj) )
                   1135:                                        c = 1;
                   1136:                                else if ( HTD(gi) < HTD(gj) )
                   1137:                                        c = -1;
                   1138:                                else
                   1139:                                        c = ndl_compare(di,dj);
                   1140:                                if ( c > 0 ) {
                   1141:                                        if ( sum )
1.14      noro     1142:                                                HCM(gj) = sum;
1.1       noro     1143:                                        else
                   1144:                                                g->body[j] = nd_remove_head(gj);
                   1145:                                        j = i;
                   1146:                                        gj = g->body[j];
                   1147:                                        dj = HDL(gj);
1.14      noro     1148:                                        sum = HCM(gj);
1.1       noro     1149:                                } else if ( c == 0 ) {
1.19      noro     1150:                                        sum = sum+HCM(gi)-mod;
1.1       noro     1151:                                        if ( sum < 0 )
1.19      noro     1152:                                                sum += mod;
1.1       noro     1153:                                        g->body[i] = nd_remove_head(gi);
                   1154:                                }
                   1155:                        }
                   1156:                }
                   1157:                if ( j < 0 )
                   1158:                        return -1;
                   1159:                else if ( sum ) {
1.14      noro     1160:                        HCM(gj) = sum;
1.1       noro     1161:                        return j;
                   1162:                } else
                   1163:                        g->body[j] = nd_remove_head(gj);
                   1164:        }
                   1165: }
                   1166:
                   1167: ND normalize_pbucket(PGeoBucket g)
                   1168: {
                   1169:        int i;
                   1170:        ND r,t;
                   1171:
                   1172:        r = 0;
                   1173:        for ( i = 0; i <= g->m; i++ )
1.19      noro     1174:                r = nd_add(mod,r,g->body[i]);
1.1       noro     1175:        return r;
                   1176: }
                   1177:
                   1178: ND nd_nf(ND g,int full)
                   1179: {
                   1180:        ND u,p,d,red;
                   1181:        NODE l;
                   1182:        NM m,mrd;
                   1183:        int sugar,psugar,n,h_reducible,h;
                   1184:        PGeoBucket bucket;
                   1185:
                   1186:        if ( !g ) {
                   1187:                return 0;
                   1188:        }
1.14      noro     1189:        sugar = SG(g);
                   1190:        n = NV(g);
1.1       noro     1191:        bucket = create_pbucket();
                   1192:        add_pbucket(bucket,g);
                   1193:        d = 0;
                   1194:        while ( 1 ) {
                   1195:                h = head_pbucket(bucket);
                   1196:                if ( h < 0 ) {
                   1197:                        if ( d )
1.14      noro     1198:                                SG(d) = sugar;
1.1       noro     1199:                        return d;
                   1200:                }
                   1201:                g = bucket->body[h];
                   1202:                red = nd_find_reducer(g);
                   1203:                if ( red ) {
                   1204:                        bucket->body[h] = nd_remove_head(g);
                   1205:                        red = nd_remove_head(red);
                   1206:                        add_pbucket(bucket,red);
1.14      noro     1207:                        sugar = MAX(sugar,SG(red));
1.1       noro     1208:                } else if ( !full ) {
                   1209:                        g = normalize_pbucket(bucket);
                   1210:                        if ( g )
1.14      noro     1211:                                SG(g) = sugar;
1.1       noro     1212:                        return g;
                   1213:                } else {
                   1214:                        m = BDY(g);
                   1215:                        if ( NEXT(m) ) {
                   1216:                                BDY(g) = NEXT(m); NEXT(m) = 0;
                   1217:                        } else {
                   1218:                                FREEND(g); g = 0;
                   1219:                        }
                   1220:                        bucket->body[h] = g;
                   1221:                        NEXT(m) = 0;
                   1222:                        if ( d ) {
                   1223:                                for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
                   1224:                                NEXT(mrd) = m;
                   1225:                        } else {
                   1226:                                MKND(n,m,d);
                   1227:                        }
                   1228:                }
                   1229:        }
                   1230: }
                   1231: #endif
                   1232:
1.23      noro     1233: NODE nd_gb(int m)
1.1       noro     1234: {
                   1235:        int i,nh,sugar,stat;
1.23      noro     1236:        NODE r,g,t;
1.1       noro     1237:        ND_pairs d;
                   1238:        ND_pairs l;
                   1239:        ND h,nf;
                   1240:
1.23      noro     1241:        g = 0; d = 0;
                   1242:        for ( i = 0; i < nd_psn; i++ ) {
1.1       noro     1243:                d = update_pairs(d,g,i);
                   1244:                g = update_base(g,i);
                   1245:        }
                   1246:        sugar = 0;
                   1247:        while ( d ) {
                   1248: again:
                   1249:                l = nd_minp(d,&d);
1.14      noro     1250:                if ( SG(l) != sugar ) {
                   1251:                        sugar = SG(l);
1.1       noro     1252:                        fprintf(asir_out,"%d",sugar);
                   1253:                }
1.16      noro     1254:                stat = nd_sp(m,l,&h);
1.1       noro     1255:                if ( !stat ) {
                   1256:                        NEXT(l) = d; d = l;
1.20      noro     1257:                        d = nd_reconstruct(m,0,d);
1.1       noro     1258:                        goto again;
                   1259:                }
1.16      noro     1260:                stat = nd_nf(m,h,!Top,&nf);
1.1       noro     1261:                if ( !stat ) {
                   1262:                        NEXT(l) = d; d = l;
1.20      noro     1263:                        d = nd_reconstruct(m,0,d);
1.1       noro     1264:                        goto again;
                   1265:                } else if ( nf ) {
                   1266:                        printf("+"); fflush(stdout);
1.16      noro     1267:                        nh = nd_newps(m,nf);
1.1       noro     1268:                        d = update_pairs(d,g,nh);
                   1269:                        g = update_base(g,nh);
                   1270:                        FREENDP(l);
                   1271:                } else {
                   1272:                        printf("."); fflush(stdout);
                   1273:                        FREENDP(l);
                   1274:                }
                   1275:        }
1.23      noro     1276:        if ( m )
                   1277:                for ( t = g; t; t = NEXT(t) )
                   1278:                        BDY(t) = (pointer)nd_ps[(int)BDY(t)];
                   1279:        else
                   1280:                for ( t = g; t; t = NEXT(t) )
                   1281:                        BDY(t) = (pointer)nd_psq[(int)BDY(t)];
1.1       noro     1282:        return g;
                   1283: }
                   1284:
1.23      noro     1285: NODE nd_gb_trace(int m)
1.20      noro     1286: {
                   1287:        int i,nh,sugar,stat;
1.23      noro     1288:        NODE r,g,t;
1.20      noro     1289:        ND_pairs d;
                   1290:        ND_pairs l;
                   1291:        ND h,nf,nfq;
                   1292:
1.23      noro     1293:        g = 0; d = 0;
                   1294:        for ( i = 0; i < nd_psn; i++ ) {
1.20      noro     1295:                d = update_pairs(d,g,i);
                   1296:                g = update_base(g,i);
                   1297:        }
                   1298:        sugar = 0;
                   1299:        while ( d ) {
                   1300: again:
                   1301:                l = nd_minp(d,&d);
                   1302:                if ( SG(l) != sugar ) {
                   1303:                        sugar = SG(l);
                   1304:                        fprintf(asir_out,"%d",sugar);
                   1305:                }
                   1306:                stat = nd_sp(m,l,&h);
                   1307:                if ( !stat ) {
                   1308:                        NEXT(l) = d; d = l;
                   1309:                        d = nd_reconstruct(m,1,d);
                   1310:                        goto again;
                   1311:                }
                   1312:                stat = nd_nf(m,h,!Top,&nf);
                   1313:                if ( !stat ) {
                   1314:                        NEXT(l) = d; d = l;
                   1315:                        d = nd_reconstruct(m,1,d);
                   1316:                        goto again;
                   1317:                } else if ( nf ) {
                   1318:                        /* overflow does not occur */
                   1319:                        nd_sp(0,l,&h);
                   1320:                        nd_nf(0,h,!Top,&nfq);
                   1321:                        if ( nfq ) {
                   1322:                                printf("+"); fflush(stdout);
                   1323:                                nh = nd_newps_trace(m,nf,nfq);
                   1324:                                d = update_pairs(d,g,nh);
                   1325:                                g = update_base(g,nh);
                   1326:                        } else {
                   1327:                                printf("*"); fflush(stdout);
                   1328:                        }
                   1329:                } else {
                   1330:                        printf("."); fflush(stdout);
                   1331:                }
                   1332:                FREENDP(l);
                   1333:        }
1.23      noro     1334:        for ( t = g; t; t = NEXT(t) )
                   1335:                BDY(t) = (pointer)nd_psq[(int)BDY(t)];
1.20      noro     1336:        return g;
                   1337: }
                   1338:
1.23      noro     1339: int ndv_compare(NDV *p1,NDV *p2)
                   1340: {
                   1341:        int td1,td2;
                   1342:
                   1343:        td1 = HTD(*p1); td2 = HTD(*p2);
                   1344:        if ( td1 > td2 ) return 1;
                   1345:        else if ( td1 < td2 ) return -1;
                   1346:        else return ndl_compare(HDL(*p1),HDL(*p2));
                   1347: }
                   1348:
                   1349: int ndv_compare_rev(NDV *p1,NDV *p2)
                   1350: {
                   1351:        return -ndv_compare(p1,p2);
                   1352: }
                   1353:
                   1354: NODE nd_reduceall(int m,NODE f)
                   1355: {
                   1356:        int i,j,n,stat;
                   1357:        NDV *w,*ps;
                   1358:        ND nf,g;
                   1359:        NODE t,a0,a;
                   1360:
                   1361:        for ( n = 0, t = f; t; t = NEXT(t), n++ );
                   1362:        ps = (NDV *)ALLOCA(n*sizeof(NDV));
                   1363:        for ( i = 0, t = f; i < n; i++, t = NEXT(t) )
                   1364:                ps[i] = (NDV)BDY(t);
                   1365:        qsort(ps,n,sizeof(NDV),(int (*)(const void *,const void *))ndv_compare);
                   1366:        w = (NDV *)ALLOCA((n-1)*sizeof(NDV));
                   1367:        for ( i = 0; i < n; i++ ) {
                   1368:                for ( j = 0; j < i; j++ )
                   1369:                        w[j] = (NDV)ps[j];
                   1370:                for ( j = i+1; j < n; j++ )
                   1371:                        w[j-1] = ps[j];
                   1372:                g = ndvtond(m,ps[i]);
                   1373:                stat = nd_nf_direct(m,g,w,n-1,1,&nf);
                   1374:                if ( !stat )
                   1375:                        nd_reconstruct_direct(m,ps,n);
                   1376:                else if ( !nf ) {
                   1377:                        printf("."); fflush(stdout);
                   1378:                        ndv_free(ps[i]);
                   1379:                        for ( j = i+1; j < n; j++ )
                   1380:                                ps[j-1] = ps[j];
                   1381:                        n--;
                   1382:                } else {
                   1383:                        printf("."); fflush(stdout);
                   1384:                        ndv_free(ps[i]);
1.24    ! noro     1385:                        nd_removecont(m,nf);
1.23      noro     1386:                        ps[i] = ndtondv(m,nf);
                   1387:                        nd_free(nf);
                   1388:                }
                   1389:        }
                   1390:        for ( a0 = 0, i = 0; i < n; i++ ) {
                   1391:                NEXTNODE(a0,a);
                   1392:                BDY(a) = (pointer)ps[i];
                   1393:        }
                   1394:        NEXT(a) = 0;
                   1395:        return a0;
                   1396: }
                   1397:
1.1       noro     1398: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t)
                   1399: {
                   1400:        ND_pairs d1,nd,cur,head,prev,remove;
                   1401:
                   1402:        if ( !g ) return d;
                   1403:        d = crit_B(d,t);
                   1404:        d1 = nd_newpairs(g,t);
                   1405:        d1 = crit_M(d1);
                   1406:        d1 = crit_F(d1);
                   1407:        prev = 0; cur = head = d1;
                   1408:        while ( cur ) {
                   1409:                if ( crit_2( cur->i1,cur->i2 ) ) {
                   1410:                        remove = cur;
                   1411:                        if ( !prev ) {
                   1412:                                head = cur = NEXT(cur);
                   1413:                        } else {
                   1414:                                cur = NEXT(prev) = NEXT(cur);
                   1415:                        }
                   1416:                        FREENDP(remove);
                   1417:                } else {
                   1418:                        prev = cur;
                   1419:                        cur = NEXT(cur);
                   1420:                }
                   1421:        }
                   1422:        if ( !d )
                   1423:                return head;
                   1424:        else {
                   1425:                nd = d;
                   1426:                while ( NEXT(nd) )
                   1427:                        nd = NEXT(nd);
                   1428:                NEXT(nd) = head;
                   1429:                return d;
                   1430:        }
                   1431: }
                   1432:
                   1433: ND_pairs nd_newpairs( NODE g, int t )
                   1434: {
                   1435:        NODE h;
                   1436:        unsigned int *dl;
                   1437:        int td,ts,s;
                   1438:        ND_pairs r,r0;
                   1439:
1.20      noro     1440:        dl = DL(nd_psh[t]);
                   1441:        td = TD(nd_psh[t]);
                   1442:        ts = SG(nd_psh[t]) - td;
1.1       noro     1443:        for ( r0 = 0, h = g; h; h = NEXT(h) ) {
                   1444:                NEXTND_pairs(r0,r);
                   1445:                r->i1 = (int)BDY(h);
                   1446:                r->i2 = t;
1.20      noro     1447:                ndl_lcm(DL(nd_psh[r->i1]),dl,r->lcm);
1.14      noro     1448:                TD(r) = ndl_td(r->lcm);
1.20      noro     1449:                s = SG(nd_psh[r->i1])-TD(nd_psh[r->i1]);
1.14      noro     1450:                SG(r) = MAX(s,ts) + TD(r);
1.1       noro     1451:        }
                   1452:        NEXT(r) = 0;
                   1453:        return r0;
                   1454: }
                   1455:
                   1456: ND_pairs crit_B( ND_pairs d, int s )
                   1457: {
                   1458:        ND_pairs cur,head,prev,remove;
                   1459:        unsigned int *t,*tl,*lcm;
                   1460:        int td,tdl;
                   1461:
                   1462:        if ( !d ) return 0;
1.20      noro     1463:        t = DL(nd_psh[s]);
1.1       noro     1464:        prev = 0;
                   1465:        head = cur = d;
                   1466:        lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                   1467:        while ( cur ) {
                   1468:                tl = cur->lcm;
                   1469:                if ( ndl_reducible(tl,t)
1.20      noro     1470:                        && (ndl_lcm(DL(nd_psh[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
                   1471:                        && (ndl_lcm(DL(nd_psh[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
1.1       noro     1472:                        remove = cur;
                   1473:                        if ( !prev ) {
                   1474:                                head = cur = NEXT(cur);
                   1475:                        } else {
                   1476:                                cur = NEXT(prev) = NEXT(cur);
                   1477:                        }
                   1478:                        FREENDP(remove);
                   1479:                } else {
                   1480:                        prev = cur;
                   1481:                        cur = NEXT(cur);
                   1482:                }
                   1483:        }
                   1484:        return head;
                   1485: }
                   1486:
                   1487: ND_pairs crit_M( ND_pairs d1 )
                   1488: {
                   1489:        ND_pairs e,d2,d3,dd,p;
                   1490:        unsigned int *id,*jd;
                   1491:        int itd,jtd;
                   1492:
                   1493:        for ( dd = 0, e = d1; e; e = d3 ) {
                   1494:                if ( !(d2 = NEXT(e)) ) {
                   1495:                        NEXT(e) = dd;
                   1496:                        return e;
                   1497:                }
                   1498:                id = e->lcm;
1.14      noro     1499:                itd = TD(e);
1.1       noro     1500:                for ( d3 = 0; d2; d2 = p ) {
                   1501:                        p = NEXT(d2),
                   1502:                        jd = d2->lcm;
1.14      noro     1503:                        jtd = TD(d2);
1.1       noro     1504:                        if ( jtd == itd  )
                   1505:                                if ( id == jd );
                   1506:                                else if ( ndl_reducible(jd,id) ) continue;
                   1507:                                else if ( ndl_reducible(id,jd) ) goto delit;
                   1508:                                else ;
                   1509:                        else if ( jtd > itd )
                   1510:                                if ( ndl_reducible(jd,id) ) continue;
                   1511:                                else ;
                   1512:                        else if ( ndl_reducible(id,jd ) ) goto delit;
                   1513:                        NEXT(d2) = d3;
                   1514:                        d3 = d2;
                   1515:                }
                   1516:                NEXT(e) = dd;
                   1517:                dd = e;
                   1518:                continue;
                   1519:                /**/
                   1520:        delit:  NEXT(d2) = d3;
                   1521:                d3 = d2;
                   1522:                for ( ; p; p = d2 ) {
                   1523:                        d2 = NEXT(p);
                   1524:                        NEXT(p) = d3;
                   1525:                        d3 = p;
                   1526:                }
                   1527:                FREENDP(e);
                   1528:        }
                   1529:        return dd;
                   1530: }
                   1531:
                   1532: ND_pairs crit_F( ND_pairs d1 )
                   1533: {
                   1534:        ND_pairs rest, head,remove;
                   1535:        ND_pairs last, p, r, w;
                   1536:        int s;
                   1537:
                   1538:        for ( head = last = 0, p = d1; NEXT(p); ) {
                   1539:                r = w = equivalent_pairs(p,&rest);
1.14      noro     1540:                s = SG(r);
1.1       noro     1541:                w = NEXT(w);
                   1542:                while ( w ) {
                   1543:                        if ( crit_2(w->i1,w->i2) ) {
                   1544:                                r = w;
                   1545:                                w = NEXT(w);
                   1546:                                while ( w ) {
                   1547:                                        remove = w;
                   1548:                                        w = NEXT(w);
                   1549:                                        FREENDP(remove);
                   1550:                                }
                   1551:                                break;
1.14      noro     1552:                        } else if ( SG(w) < s ) {
1.1       noro     1553:                                FREENDP(r);
                   1554:                                r = w;
1.14      noro     1555:                                s = SG(r);
1.1       noro     1556:                                w = NEXT(w);
                   1557:                        } else {
                   1558:                                remove = w;
                   1559:                                w = NEXT(w);
                   1560:                                FREENDP(remove);
                   1561:                        }
                   1562:                }
                   1563:                if ( last ) NEXT(last) = r;
                   1564:                else head = r;
                   1565:                NEXT(last = r) = 0;
                   1566:                p = rest;
                   1567:                if ( !p ) return head;
                   1568:        }
                   1569:        if ( !last ) return p;
                   1570:        NEXT(last) = p;
                   1571:        return head;
                   1572: }
                   1573:
                   1574: int crit_2( int dp1, int dp2 )
                   1575: {
1.20      noro     1576:        return ndl_disjoint(DL(nd_psh[dp1]),DL(nd_psh[dp2]));
1.1       noro     1577: }
                   1578:
                   1579: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest )
                   1580: {
                   1581:        ND_pairs w,p,r,s;
                   1582:        unsigned int *d;
                   1583:        int td;
                   1584:
                   1585:        w = d1;
                   1586:        d = w->lcm;
1.14      noro     1587:        td = TD(w);
1.1       noro     1588:        s = NEXT(w);
                   1589:        NEXT(w) = 0;
                   1590:        for ( r = 0; s; s = p ) {
                   1591:                p = NEXT(s);
1.14      noro     1592:                if ( td == TD(s) && ndl_equal(d,s->lcm) ) {
1.1       noro     1593:                        NEXT(s) = w;
                   1594:                        w = s;
                   1595:                } else {
                   1596:                        NEXT(s) = r;
                   1597:                        r = s;
                   1598:                }
                   1599:        }
                   1600:        *prest = r;
                   1601:        return w;
                   1602: }
                   1603:
                   1604: NODE update_base(NODE nd,int ndp)
                   1605: {
                   1606:        unsigned int *dl, *dln;
                   1607:        NODE last, p, head;
                   1608:        int td,tdn;
                   1609:
1.20      noro     1610:        dl = DL(nd_psh[ndp]);
                   1611:        td = TD(nd_psh[ndp]);
1.1       noro     1612:        for ( head = last = 0, p = nd; p; ) {
1.20      noro     1613:                dln = DL(nd_psh[(int)BDY(p)]);
                   1614:                tdn = TD(nd_psh[(int)BDY(p)]);
1.1       noro     1615:                if ( tdn >= td && ndl_reducible( dln, dl ) ) {
                   1616:                        p = NEXT(p);
                   1617:                        if ( last ) NEXT(last) = p;
                   1618:                } else {
                   1619:                        if ( !last ) head = p;
                   1620:                        p = NEXT(last = p);
                   1621:                }
                   1622:        }
                   1623:        head = append_one(head,ndp);
                   1624:        return head;
                   1625: }
                   1626:
                   1627: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
                   1628: {
                   1629:        ND_pairs m,ml,p,l;
                   1630:        unsigned int *lcm;
                   1631:        int s,td,len,tlen,c;
                   1632:
                   1633:        if ( !(p = NEXT(m = d)) ) {
                   1634:                *prest = p;
                   1635:                NEXT(m) = 0;
                   1636:                return m;
                   1637:        }
                   1638:        lcm = m->lcm;
1.14      noro     1639:        s = SG(m);
                   1640:        td = TD(m);
1.6       noro     1641:        len = nd_psl[m->i1]+nd_psl[m->i2];
1.1       noro     1642:        for ( ml = 0, l = m; p; p = NEXT(l = p) ) {
1.14      noro     1643:                if (SG(p) < s)
1.1       noro     1644:                        goto find;
1.14      noro     1645:                else if ( SG(p) == s ) {
                   1646:                        if ( TD(p) < td )
1.1       noro     1647:                                goto find;
1.14      noro     1648:                        else if ( TD(p) == td ) {
1.1       noro     1649:                                c = ndl_compare(p->lcm,lcm);
                   1650:                                if ( c < 0 )
                   1651:                                        goto find;
1.10      noro     1652: #if 0
1.1       noro     1653:                                else if ( c == 0 ) {
1.6       noro     1654:                                        tlen = nd_psl[p->i1]+nd_psl[p->i2];
1.1       noro     1655:                                        if ( tlen < len )
                   1656:                                                goto find;
                   1657:                                }
1.10      noro     1658: #endif
1.1       noro     1659:                        }
                   1660:                }
                   1661:                continue;
                   1662: find:
                   1663:                ml = l;
                   1664:                m = p;
                   1665:                lcm = m->lcm;
1.14      noro     1666:                s = SG(m);
                   1667:                td = TD(m);
1.1       noro     1668:                len = tlen;
                   1669:        }
                   1670:        if ( !ml ) *prest = NEXT(m);
                   1671:        else {
                   1672:                NEXT(ml) = NEXT(m);
                   1673:                *prest = d;
                   1674:        }
                   1675:        NEXT(m) = 0;
                   1676:        return m;
                   1677: }
                   1678:
1.16      noro     1679: int nd_newps(int mod,ND a)
1.1       noro     1680: {
1.3       noro     1681:        int len;
1.13      noro     1682:        RHist r;
1.20      noro     1683:        NDV b;
1.3       noro     1684:
1.1       noro     1685:        if ( nd_psn == nd_pslen ) {
                   1686:                nd_pslen *= 2;
1.6       noro     1687:                nd_psl = (int *)REALLOC((char *)nd_psl,nd_pslen*sizeof(int));
1.11      noro     1688:                nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV));
1.20      noro     1689:                nd_psq = (NDV *)REALLOC((char *)nd_psq,nd_pslen*sizeof(NDV));
1.13      noro     1690:                nd_psh = (RHist *)REALLOC((char *)nd_psh,nd_pslen*sizeof(RHist));
1.1       noro     1691:                nd_bound = (unsigned int **)
                   1692:                        REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
                   1693:        }
1.20      noro     1694:        nd_removecont(mod,a);
                   1695:        nd_bound[nd_psn] = nd_compute_bound(a);
                   1696:        NEWRHist(r); SG(r) = SG(a); TD(r) = HTD(a); ndl_copy(HDL(a),DL(r));
                   1697:        nd_psh[nd_psn] = r;
                   1698:        b = ndtondv(mod,a);
                   1699:        len = LEN(b);
1.16      noro     1700:        if ( mod )
1.20      noro     1701:                nd_ps[nd_psn] = b;
1.16      noro     1702:        else
1.20      noro     1703:                nd_psq[nd_psn] = b;
                   1704:        nd_psl[nd_psn] = len;
1.13      noro     1705:        nd_free(a);
1.20      noro     1706:        if ( len > nmv_len ) {
                   1707:                nmv_len = 2*len;
                   1708:                BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
                   1709:        }
                   1710:        return nd_psn++;
                   1711: }
                   1712:
                   1713: int nd_newps_trace(int mod,ND nf,ND nfq)
                   1714: {
                   1715:        int len;
                   1716:        RHist r;
                   1717:        NDV b;
                   1718:
                   1719:        if ( nd_psn == nd_pslen ) {
                   1720:                nd_pslen *= 2;
                   1721:                nd_psl = (int *)REALLOC((char *)nd_psl,nd_pslen*sizeof(int));
                   1722:                nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV));
                   1723:                nd_psq = (NDV *)REALLOC((char *)nd_psq,nd_pslen*sizeof(NDV));
                   1724:                nd_psh = (RHist *)REALLOC((char *)nd_psh,nd_pslen*sizeof(RHist));
                   1725:                nd_bound = (unsigned int **)
                   1726:                        REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
                   1727:        }
                   1728:        nd_removecont(mod,nf);
                   1729:        nd_ps[nd_psn] = ndtondv(mod,nf);
                   1730:
                   1731:        nd_removecont(0,nfq);
                   1732:        nd_psq[nd_psn] = ndtondv(0,nfq);
                   1733:
                   1734:        nd_bound[nd_psn] = nd_compute_bound(nfq);
                   1735:        NEWRHist(r); SG(r) = SG(nf); TD(r) = HTD(nf); ndl_copy(HDL(nf),DL(r));
                   1736:        nd_psh[nd_psn] = r;
                   1737:
                   1738:        len = LEN(nd_psq[nd_psn]);
                   1739:        nd_psl[nd_psn] = len;
                   1740:
                   1741:        nd_free(nf); nd_free(nfq);
1.3       noro     1742:        if ( len > nmv_len ) {
                   1743:                nmv_len = 2*len;
                   1744:                BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
                   1745:        }
1.1       noro     1746:        return nd_psn++;
                   1747: }
                   1748:
1.23      noro     1749: void nd_setup(int mod,NODE f)
1.1       noro     1750: {
1.5       noro     1751:        int i,j,td,len,max;
1.1       noro     1752:        NODE s,s0,f0;
1.5       noro     1753:        unsigned int *d;
1.13      noro     1754:        RHist r;
1.20      noro     1755:        NDV a;
1.11      noro     1756:
                   1757:        nd_found = 0; nd_notfirst = 0; nd_create = 0;
1.1       noro     1758:
                   1759:        nd_psn = length(f); nd_pslen = 2*nd_psn;
1.6       noro     1760:        nd_psl = (int *)MALLOC(nd_pslen*sizeof(int));
1.11      noro     1761:        nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
1.20      noro     1762:        nd_psq = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
1.13      noro     1763:        nd_psh = (RHist *)MALLOC(nd_pslen*sizeof(RHist));
1.1       noro     1764:        nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
1.5       noro     1765:        for ( max = 0, i = 0, s = f; i < nd_psn; i++, s = NEXT(s) ) {
                   1766:                nd_bound[i] = d = dp_compute_bound((DP)BDY(s));
                   1767:                for ( j = 0; j < nd_nvar; j++ )
                   1768:                        max = MAX(d[j],max);
                   1769:        }
1.11      noro     1770:        if ( !nd_red )
1.13      noro     1771:                nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));
                   1772:        bzero(nd_red,REDTAB_LEN*sizeof(RHist));
1.5       noro     1773:
                   1774:        if ( max < 2 )
                   1775:                nd_bpe = 2;
                   1776:        else if ( max < 4 )
                   1777:                nd_bpe = 4;
                   1778:        else if ( max < 64 )
                   1779:                nd_bpe = 6;
                   1780:        else if ( max < 256 )
                   1781:                nd_bpe = 8;
                   1782:        else if ( max < 65536 )
                   1783:                nd_bpe = 16;
                   1784:        else
                   1785:                nd_bpe = 32;
1.13      noro     1786:
1.1       noro     1787:        nd_setup_parameters();
                   1788:        nd_free_private_storage();
1.3       noro     1789:        len = 0;
1.1       noro     1790:        for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
1.20      noro     1791:                NEWRHist(r);
                   1792:                a = dptondv(mod,(DP)BDY(f));
                   1793:                ndv_removecont(mod,a);
                   1794:                len = MAX(len,LEN(a));
                   1795:                SG(r) = HTD(a); TD(r) = HTD(a); ndl_copy(HDL(a),DL(r));
1.16      noro     1796:                if ( mod )
1.20      noro     1797:                        nd_ps[i] = a;
1.16      noro     1798:                else
1.20      noro     1799:                        nd_psq[i] = a;
1.13      noro     1800:                nd_psh[i] = r;
1.1       noro     1801:        }
1.3       noro     1802:        nmv_len = 16*len;
                   1803:        NEWNDV(ndv_red);
1.16      noro     1804:        if ( mod )
                   1805:                BDY(ndv_red) = (NMV)MALLOC_ATOMIC(nmv_len*nmv_adv);
                   1806:        else
                   1807:                BDY(ndv_red) = (NMV)MALLOC(nmv_len*nmv_adv);
1.1       noro     1808: }
                   1809:
1.23      noro     1810: void nd_setup_trace(int mod,NODE f)
1.20      noro     1811: {
                   1812:        int i,j,td,len,max;
                   1813:        NODE s,s0,f0;
                   1814:        unsigned int *d;
                   1815:        RHist r;
                   1816:        NDV a;
                   1817:
                   1818:        nd_found = 0; nd_notfirst = 0; nd_create = 0;
                   1819:
                   1820:        nd_psn = length(f); nd_pslen = 2*nd_psn;
                   1821:        nd_psl = (int *)MALLOC(nd_pslen*sizeof(int));
                   1822:        nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
                   1823:        nd_psq = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
                   1824:        nd_psh = (RHist *)MALLOC(nd_pslen*sizeof(RHist));
                   1825:        nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
                   1826:        for ( max = 0, i = 0, s = f; i < nd_psn; i++, s = NEXT(s) ) {
                   1827:                nd_bound[i] = d = dp_compute_bound((DP)BDY(s));
                   1828:                for ( j = 0; j < nd_nvar; j++ )
                   1829:                        max = MAX(d[j],max);
                   1830:        }
                   1831:        if ( !nd_red )
                   1832:                nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));
                   1833:        bzero(nd_red,REDTAB_LEN*sizeof(RHist));
                   1834:
                   1835:        if ( max < 2 )
                   1836:                nd_bpe = 2;
                   1837:        else if ( max < 4 )
                   1838:                nd_bpe = 4;
                   1839:        else if ( max < 64 )
                   1840:                nd_bpe = 6;
                   1841:        else if ( max < 256 )
                   1842:                nd_bpe = 8;
                   1843:        else if ( max < 65536 )
                   1844:                nd_bpe = 16;
                   1845:        else
                   1846:                nd_bpe = 32;
                   1847:
                   1848:        nd_setup_parameters();
                   1849:        nd_free_private_storage();
                   1850:        len = 0;
                   1851:        for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
                   1852:                a = dptondv(mod,(DP)BDY(f)); ndv_removecont(mod,a); nd_ps[i] = a;
                   1853:                a = dptondv(0,(DP)BDY(f)); ndv_removecont(0,a); nd_psq[i] = a;
                   1854:                NEWRHist(r);
                   1855:                len = MAX(len,LEN(a));
                   1856:                SG(r) = HTD(a); TD(r) = HTD(a); ndl_copy(HDL(a),DL(r));
                   1857:                nd_psh[i] = r;
                   1858:        }
                   1859:        nmv_len = 16*len;
                   1860:        NEWNDV(ndv_red);
                   1861:        BDY(ndv_red) = (NMV)MALLOC(nmv_len*nmv_adv);
                   1862: }
                   1863:
1.1       noro     1864: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
                   1865: {
                   1866:        struct order_spec ord1;
                   1867:        VL fv,vv,vc;
                   1868:        NODE fd,fd0,r,r0,t,x,s,xx;
                   1869:        DP a,b,c;
                   1870:
                   1871:        get_vars((Obj)f,&fv); pltovl(v,&vv);
                   1872:        nd_nvar = length(vv);
                   1873:        if ( ord->id )
                   1874:                error("nd_gr : unsupported order");
                   1875:        switch ( ord->ord.simple ) {
                   1876:                case 0:
                   1877:                        is_rlex = 1;
                   1878:                        break;
                   1879:                case 1:
                   1880:                        is_rlex = 0;
                   1881:                        break;
                   1882:                default:
                   1883:                        error("nd_gr : unsupported order");
                   1884:        }
                   1885:        initd(ord);
                   1886:        for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                   1887:                ptod(CO,vv,(P)BDY(t),&b);
1.20      noro     1888:                NEXTNODE(fd0,fd); BDY(fd) = (pointer)b;
1.1       noro     1889:        }
                   1890:        if ( fd0 ) NEXT(fd) = 0;
1.23      noro     1891:        nd_setup(m,fd0);
                   1892:        x = nd_gb(m);
                   1893:        x = nd_reduceall(m,x);
                   1894:        for ( r0 = 0, t = x; t; t = NEXT(t) ) {
1.1       noro     1895:                NEXTNODE(r0,r);
1.20      noro     1896:                if ( m ) {
1.23      noro     1897:                        a = ndvtodp(m,BDY(t));
1.16      noro     1898:                        _dtop_mod(CO,vv,a,(P *)&BDY(r));
1.20      noro     1899:                } else {
1.23      noro     1900:                        a = ndvtodp(m,BDY(t));
1.16      noro     1901:                        dtop(CO,vv,a,(P *)&BDY(r));
1.20      noro     1902:                }
                   1903:        }
                   1904:        if ( r0 ) NEXT(r) = 0;
                   1905:        MKLIST(*rp,r0);
                   1906:        fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
                   1907:                nd_found,nd_notfirst,nd_create);
                   1908: }
                   1909:
1.23      noro     1910: void nd_gr_trace(LIST f,LIST v,int m,int homo,struct order_spec *ord,LIST *rp)
1.20      noro     1911: {
                   1912:        struct order_spec ord1;
                   1913:        VL fv,vv,vc;
1.23      noro     1914:        NODE fd,fd0,r,r0,t,x,s,xx,t0;
                   1915:        DP a,b,c,h;
                   1916:        NDV *w;
                   1917:        int len,i,j;
1.20      noro     1918:
                   1919:        get_vars((Obj)f,&fv); pltovl(v,&vv);
                   1920:        nd_nvar = length(vv);
                   1921:        if ( ord->id )
                   1922:                error("nd_gr : unsupported order");
                   1923:        initd(ord);
1.23      noro     1924:        if ( homo ) {
                   1925:                homogenize_order(ord,nd_nvar,&ord1);
                   1926:                switch ( ord1.ord.simple ) {
                   1927:                        case 0: is_rlex = 1; break;
                   1928:                        case 1: is_rlex = 0; break;
                   1929:                        default: error("nd_gr : unsupported order");
                   1930:                }
                   1931:                for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                   1932:                        ptod(CO,vv,(P)BDY(t),&c);
                   1933:                        if ( c ) {
                   1934:                                dp_homo(c,&h); NEXTNODE(fd0,fd); BDY(fd) = (pointer)h;
                   1935:                        }
                   1936:                }
                   1937:                if ( fd0 ) NEXT(fd) = 0;
                   1938:                initd(&ord1);
                   1939:                nd_nvar++;
                   1940:        } else {
                   1941:                switch ( ord->ord.simple ) {
                   1942:                        case 0: is_rlex = 1; break;
                   1943:                        case 1: is_rlex = 0; break;
                   1944:                        default: error("nd_gr : unsupported order");
                   1945:                }
                   1946:                for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                   1947:                        ptod(CO,vv,(P)BDY(t),&c);
                   1948:                        if ( c ) {
                   1949:                                NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
                   1950:                        }
1.20      noro     1951:                }
1.23      noro     1952:                if ( fd0 ) NEXT(fd) = 0;
1.20      noro     1953:        }
                   1954:        /* setup over GF(m) */
1.23      noro     1955:        nd_setup_trace(m,fd0);
                   1956:        x = nd_gb_trace(m);
                   1957:        if ( homo ) {
                   1958:                /* dehomogenization */
                   1959:                for ( t = x; t; t = NEXT(t) )
                   1960:                        ndv_dehomogenize((NDV)BDY(t));
                   1961:                nd_nvar--;
                   1962:                nd_setup_parameters();
                   1963:                initd(ord);
                   1964:                len = length(x);
                   1965:                w = (NDV *)ALLOCA(len*sizeof(NDV));
                   1966:                for ( i = 0, t = x; i < len; i++, t = NEXT(t) )
                   1967:                        w[i] = BDY(t);
                   1968:                for ( i = 0; i < len; i++ ) {
                   1969:                        for ( j = 0; j < i; j++ ) {
                   1970:                                if ( w[i] && w[j] )
                   1971:                                        if ( ndl_reducible(HDL(w[i]),HDL(w[j])) )
                   1972:                                                w[i] = 0;
                   1973:                                        else if ( ndl_reducible(HDL(w[j]),HDL(w[i])) )
                   1974:                                                w[j] = 0;
                   1975:                        }
                   1976:                }
                   1977:                for ( i = len-1, t0 = 0; i >= 0; i-- ) {
                   1978:                        if ( w[i] ) {
                   1979:                                NEXTNODE(t0,t);
                   1980:                                BDY(t) = (pointer)w[i];
                   1981:                        }
                   1982:                }
                   1983:                NEXT(t) = 0;
                   1984:                x = t0;
                   1985:        }
                   1986:        fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
                   1987:                nd_found,nd_notfirst,nd_create);
                   1988:        x = nd_reduceall(0,x);
1.20      noro     1989:        for ( r0 = 0; x; x = NEXT(x) ) {
                   1990:                NEXTNODE(r0,r);
1.23      noro     1991:                a = ndvtodp(0,(NDV)BDY(x));
1.20      noro     1992:                dtop(CO,vv,a,(P *)&BDY(r));
1.1       noro     1993:        }
                   1994:        if ( r0 ) NEXT(r) = 0;
                   1995:        MKLIST(*rp,r0);
                   1996: }
                   1997:
                   1998: void dltondl(int n,DL dl,unsigned int *r)
                   1999: {
                   2000:        unsigned int *d;
                   2001:        int i;
                   2002:
                   2003:        d = dl->d;
                   2004:        bzero(r,nd_wpd*sizeof(unsigned int));
                   2005:        if ( is_rlex )
                   2006:                for ( i = 0; i < n; i++ )
                   2007:                        r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe));
                   2008:        else
                   2009:                for ( i = 0; i < n; i++ )
                   2010:                        r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe);
                   2011: }
                   2012:
                   2013: DL ndltodl(int n,int td,unsigned int *ndl)
                   2014: {
                   2015:        DL dl;
                   2016:        int *d;
                   2017:        int i;
                   2018:
                   2019:        NEWDL(dl,n);
1.14      noro     2020:        TD(dl) = td;
1.1       noro     2021:        d = dl->d;
                   2022:        if ( is_rlex )
                   2023:                for ( i = 0; i < n; i++ )
                   2024:                        d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
                   2025:                                &((1<<nd_bpe)-1);
                   2026:        else
                   2027:                for ( i = 0; i < n; i++ )
                   2028:                        d[i] = (ndl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
                   2029:                                &((1<<nd_bpe)-1);
                   2030:        return dl;
                   2031: }
                   2032:
1.17      noro     2033: ND dptond(int mod,DP p)
1.1       noro     2034: {
                   2035:        ND d;
                   2036:        NM m0,m;
                   2037:        MP t;
                   2038:        int n;
                   2039:
                   2040:        if ( !p )
                   2041:                return 0;
                   2042:        n = NV(p);
                   2043:        m0 = 0;
                   2044:        for ( t = BDY(p); t; t = NEXT(t) ) {
                   2045:                NEXTNM(m0,m);
1.17      noro     2046:                if ( mod )
                   2047:                        CM(m) = ITOS(C(t));
                   2048:                else
                   2049:                        CQ(m) = (Q)C(t);
1.14      noro     2050:                TD(m) = TD(DL(t));
                   2051:                dltondl(n,DL(t),DL(m));
1.1       noro     2052:        }
                   2053:        NEXT(m) = 0;
                   2054:        MKND(n,m0,d);
1.14      noro     2055:        NV(d) = n;
                   2056:        SG(d) = SG(p);
1.1       noro     2057:        return d;
                   2058: }
                   2059:
1.17      noro     2060: DP ndtodp(int mod,ND p)
1.1       noro     2061: {
                   2062:        DP d;
                   2063:        MP m0,m;
                   2064:        NM t;
                   2065:        int n;
                   2066:
                   2067:        if ( !p )
                   2068:                return 0;
                   2069:        n = NV(p);
                   2070:        m0 = 0;
                   2071:        for ( t = BDY(p); t; t = NEXT(t) ) {
                   2072:                NEXTMP(m0,m);
1.17      noro     2073:                if ( mod )
                   2074:                        C(m) = STOI(CM(t));
                   2075:                else
                   2076:                        C(m) = (P)CQ(t);
1.14      noro     2077:                DL(m) = ndltodl(n,TD(t),DL(t));
1.1       noro     2078:        }
                   2079:        NEXT(m) = 0;
                   2080:        MKDP(n,m0,d);
1.14      noro     2081:        SG(d) = SG(p);
1.1       noro     2082:        return d;
                   2083: }
                   2084:
                   2085: void ndl_print(unsigned int *dl)
                   2086: {
                   2087:        int n;
                   2088:        int i;
                   2089:
                   2090:        n = nd_nvar;
                   2091:        printf("<<");
                   2092:        if ( is_rlex )
                   2093:                for ( i = 0; i < n; i++ )
                   2094:                        printf(i==n-1?"%d":"%d,",
                   2095:                                (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
                   2096:                                        &((1<<nd_bpe)-1));
                   2097:        else
                   2098:                for ( i = 0; i < n; i++ )
                   2099:                        printf(i==n-1?"%d":"%d,",
                   2100:                                (dl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
                   2101:                                        &((1<<nd_bpe)-1));
                   2102:        printf(">>");
                   2103: }
                   2104:
                   2105: void nd_print(ND p)
                   2106: {
                   2107:        NM m;
                   2108:
                   2109:        if ( !p )
                   2110:                printf("0\n");
                   2111:        else {
                   2112:                for ( m = BDY(p); m; m = NEXT(m) ) {
1.14      noro     2113:                        printf("+%d*",CM(m));
                   2114:                        ndl_print(DL(m));
1.1       noro     2115:                }
                   2116:                printf("\n");
                   2117:        }
                   2118: }
                   2119:
1.16      noro     2120: void nd_print_q(ND p)
                   2121: {
                   2122:        NM m;
                   2123:
                   2124:        if ( !p )
                   2125:                printf("0\n");
                   2126:        else {
                   2127:                for ( m = BDY(p); m; m = NEXT(m) ) {
                   2128:                        printf("+");
                   2129:                        printexpr(CO,CQ(m));
                   2130:                        printf("*");
                   2131:                        ndl_print(DL(m));
                   2132:                }
                   2133:                printf("\n");
                   2134:        }
                   2135: }
                   2136:
1.1       noro     2137: void ndp_print(ND_pairs d)
                   2138: {
                   2139:        ND_pairs t;
                   2140:
                   2141:        for ( t = d; t; t = NEXT(t) ) {
                   2142:                printf("%d,%d ",t->i1,t->i2);
                   2143:        }
                   2144:        printf("\n");
                   2145: }
                   2146:
1.20      noro     2147: void nd_removecont(int mod,ND p)
1.16      noro     2148: {
                   2149:        int i,n;
                   2150:        Q *w;
                   2151:        Q dvr,t;
                   2152:        NM m;
1.21      noro     2153:        struct oVECT v;
                   2154:        N q,r;
1.16      noro     2155:
1.20      noro     2156:        if ( mod )
                   2157:                nd_mul_c(mod,p,invm(HCM(p),mod));
                   2158:        else {
                   2159:                for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
                   2160:                w = (Q *)ALLOCA(n*sizeof(Q));
1.21      noro     2161:                v.len = n;
                   2162:                v.body = (pointer *)w;
1.20      noro     2163:                for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ )
                   2164:                        w[i] = CQ(m);
1.21      noro     2165:                removecont_array(w,n);
                   2166:                for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ ) CQ(m) = w[i];
1.16      noro     2167:        }
                   2168: }
                   2169:
1.21      noro     2170: void nd_removecont2(ND p1,ND p2)
                   2171: {
                   2172:        int i,n1,n2,n;
                   2173:        Q *w;
                   2174:        Q dvr,t;
                   2175:        NM m;
                   2176:        struct oVECT v;
                   2177:        N q,r;
                   2178:
                   2179:        if ( !p1 ) {
                   2180:                nd_removecont(0,p2); return;
                   2181:        } else if ( !p2 ) {
                   2182:                nd_removecont(0,p1); return;
                   2183:        }
                   2184:        n1 = nd_length(p1);
                   2185:        n2 = nd_length(p2);
                   2186:        n = n1+n2;
                   2187:        w = (Q *)ALLOCA(n*sizeof(Q));
                   2188:        v.len = n;
                   2189:        v.body = (pointer *)w;
                   2190:        for ( m = BDY(p1), i = 0; i < n1; m = NEXT(m), i++ )
                   2191:                w[i] = CQ(m);
                   2192:        for ( m = BDY(p2); i < n; m = NEXT(m), i++ )
                   2193:                w[i] = CQ(m);
                   2194:        removecont_array(w,n);
                   2195:        for ( m = BDY(p1), i = 0; i < n1; m = NEXT(m), i++ ) CQ(m) = w[i];
                   2196:        for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) CQ(m) = w[i];
                   2197: }
                   2198:
1.20      noro     2199: void ndv_removecont(int mod,NDV p)
1.16      noro     2200: {
                   2201:        int i,len;
                   2202:        Q *w;
                   2203:        Q dvr,t;
                   2204:        NMV m;
                   2205:
1.20      noro     2206:        if ( mod )
                   2207:                ndv_mul_c(mod,p,invm(HCM(p),mod));
                   2208:        else {
                   2209:                len = p->len;
                   2210:                w = (Q *)ALLOCA(len*sizeof(Q));
                   2211:                for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ )
                   2212:                        w[i] = CQ(m);
                   2213:                sortbynm(w,len);
                   2214:                qltozl(w,len,&dvr);
                   2215:                for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ ) {
                   2216:                        divq(CQ(m),dvr,&t); CQ(m) = t;
                   2217:                }
1.16      noro     2218:        }
1.21      noro     2219: }
                   2220:
1.23      noro     2221: void ndv_dehomogenize(NDV p)
                   2222: {
                   2223:        int i,len,newnvar,newwpd,newadv;
                   2224:        Q *w;
                   2225:        Q dvr,t;
                   2226:        NMV m,r;
                   2227:        unsigned int *d;
                   2228: #define NEWADV(m) (m = (NMV)(((char *)m)+newadv))
                   2229:
                   2230:        len = p->len;
                   2231:        newnvar = nd_nvar-1;
                   2232:        newwpd = newnvar/nd_epw+(newnvar%nd_epw?1:0);
                   2233:        for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ )
                   2234:                TD(m) -= ndl_dehomogenize(DL(m));
                   2235:        if ( newwpd != nd_wpd ) {
                   2236:                d = (unsigned int *)ALLOCA(newwpd*sizeof(unsigned int));
                   2237:                newadv = sizeof(struct oNMV)+(newwpd-1)*sizeof(unsigned int);
                   2238:                for ( m = r = BDY(p), i = 0; i < len; NMV_ADV(m), NEWADV(r), i++ ) {
                   2239:                        CQ(r) = CQ(m); TD(r) = TD(m); ndl_copy(DL(m),d); ndl_copy(d,DL(r));
                   2240:                }
                   2241:        }
                   2242:        NV(p)--;
                   2243: }
                   2244:
1.21      noro     2245: void removecont_array(Q *c,int n)
                   2246: {
                   2247:        struct oVECT v;
                   2248:        Q d0,d1,a,u,u1,gcd;
                   2249:        int i;
                   2250:        N qn,rn,gn;
                   2251:        Q *q,*r;
                   2252:
                   2253:        q = (Q *)ALLOCA(n*sizeof(Q));
                   2254:        r = (Q *)ALLOCA(n*sizeof(Q));
                   2255:        v.id = O_VECT; v.len = n; v.body = (pointer *)c;
                   2256:        igcdv_estimate(&v,&d0);
                   2257:        for ( i = 0; i < n; i++ ) {
                   2258:                divn(NM(c[i]),NM(d0),&qn,&rn);
                   2259:                NTOQ(qn,SGN(c[i])*SGN(d0),q[i]);
                   2260:                NTOQ(rn,SGN(c[i]),r[i]);
                   2261:        }
                   2262:        for ( i = 0; i < n; i++ )
                   2263:                if ( r[i] )
                   2264:                        break;
                   2265:        if ( i < n ) {
                   2266:                v.id = O_VECT; v.len = n; v.body = (pointer *)r;
                   2267:                igcdv(&v,&d1);
                   2268:                gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
                   2269:                divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
                   2270:                for ( i = 0; i < n; i++ ) {
                   2271:                        mulq(a,q[i],&u);
                   2272:                        if ( r[i] ) {
                   2273:                                divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
                   2274:                                addq(u,u1,&q[i]);
                   2275:                        } else
                   2276:                                q[i] = u;
                   2277:                }
                   2278:        }
                   2279:        for ( i = 0; i < n; i++ )
                   2280:                c[i] = q[i];
1.16      noro     2281: }
                   2282:
1.19      noro     2283: void nd_mul_c(int mod,ND p,int mul)
1.1       noro     2284: {
                   2285:        NM m;
                   2286:        int c,c1;
                   2287:
                   2288:        if ( !p )
                   2289:                return;
                   2290:        for ( m = BDY(p); m; m = NEXT(m) ) {
1.14      noro     2291:                c1 = CM(m);
1.19      noro     2292:                DMAR(c1,mul,0,mod,c);
1.14      noro     2293:                CM(m) = c;
1.1       noro     2294:        }
                   2295: }
                   2296:
1.16      noro     2297: void nd_mul_c_q(ND p,Q mul)
                   2298: {
                   2299:        NM m;
                   2300:        Q c;
                   2301:
                   2302:        if ( !p )
                   2303:                return;
                   2304:        for ( m = BDY(p); m; m = NEXT(m) ) {
                   2305:                mulq(CQ(m),mul,&c); CQ(m) = c;
                   2306:        }
                   2307: }
                   2308:
1.1       noro     2309: void nd_free(ND p)
                   2310: {
                   2311:        NM t,s;
                   2312:
                   2313:        if ( !p )
                   2314:                return;
                   2315:        t = BDY(p);
                   2316:        while ( t ) {
                   2317:                s = NEXT(t);
                   2318:                FREENM(t);
                   2319:                t = s;
                   2320:        }
                   2321:        FREEND(p);
                   2322: }
                   2323:
1.23      noro     2324: void ndv_free(NDV p)
                   2325: {
                   2326:        GC_free(BDY(p));
                   2327: }
                   2328:
1.1       noro     2329: void nd_append_red(unsigned int *d,int td,int i)
                   2330: {
1.13      noro     2331:        RHist m,m0;
1.1       noro     2332:        int h;
                   2333:
1.13      noro     2334:        NEWRHist(m);
1.1       noro     2335:        h = ndl_hash_value(td,d);
1.13      noro     2336:        m->index = i;
1.14      noro     2337:        TD(m) = td;
                   2338:        ndl_copy(d,DL(m));
1.1       noro     2339:        NEXT(m) = nd_red[h];
                   2340:        nd_red[h] = m;
                   2341: }
                   2342:
1.5       noro     2343: unsigned int *dp_compute_bound(DP p)
                   2344: {
                   2345:        unsigned int *d,*d1,*d2,*t;
                   2346:        MP m;
1.7       noro     2347:        int i,l;
1.5       noro     2348:
                   2349:        if ( !p )
                   2350:                return 0;
                   2351:        d1 = (unsigned int *)ALLOCA(nd_nvar*sizeof(unsigned int));
                   2352:        d2 = (unsigned int *)ALLOCA(nd_nvar*sizeof(unsigned int));
                   2353:        m = BDY(p);
1.14      noro     2354:        bcopy(DL(m)->d,d1,nd_nvar*sizeof(unsigned int));
1.5       noro     2355:        for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
1.14      noro     2356:                d = DL(m)->d;
1.5       noro     2357:                for ( i = 0; i < nd_nvar; i++ )
                   2358:                        d2[i] = d[i] > d1[i] ? d[i] : d1[i];
                   2359:                t = d1; d1 = d2; d2 = t;
                   2360:        }
1.13      noro     2361:        l = (nd_nvar+31);
1.7       noro     2362:        t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
                   2363:        bzero(t,l*sizeof(unsigned int));
1.5       noro     2364:        bcopy(d1,t,nd_nvar*sizeof(unsigned int));
                   2365:        return t;
                   2366: }
                   2367:
1.1       noro     2368: unsigned int *nd_compute_bound(ND p)
                   2369: {
                   2370:        unsigned int *d1,*d2,*t;
1.9       noro     2371:        int i,l;
1.1       noro     2372:        NM m;
                   2373:
                   2374:        if ( !p )
                   2375:                return 0;
                   2376:        d1 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                   2377:        d2 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                   2378:        bcopy(HDL(p),d1,nd_wpd*sizeof(unsigned int));
                   2379:        for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
1.14      noro     2380:                ndl_lcm(DL(m),d1,d2);
1.1       noro     2381:                t = d1; d1 = d2; d2 = t;
                   2382:        }
1.12      noro     2383:        l = nd_nvar+31;
1.9       noro     2384:        t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
                   2385:        bzero(t,l*sizeof(unsigned int));
1.5       noro     2386:        for ( i = 0; i < nd_nvar; i++ )
                   2387:                t[i] = (d1[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))&nd_mask0;
1.1       noro     2388:        return t;
                   2389: }
                   2390:
                   2391: void nd_setup_parameters() {
                   2392:        int i;
                   2393:
                   2394:        nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
                   2395:        nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
                   2396:        if ( nd_bpe < 32 ) {
                   2397:                nd_mask0 = (1<<nd_bpe)-1;
                   2398:        } else {
                   2399:                nd_mask0 = 0xffffffff;
                   2400:        }
                   2401:        bzero(nd_mask,sizeof(nd_mask));
                   2402:        nd_mask1 = 0;
                   2403:        for ( i = 0; i < nd_epw; i++ ) {
                   2404:                nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
                   2405:                nd_mask1 |= (1<<(nd_bpe-1))<<(i*nd_bpe);
                   2406:        }
1.13      noro     2407:        nm_adv = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);
1.3       noro     2408:        nmv_adv = sizeof(struct oNMV)+(nd_wpd-1)*sizeof(unsigned int);
1.1       noro     2409: }
                   2410:
1.20      noro     2411: /* mod < 0 => realloc nd_ps and pd_psq */
                   2412:
                   2413: ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d)
1.1       noro     2414: {
1.11      noro     2415:        int i,obpe,oadv;
1.13      noro     2416:        NM prev_nm_free_list;
                   2417:        RHist mr0,mr;
                   2418:        RHist r;
1.1       noro     2419:        ND_pairs s0,s,t,prev_ndp_free_list;
1.15      noro     2420:
1.1       noro     2421:        obpe = nd_bpe;
1.11      noro     2422:        oadv = nmv_adv;
1.5       noro     2423:        if ( obpe < 4 )
                   2424:                nd_bpe = 4;
                   2425:        else if ( obpe < 6 )
                   2426:                nd_bpe = 6;
                   2427:        else if ( obpe < 8 )
                   2428:                nd_bpe = 8;
                   2429:        else if ( obpe < 16 )
                   2430:                nd_bpe = 16;
                   2431:        else if ( obpe < 32 )
                   2432:                nd_bpe = 32;
                   2433:        else
                   2434:                error("nd_reconstruct : exponent too large");
                   2435:
1.1       noro     2436:        nd_setup_parameters();
                   2437:        prev_nm_free_list = _nm_free_list;
                   2438:        prev_ndp_free_list = _ndp_free_list;
                   2439:        _nm_free_list = 0;
                   2440:        _ndp_free_list = 0;
1.20      noro     2441:        if ( mod != 0 )
                   2442:                for ( i = nd_psn-1; i >= 0; i-- )
                   2443:                        ndv_realloc(nd_ps[i],obpe,oadv);
                   2444:        if ( !mod || trace )
                   2445:                for ( i = nd_psn-1; i >= 0; i-- )
                   2446:                        ndv_realloc(nd_psq[i],obpe,oadv);
1.1       noro     2447:        s0 = 0;
                   2448:        for ( t = d; t; t = NEXT(t) ) {
                   2449:                NEXTND_pairs(s0,s);
                   2450:                s->i1 = t->i1;
                   2451:                s->i2 = t->i2;
1.14      noro     2452:                TD(s) = TD(t);
                   2453:                SG(s) = SG(t);
1.1       noro     2454:                ndl_dup(obpe,t->lcm,s->lcm);
                   2455:        }
1.6       noro     2456:        for ( i = 0; i < REDTAB_LEN; i++ ) {
1.13      noro     2457:                for ( mr0 = 0, r = nd_red[i]; r; r = NEXT(r) ) {
1.16      noro     2458:                        NEXTRHist(mr0,mr);
1.13      noro     2459:                        mr->index = r->index;
1.20      noro     2460:                        SG(mr) = SG(r);
1.14      noro     2461:                        TD(mr) = TD(r);
                   2462:                        ndl_dup(obpe,DL(r),DL(mr));
1.6       noro     2463:                }
                   2464:                if ( mr0 ) NEXT(mr) = 0;
                   2465:                nd_red[i] = mr0;
                   2466:        }
1.11      noro     2467:        for ( i = 0; i < nd_psn; i++ ) {
1.20      noro     2468:                NEWRHist(r); SG(r) = SG(nd_psh[i]);
                   2469:                TD(r) = TD(nd_psh[i]); ndl_dup(obpe,DL(nd_psh[i]),DL(r));
1.13      noro     2470:                nd_psh[i] = r;
1.11      noro     2471:        }
1.1       noro     2472:        if ( s0 ) NEXT(s) = 0;
                   2473:        prev_nm_free_list = 0;
                   2474:        prev_ndp_free_list = 0;
1.3       noro     2475:        BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
1.1       noro     2476:        GC_gcollect();
                   2477:        return s0;
                   2478: }
                   2479:
1.23      noro     2480: void nd_reconstruct_direct(int mod,NDV *ps,int len)
                   2481: {
                   2482:        int i,obpe,oadv;
                   2483:        NM prev_nm_free_list;
                   2484:        RHist mr0,mr;
                   2485:        RHist r;
                   2486:        ND_pairs s0,s,t,prev_ndp_free_list;
                   2487:
                   2488:        obpe = nd_bpe;
                   2489:        oadv = nmv_adv;
                   2490:        if ( obpe < 4 )
                   2491:                nd_bpe = 4;
                   2492:        else if ( obpe < 6 )
                   2493:                nd_bpe = 6;
                   2494:        else if ( obpe < 8 )
                   2495:                nd_bpe = 8;
                   2496:        else if ( obpe < 16 )
                   2497:                nd_bpe = 16;
                   2498:        else if ( obpe < 32 )
                   2499:                nd_bpe = 32;
                   2500:        else
                   2501:                error("nd_reconstruct_direct : exponent too large");
                   2502:
                   2503:        nd_setup_parameters();
                   2504:        prev_nm_free_list = _nm_free_list;
                   2505:        prev_ndp_free_list = _ndp_free_list;
                   2506:        _nm_free_list = 0;
                   2507:        _ndp_free_list = 0;
                   2508:        if ( mod != 0 )
                   2509:                for ( i = len-1; i >= 0; i-- )
                   2510:                        ndv_realloc(ps[i],obpe,oadv);
                   2511:        prev_nm_free_list = 0;
                   2512:        prev_ndp_free_list = 0;
                   2513:        BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
                   2514:        GC_gcollect();
                   2515: }
                   2516:
1.1       noro     2517: void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
                   2518: {
                   2519:        int n,i,ei,oepw,cepw,cbpe;
                   2520:
                   2521:        n = nd_nvar;
                   2522:        oepw = (sizeof(unsigned int)*8)/obpe;
                   2523:        cepw = nd_epw;
                   2524:        cbpe = nd_bpe;
1.15      noro     2525:        for ( i = 0; i < nd_wpd; i++ )
                   2526:                r[i] = 0;
1.1       noro     2527:        if ( is_rlex )
                   2528:                for ( i = 0; i < n; i++ ) {
                   2529:                        ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe))
                   2530:                                &((1<<obpe)-1);
                   2531:                        r[(n-1-i)/cepw] |= (ei<<((cepw-((n-1-i)%cepw)-1)*cbpe));
                   2532:                }
                   2533:        else
                   2534:                for ( i = 0; i < n; i++ ) {
                   2535:                        ei = (d[i/oepw]>>((oepw-(i%oepw)-1)*obpe))
                   2536:                                &((1<<obpe)-1);
                   2537:                        r[i/cepw] |= (ei<<((cepw-(i%cepw)-1)*cbpe));
                   2538:                }
                   2539: }
                   2540:
1.11      noro     2541: void nd_realloc(ND p,int obpe)
1.1       noro     2542: {
                   2543:        NM m,mr,mr0;
                   2544:
1.11      noro     2545:        if ( p ) {
                   2546:                m = BDY(p);
1.1       noro     2547:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   2548:                        NEXTNM(mr0,mr);
1.14      noro     2549:                        CM(mr) = CM(m);
                   2550:                        TD(mr) = TD(m);
                   2551:                        ndl_dup(obpe,DL(m),DL(mr));
1.1       noro     2552:                }
                   2553:                NEXT(mr) = 0;
1.11      noro     2554:                BDY(p) = mr0;
1.1       noro     2555:        }
                   2556: }
1.3       noro     2557:
1.6       noro     2558: ND nd_copy(ND p)
                   2559: {
                   2560:        NM m,mr,mr0;
                   2561:        int c,n,s;
                   2562:        ND r;
                   2563:
                   2564:        if ( !p )
                   2565:                return 0;
                   2566:        else {
                   2567:                s = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);
                   2568:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                   2569:                        NEXTNM(mr0,mr);
1.14      noro     2570:                        CM(mr) = CM(m);
                   2571:                        TD(mr) = TD(m);
                   2572:                        ndl_copy(DL(m),DL(mr));
1.6       noro     2573:                }
                   2574:                NEXT(mr) = 0;
                   2575:                MKND(NV(p),mr0,r);
1.14      noro     2576:                SG(r) = SG(p);
1.6       noro     2577:                return r;
                   2578:        }
                   2579: }
                   2580:
1.16      noro     2581: int nd_sp(int mod,ND_pairs p,ND *rp)
1.11      noro     2582: {
                   2583:        NM m;
                   2584:        NDV p1,p2;
                   2585:        ND t1,t2;
                   2586:        unsigned int *lcm;
                   2587:        int td;
                   2588:
1.20      noro     2589:        if ( mod ) {
                   2590:                p1 = nd_ps[p->i1]; p2 = nd_ps[p->i2];
                   2591:        } else {
                   2592:                p1 = nd_psq[p->i1]; p2 = nd_psq[p->i2];
                   2593:        }
1.11      noro     2594:        lcm = p->lcm;
1.14      noro     2595:        td = TD(p);
1.11      noro     2596:        NEWNM(m);
1.20      noro     2597:        CQ(m) = HCQ(p2);
1.16      noro     2598:        TD(m) = td-HTD(p1); ndl_sub(lcm,HDL(p1),DL(m));
1.14      noro     2599:        if ( ndl_check_bound2(p->i1,DL(m)) )
1.11      noro     2600:                return 0;
1.16      noro     2601:        t1 = ndv_mul_nm_create(mod,p1,m);
                   2602:        if ( mod )
                   2603:                CM(m) = mod-HCM(p1);
                   2604:        else
                   2605:                chsgnq(HCQ(p1),&CQ(m));
                   2606:        TD(m) = td-HTD(p2); ndl_sub(lcm,HDL(p2),DL(m));
1.14      noro     2607:        if ( ndl_check_bound2(p->i2,DL(m)) ) {
1.11      noro     2608:                nd_free(t1);
                   2609:                return 0;
                   2610:        }
1.16      noro     2611:        ndv_mul_nm(mod,p2,m,ndv_red);
1.11      noro     2612:        FREENM(m);
1.20      noro     2613:        *rp = ndv_add(mod,t1,ndv_red);
1.11      noro     2614:        return 1;
                   2615: }
                   2616:
1.19      noro     2617: void ndv_mul_c(int mod,NDV p,int mul)
1.11      noro     2618: {
                   2619:        NMV m;
                   2620:        int c,c1,len,i;
                   2621:
                   2622:        if ( !p )
                   2623:                return;
1.14      noro     2624:        len = LEN(p);
1.11      noro     2625:        for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
1.14      noro     2626:                c1 = CM(m);
1.19      noro     2627:                DMAR(c1,mul,0,mod,c);
1.14      noro     2628:                CM(m) = c;
1.11      noro     2629:        }
                   2630: }
                   2631:
1.16      noro     2632: void ndv_mul_c_q(NDV p,Q mul)
                   2633: {
                   2634:        NMV m;
                   2635:        Q c;
                   2636:        int len,i;
                   2637:
                   2638:        if ( !p )
                   2639:                return;
                   2640:        len = LEN(p);
                   2641:        for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
                   2642:                mulq(CQ(m),mul,&c); CQ(m) = c;
                   2643:        }
                   2644: }
                   2645:
                   2646: void ndv_mul_nm(int mod,NDV p,NM m0,NDV r)
1.4       noro     2647: {
                   2648:        NMV m,mr,mr0;
                   2649:        unsigned int *d,*dt,*dm;
                   2650:        int c,n,td,i,c1,c2,len;
1.16      noro     2651:        Q q;
1.4       noro     2652:
                   2653:        if ( !p )
                   2654:                /* XXX */
1.14      noro     2655:                LEN(r) = 0;
1.4       noro     2656:        else {
1.14      noro     2657:                n = NV(p); m = BDY(p); len = LEN(p);
1.16      noro     2658:                d = DL(m0); td = TD(m0);
1.4       noro     2659:                mr = BDY(r);
1.16      noro     2660:                if ( mod ) {
                   2661:                        c = CM(m0);
                   2662:                        for ( ; len > 0; len--, NMV_ADV(m), NMV_ADV(mr) ) {
1.19      noro     2663:                                c1 = CM(m); DMAR(c1,c,0,mod,c2); CM(mr) = c2;
1.16      noro     2664:                                TD(mr) = TD(m)+td; ndl_add(DL(m),d,DL(mr));
                   2665:                        }
                   2666:                } else {
                   2667:                        q = CQ(m0);
                   2668:                        for ( ; len > 0; len--, NMV_ADV(m), NMV_ADV(mr) ) {
                   2669:                                mulq(CQ(m),q,&CQ(mr));
                   2670:                                TD(mr) = TD(m)+td; ndl_add(DL(m),d,DL(mr));
                   2671:                        }
1.9       noro     2672:                }
                   2673:                NV(r) = NV(p);
1.14      noro     2674:                LEN(r) = LEN(p);
                   2675:                SG(r) = SG(p) + td;
1.9       noro     2676:        }
                   2677: }
                   2678:
1.16      noro     2679: ND ndv_mul_nm_create(int mod,NDV p,NM m0)
1.9       noro     2680: {
                   2681:        NM mr,mr0;
                   2682:        NMV m;
                   2683:        unsigned int *d,*dt,*dm;
                   2684:        int c,n,td,i,c1,c2,len;
1.16      noro     2685:        Q q;
1.9       noro     2686:        ND r;
                   2687:
                   2688:        if ( !p )
                   2689:                return 0;
                   2690:        else {
                   2691:                n = NV(p); m = BDY(p);
1.16      noro     2692:                d = DL(m0); td = TD(m0);
1.14      noro     2693:                len = LEN(p);
1.9       noro     2694:                mr0 = 0;
1.16      noro     2695:                if ( mod ) {
                   2696:                        c = CM(m0);
                   2697:                        for ( i = 0; i < len; i++, NMV_ADV(m) ) {
                   2698:                                NEXTNM(mr0,mr);
                   2699:                                c1 = CM(m);
1.19      noro     2700:                                DMAR(c1,c,0,mod,c2);
1.16      noro     2701:                                CM(mr) = c2;
                   2702:                                TD(mr) = TD(m)+td;
                   2703:                                ndl_add(DL(m),d,DL(mr));
                   2704:                        }
                   2705:                } else {
                   2706:                        q = CQ(m0);
                   2707:                        for ( i = 0; i < len; i++, NMV_ADV(m) ) {
                   2708:                                NEXTNM(mr0,mr);
                   2709:                                mulq(CQ(m),q,&CQ(mr));
                   2710:                                TD(mr) = TD(m)+td;
                   2711:                                ndl_add(DL(m),d,DL(mr));
                   2712:                        }
1.4       noro     2713:                }
1.9       noro     2714:                NEXT(mr) = 0;
                   2715:                MKND(NV(p),mr0,r);
1.14      noro     2716:                SG(r) = SG(p) + td;
1.9       noro     2717:                return r;
1.4       noro     2718:        }
                   2719: }
                   2720:
1.19      noro     2721: ND ndv_add(int mod,ND p1,NDV p2)
1.4       noro     2722: {
1.9       noro     2723:        register NM prev,cur,new;
1.4       noro     2724:        int c,c1,c2,t,td,td2,mul,len,i;
1.9       noro     2725:        NM head;
1.4       noro     2726:        unsigned int *d;
                   2727:        NMV m2;
1.16      noro     2728:        Q q;
1.4       noro     2729:
                   2730:        if ( !p1 )
                   2731:                return 0;
1.20      noro     2732:        else if ( !mod )
                   2733:                return ndv_add_q(p1,p2);
1.4       noro     2734:        else {
                   2735:                prev = 0; head = cur = BDY(p1);
1.14      noro     2736:                NEWNM(new); len = LEN(p2);
1.9       noro     2737:                for ( m2 = BDY(p2), i = 0; cur && i < len; ) {
1.14      noro     2738:                        td2 = TD(new) = TD(m2);
                   2739:                        if ( TD(cur) > td2 ) {
1.13      noro     2740:                                prev = cur; cur = NEXT(cur);
                   2741:                                continue;
1.14      noro     2742:                        } else if ( TD(cur) < td2 ) c = -1;
1.13      noro     2743:                        else if ( nd_wpd == 1 ) {
1.14      noro     2744:                                if ( DL(cur)[0] > DL(m2)[0] ) c = is_rlex ? -1 : 1;
                   2745:                                else if ( DL(cur)[0] < DL(m2)[0] ) c = is_rlex ? 1 : -1;
1.13      noro     2746:                                else c = 0;
                   2747:                        }
1.14      noro     2748:                        else c = ndl_compare(DL(cur),DL(m2));
1.4       noro     2749:                        switch ( c ) {
                   2750:                                case 0:
1.19      noro     2751:                                        t = CM(m2)+CM(cur)-mod;
                   2752:                                        if ( t < 0 ) t += mod;
1.14      noro     2753:                                        if ( t ) CM(cur) = t;
1.4       noro     2754:                                        else if ( !prev ) {
1.9       noro     2755:                                                head = NEXT(cur); FREENM(cur); cur = head;
1.4       noro     2756:                                        } else {
1.9       noro     2757:                                                NEXT(prev) = NEXT(cur); FREENM(cur); cur = NEXT(prev);
1.4       noro     2758:                                        }
                   2759:                                        NMV_ADV(m2); i++;
                   2760:                                        break;
                   2761:                                case 1:
1.9       noro     2762:                                        prev = cur; cur = NEXT(cur);
1.4       noro     2763:                                        break;
                   2764:                                case -1:
1.14      noro     2765:                                        ndl_copy(DL(m2),DL(new));
1.16      noro     2766:                                        CQ(new) = CQ(m2);
                   2767:                                        if ( !prev ) {
                   2768:                                                /* cur = head */
                   2769:                                                prev = new; NEXT(prev) = head; head = prev;
                   2770:                                        } else {
                   2771:                                                NEXT(prev) = new; NEXT(new) = cur; prev = new;
                   2772:                                        }
                   2773:                                        NEWNM(new); NMV_ADV(m2); i++;
                   2774:                                        break;
                   2775:                        }
                   2776:                }
                   2777:                for ( ; i < len; i++, NMV_ADV(m2) ) {
                   2778:                        td2 = TD(new) = TD(m2); CQ(new) = CQ(m2); ndl_copy(DL(m2),DL(new));
                   2779:                        if ( !prev ) {
                   2780:                                prev = new; NEXT(prev) = 0; head = prev;
                   2781:                        } else {
                   2782:                                NEXT(prev) = new; NEXT(new) = 0; prev = new;
                   2783:                        }
                   2784:                        NEWNM(new);
                   2785:                }
                   2786:                FREENM(new);
                   2787:                if ( head ) {
                   2788:                        BDY(p1) = head; SG(p1) = MAX(SG(p1),SG(p2));
                   2789:                        return p1;
                   2790:                } else {
                   2791:                        FREEND(p1);
                   2792:                        return 0;
                   2793:                }
                   2794:
                   2795:        }
                   2796: }
                   2797:
                   2798: ND ndv_add_q(ND p1,NDV p2)
                   2799: {
                   2800:        register NM prev,cur,new;
                   2801:        int c,c1,c2,t,td,td2,mul,len,i;
                   2802:        NM head;
                   2803:        unsigned int *d;
                   2804:        NMV m2;
                   2805:        Q q;
                   2806:
                   2807:        if ( !p1 )
                   2808:                return 0;
                   2809:        else {
                   2810:                prev = 0; head = cur = BDY(p1);
                   2811:                NEWNM(new); len = LEN(p2);
                   2812:                for ( m2 = BDY(p2), i = 0; cur && i < len; ) {
                   2813:                        td2 = TD(new) = TD(m2);
                   2814:                        if ( TD(cur) > td2 ) {
                   2815:                                prev = cur; cur = NEXT(cur);
                   2816:                                continue;
                   2817:                        } else if ( TD(cur) < td2 ) c = -1;
                   2818:                        else if ( nd_wpd == 1 ) {
                   2819:                                if ( DL(cur)[0] > DL(m2)[0] ) c = is_rlex ? -1 : 1;
                   2820:                                else if ( DL(cur)[0] < DL(m2)[0] ) c = is_rlex ? 1 : -1;
                   2821:                                else c = 0;
                   2822:                        }
                   2823:                        else c = ndl_compare(DL(cur),DL(m2));
                   2824:                        switch ( c ) {
                   2825:                                case 0:
                   2826:                                        addq(CQ(cur),CQ(m2),&q);
                   2827:                                        if ( q )
                   2828:                                                CQ(cur) = q;
                   2829:                                        else if ( !prev ) {
                   2830:                                                head = NEXT(cur); FREENM(cur); cur = head;
                   2831:                                        } else {
                   2832:                                                NEXT(prev) = NEXT(cur); FREENM(cur); cur = NEXT(prev);
                   2833:                                        }
                   2834:                                        NMV_ADV(m2); i++;
                   2835:                                        break;
                   2836:                                case 1:
                   2837:                                        prev = cur; cur = NEXT(cur);
                   2838:                                        break;
                   2839:                                case -1:
                   2840:                                        ndl_copy(DL(m2),DL(new));
                   2841:                                        CQ(new) = CQ(m2);
1.4       noro     2842:                                        if ( !prev ) {
                   2843:                                                /* cur = head */
1.9       noro     2844:                                                prev = new; NEXT(prev) = head; head = prev;
1.4       noro     2845:                                        } else {
1.9       noro     2846:                                                NEXT(prev) = new; NEXT(new) = cur; prev = new;
1.4       noro     2847:                                        }
1.9       noro     2848:                                        NEWNM(new); NMV_ADV(m2); i++;
1.4       noro     2849:                                        break;
                   2850:                        }
                   2851:                }
1.9       noro     2852:                for ( ; i < len; i++, NMV_ADV(m2) ) {
1.16      noro     2853:                        td2 = TD(new) = TD(m2); CQ(new) = CQ(m2); ndl_copy(DL(m2),DL(new));
1.9       noro     2854:                        if ( !prev ) {
                   2855:                                prev = new; NEXT(prev) = 0; head = prev;
                   2856:                        } else {
                   2857:                                NEXT(prev) = new; NEXT(new) = 0; prev = new;
                   2858:                        }
                   2859:                        NEWNM(new);
                   2860:                }
1.4       noro     2861:                FREENM(new);
                   2862:                if ( head ) {
1.14      noro     2863:                        BDY(p1) = head; SG(p1) = MAX(SG(p1),SG(p2));
1.4       noro     2864:                        return p1;
                   2865:                } else {
                   2866:                        FREEND(p1);
                   2867:                        return 0;
                   2868:                }
                   2869:
                   2870:        }
                   2871: }
                   2872:
1.11      noro     2873: void ndv_realloc(NDV p,int obpe,int oadv)
                   2874: {
1.13      noro     2875:        NMV m,mr,mr0,t;
                   2876:        int len,i,k;
1.11      noro     2877:
1.13      noro     2878: #define NMV_OPREV(m) (m = (NMV)(((char *)m)-oadv))
                   2879: #define NMV_PREV(m) (m = (NMV)(((char *)m)-nmv_adv))
1.11      noro     2880:
                   2881:        if ( p ) {
1.14      noro     2882:                m = BDY(p); len = LEN(p);
1.15      noro     2883:                if ( nmv_adv > oadv )
                   2884:                        mr0 = (NMV)REALLOC(BDY(p),len*nmv_adv);
                   2885:                else
                   2886:                        mr0 = BDY(p);
1.13      noro     2887:                m = (NMV)((char *)mr0+(len-1)*oadv);
                   2888:                mr = (NMV)((char *)mr0+(len-1)*nmv_adv);
                   2889:                t = (NMV)ALLOCA(nmv_adv);
                   2890:                for ( i = 0; i < len; i++, NMV_OPREV(m), NMV_PREV(mr) ) {
1.16      noro     2891:                        CQ(t) = CQ(m);
1.14      noro     2892:                        TD(t) = TD(m);
                   2893:                        for ( k = 0; k < nd_wpd; k++ ) DL(t)[k] = 0;
                   2894:                        ndl_dup(obpe,DL(m),DL(t));
1.16      noro     2895:                        CQ(mr) = CQ(t);
1.14      noro     2896:                        TD(mr) = TD(t);
                   2897:                        ndl_copy(DL(t),DL(mr));
1.11      noro     2898:                }
                   2899:                BDY(p) = mr0;
                   2900:        }
                   2901: }
                   2902:
1.16      noro     2903: NDV ndtondv(int mod,ND p)
1.3       noro     2904: {
                   2905:        NDV d;
                   2906:        NMV m,m0;
                   2907:        NM t;
                   2908:        int i,len;
                   2909:
                   2910:        if ( !p )
                   2911:                return 0;
                   2912:        len = nd_length(p);
1.16      noro     2913:        if ( mod )
                   2914:                m0 = m = (NMV)MALLOC_ATOMIC(len*nmv_adv);
                   2915:        else
                   2916:                m0 = m = (NMV)MALLOC(len*nmv_adv);
1.3       noro     2917:        for ( t = BDY(p), i = 0; t; t = NEXT(t), i++, NMV_ADV(m) ) {
1.14      noro     2918:                TD(m) = TD(t);
                   2919:                ndl_copy(DL(t),DL(m));
1.16      noro     2920:                CQ(m) = CQ(t);
1.3       noro     2921:        }
                   2922:        MKNDV(NV(p),m0,len,d);
1.23      noro     2923:        SG(d) = SG(p);
                   2924:        return d;
                   2925: }
                   2926:
                   2927: ND ndvtond(int mod,NDV p)
                   2928: {
                   2929:        ND d;
                   2930:        NM m,m0;
                   2931:        NMV t;
                   2932:        int i,len;
                   2933:
                   2934:        if ( !p )
                   2935:                return 0;
                   2936:        m0 = 0;
                   2937:        len = p->len;
                   2938:        for ( t = BDY(p), i = 0; i < len; NMV_ADV(t), i++ ) {
                   2939:                NEXTNM(m0,m);
                   2940:                TD(m) = TD(t);
                   2941:                ndl_copy(DL(t),DL(m));
                   2942:                CQ(m) = CQ(t);
                   2943:        }
                   2944:        NEXT(m) = 0;
                   2945:        MKND(NV(p),m0,d);
1.14      noro     2946:        SG(d) = SG(p);
1.3       noro     2947:        return d;
                   2948: }
                   2949:
1.16      noro     2950: NDV dptondv(int mod,DP p)
1.11      noro     2951: {
                   2952:        NDV d;
                   2953:        NMV m0,m;
                   2954:        MP t;
1.20      noro     2955:        DP q;
1.11      noro     2956:        int l,i,n;
                   2957:
                   2958:        if ( !p )
                   2959:                return 0;
                   2960:        for ( t = BDY(p), l = 0; t; t = NEXT(t), l++ );
1.20      noro     2961:        if ( mod ) {
                   2962:                _dp_mod(p,mod,0,&q); p = q;
1.16      noro     2963:                m0 = m = (NMV)MALLOC_ATOMIC(l*nmv_adv);
1.20      noro     2964:        } else
1.16      noro     2965:                m0 = m = (NMV)MALLOC(l*nmv_adv);
1.11      noro     2966:        n = NV(p);
1.17      noro     2967:        for ( t = BDY(p), i = 0; i < l; i++, t = NEXT(t), NMV_ADV(m) ) {
                   2968:                if ( mod )
1.16      noro     2969:                        CM(m) = ITOS(C(t));
1.17      noro     2970:                else
1.16      noro     2971:                        CQ(m) = (Q)C(t);
1.17      noro     2972:                TD(m) = TD(DL(t));
                   2973:                dltondl(n,DL(t),DL(m));
1.11      noro     2974:        }
                   2975:        MKNDV(n,m0,l,d);
1.14      noro     2976:        SG(d) = SG(p);
1.11      noro     2977:        return d;
                   2978: }
                   2979:
1.16      noro     2980: DP ndvtodp(int mod,NDV p)
1.11      noro     2981: {
                   2982:        DP d;
                   2983:        MP m0,m;
                   2984:        NMV t;
                   2985:        int len,i,n;
                   2986:
                   2987:        if ( !p )
                   2988:                return 0;
                   2989:        m0 = 0;
1.14      noro     2990:        len = LEN(p);
1.11      noro     2991:        n = NV(p);
1.17      noro     2992:        for ( t = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) {
                   2993:                NEXTMP(m0,m);
                   2994:                if ( mod )
1.16      noro     2995:                        C(m) = STOI(CM(t));
1.17      noro     2996:                else
1.16      noro     2997:                        C(m) = (P)CQ(t);
1.17      noro     2998:                DL(m) = ndltodl(n,TD(t),DL(t));
1.11      noro     2999:        }
                   3000:        NEXT(m) = 0;
                   3001:        MKDP(NV(p),m0,d);
1.14      noro     3002:        SG(d) = SG(p);
1.11      noro     3003:        return d;
                   3004: }
                   3005:
1.3       noro     3006: void ndv_print(NDV p)
                   3007: {
                   3008:        NMV m;
                   3009:        int i,len;
                   3010:
                   3011:        if ( !p )
                   3012:                printf("0\n");
                   3013:        else {
1.14      noro     3014:                len = LEN(p);
1.3       noro     3015:                for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
1.14      noro     3016:                        printf("+%d*",CM(m));
1.16      noro     3017:                        ndl_print(DL(m));
                   3018:                }
                   3019:                printf("\n");
                   3020:        }
                   3021: }
                   3022:
                   3023: void ndv_print_q(NDV p)
                   3024: {
                   3025:        NMV m;
                   3026:        int i,len;
                   3027:
                   3028:        if ( !p )
                   3029:                printf("0\n");
                   3030:        else {
                   3031:                len = LEN(p);
                   3032:                for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
                   3033:                        printf("+");
                   3034:                        printexpr(CO,CQ(m));
                   3035:                        printf("*");
1.14      noro     3036:                        ndl_print(DL(m));
1.3       noro     3037:                }
                   3038:                printf("\n");
                   3039:        }
1.11      noro     3040: }

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