[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.20

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

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