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

1.3     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.2 2003/07/23 01:23:34 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;
                     36:        int td;
                     37:        int c;
                     38:        unsigned int dl[1];
                     39: } *NM;
                     40:
1.3     ! noro       41: typedef struct oNMV {
        !            42:        int td;
        !            43:        int c;
        !            44:        unsigned int dl[1];
        !            45: } *NMV;
        !            46:
1.1       noro       47: typedef struct oND_pairs {
                     48:        struct oND_pairs *next;
                     49:        int i1,i2;
                     50:        int td,sugar;
                     51:        unsigned int lcm[1];
                     52: } *ND_pairs;
                     53:
                     54: static ND *nd_ps;
1.3     ! noro       55: static NDV *nd_psv;
1.1       noro       56: static unsigned int **nd_bound;
                     57: int nd_mod,nd_nvar;
                     58: int is_rlex;
                     59: int nd_epw,nd_bpe,nd_wpd;
                     60: unsigned int nd_mask[32];
                     61: unsigned int nd_mask0,nd_mask1;
                     62: NM _nm_free_list;
                     63: ND _nd_free_list;
                     64: ND_pairs _ndp_free_list;
                     65: NM *nd_red;
                     66: int nd_red_len;
1.3     ! noro       67: int nmv_adv;
        !            68: int nmv_len;
        !            69: NDV ndv_red;
1.1       noro       70:
                     71: extern int Top,Reverse;
                     72: int nd_psn,nd_pslen;
                     73: int nd_found,nd_create,nd_notfirst;
                     74:
                     75: void GC_gcollect();
                     76: NODE append_one(NODE,int);
                     77:
                     78: #define HTD(d) ((d)->body->td)
                     79: #define HDL(d) ((d)->body->dl)
                     80: #define HC(d) ((d)->body->c)
1.3     ! noro       81: #define NMV_ADV(m) (m = (NMV)(((char *)m)+nmv_adv))
1.1       noro       82:
                     83: #define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list)
                     84: #define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list)
                     85: #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)
                     86:
1.3     ! noro       87: #define NEWNDV(d) ((d)=(NDV)MALLOC(sizeof(struct oNDV)))
        !            88: #define MKNDV(n,m,l,d) NEWNDV(d); (d)->nv=(n); BDY(d)=(m); (d)->len = l;
        !            89:
1.1       noro       90: #define NEXTNM(r,c) \
                     91: if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);}
                     92: #define NEXTNM2(r,c,s) \
                     93: if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
                     94: #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)
                     95: #define FREENDP(m) NEXT(m)=_ndp_free_list; _ndp_free_list=(m)
                     96: #define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m)
                     97:
                     98: #define NEXTND_pairs(r,c) \
                     99: if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
                    100:
                    101: ND_pairs crit_B( ND_pairs d, int s );
                    102: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp);
                    103: NODE nd_setup(NODE f);
                    104: int nd_newps(ND a);
                    105: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest );
                    106: NODE update_base(NODE nd,int ndp);
                    107: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest );
                    108: int crit_2( int dp1, int dp2 );
                    109: ND_pairs crit_F( ND_pairs d1 );
                    110: ND_pairs crit_M( ND_pairs d1 );
                    111: ND_pairs nd_newpairs( NODE g, int t );
                    112: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t);
                    113: NODE nd_gb(NODE f);
                    114: void nd_free_private_storage();
                    115: void _NM_alloc();
                    116: void _ND_alloc();
                    117: int ndl_td(unsigned int *d);
                    118: ND nd_add(ND p1,ND p2);
1.3     ! noro      119: ND ndv_add(ND p1,NDV p2);
1.1       noro      120: ND nd_mul_nm(ND p,NM m0);
                    121: ND nd_mul_term(ND p,int td,unsigned int *d);
                    122: int nd_sp(ND_pairs p,ND *nf);
                    123: int nd_find_reducer(ND g,ND *red);
1.3     ! noro      124: int ndv_find_reducer(ND g,NDV red);
1.1       noro      125: int nd_nf(ND g,int full,ND *nf);
                    126: ND nd_reduce(ND p1,ND p2);
                    127: ND nd_reduce_special(ND p1,ND p2);
                    128: void nd_free(ND p);
                    129: void ndl_print(unsigned int *dl);
                    130: void nd_print(ND p);
                    131: void ndp_print(ND_pairs d);
                    132: int nd_length(ND p);
                    133: void nd_monic(ND p);
                    134: void nd_mul_c(ND p,int mul);
                    135: void nd_free_redlist();
                    136: void nd_append_red(unsigned int *d,int td,int i);
                    137: unsigned int *nd_compute_bound(ND p);
                    138: ND_pairs nd_reconstruct(ND_pairs);
                    139: void nd_setup_parameters();
                    140: ND nd_dup(ND p,int obpe);
                    141: void ndl_dup(int obpe,unsigned int *d,unsigned int *r);
1.3     ! noro      142: NDV ndtondv(ND p);
        !           143: void ndv_mul_nm(NDV pv,NM m,NDV r);
1.1       noro      144:
                    145: void nd_free_private_storage()
                    146: {
                    147:        _nd_free_list = 0;
                    148:        _nm_free_list = 0;
                    149:        nd_red = 0;
                    150:        GC_gcollect();
                    151: }
                    152:
                    153: void _NM_alloc()
                    154: {
                    155:        NM p;
                    156:        int i;
                    157:
                    158:        for ( i = 0; i < 16; i++ ) {
                    159:                p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
                    160:                p->next = _nm_free_list; _nm_free_list = p;
                    161:        }
                    162: }
                    163:
                    164: void _ND_alloc()
                    165: {
                    166:        ND p;
                    167:        int i;
                    168:
                    169:        for ( i = 0; i < 1024; i++ ) {
                    170:                p = (ND)GC_malloc(sizeof(struct oND));
                    171:                p->body = (NM)_nd_free_list; _nd_free_list = p;
                    172:        }
                    173: }
                    174:
                    175: void _NDP_alloc()
                    176: {
                    177:        ND_pairs p;
                    178:        int i;
                    179:
                    180:        for ( i = 0; i < 10240; i++ ) {
                    181:                p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs)
                    182:                        +(nd_wpd-1)*sizeof(unsigned int));
                    183:                p->next = _ndp_free_list; _ndp_free_list = p;
                    184:        }
                    185: }
                    186:
                    187: INLINE nd_length(ND p)
                    188: {
                    189:        NM m;
                    190:        int i;
                    191:
                    192:        if ( !p )
                    193:                return 0;
                    194:        else {
                    195:                for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
                    196:                return i;
                    197:        }
                    198: }
                    199:
                    200: INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2)
                    201: {
                    202:        unsigned int u1,u2;
                    203:        int i,j;
                    204:
                    205:        switch ( nd_bpe ) {
                    206:                case 4:
                    207:                        for ( i = 0; i < nd_wpd; i++ ) {
                    208:                                u1 = d1[i]; u2 = d2[i];
                    209:                                if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0;
                    210:                                if ( (u1&0xf000000) < (u2&0xf000000) ) return 0;
                    211:                                if ( (u1&0xf00000) < (u2&0xf00000) ) return 0;
                    212:                                if ( (u1&0xf0000) < (u2&0xf0000) ) return 0;
                    213:                                if ( (u1&0xf000) < (u2&0xf000) ) return 0;
                    214:                                if ( (u1&0xf00) < (u2&0xf00) ) return 0;
                    215:                                if ( (u1&0xf0) < (u2&0xf0) ) return 0;
                    216:                                if ( (u1&0xf) < (u2&0xf) ) return 0;
                    217:                        }
                    218:                        return 1;
                    219:                        break;
                    220:                case 6:
                    221:                        for ( i = 0; i < nd_wpd; i++ ) {
                    222:                                u1 = d1[i]; u2 = d2[i];
                    223:                                if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0;
                    224:                                if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0;
                    225:                                if ( (u1&0x3f000) < (u2&0x3f000) ) return 0;
                    226:                                if ( (u1&0xfc0) < (u2&0xfc0) ) return 0;
                    227:                                if ( (u1&0x3f) < (u2&0x3f) ) return 0;
                    228:                        }
                    229:                        return 1;
                    230:                        break;
                    231:                case 8:
                    232:                        for ( i = 0; i < nd_wpd; i++ ) {
                    233:                                u1 = d1[i]; u2 = d2[i];
                    234:                                if ( (u1&0xff000000) < (u2&0xff000000) ) return 0;
                    235:                                if ( (u1&0xff0000) < (u2&0xff0000) ) return 0;
                    236:                                if ( (u1&0xff00) < (u2&0xff00) ) return 0;
                    237:                                if ( (u1&0xff) < (u2&0xff) ) return 0;
                    238:                        }
                    239:                        return 1;
                    240:                        break;
                    241:                case 16:
                    242:                        for ( i = 0; i < nd_wpd; i++ ) {
                    243:                                u1 = d1[i]; u2 = d2[i];
                    244:                                if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0;
                    245:                                if ( (u1&0xffff) < (u2&0xffff) ) return 0;
                    246:                        }
                    247:                        return 1;
                    248:                        break;
                    249:                case 32:
                    250:                        for ( i = 0; i < nd_wpd; i++ )
                    251:                                if ( d1[i] < d2[i] ) return 0;
                    252:                        return 1;
                    253:                        break;
                    254:                default:
                    255:                        for ( i = 0; i < nd_wpd; i++ ) {
                    256:                                u1 = d1[i]; u2 = d2[i];
                    257:                                for ( j = 0; j < nd_epw; j++ )
                    258:                                        if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0;
                    259:                        }
                    260:                        return 1;
                    261:        }
                    262: }
                    263:
                    264: void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
                    265: {
                    266:        unsigned int t1,t2,u,u1,u2;
                    267:        int i,j;
                    268:
                    269:        switch ( nd_bpe ) {
                    270:                case 4:
                    271:                        for ( i = 0; i < nd_wpd; i++ ) {
                    272:                                u1 = d1[i]; u2 = d2[i];
                    273:                                t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2;
                    274:                                t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2;
                    275:                                t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2;
                    276:                                t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2;
                    277:                                t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2;
                    278:                                t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2;
                    279:                                t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2;
                    280:                                t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2;
                    281:                                d[i] = u;
                    282:                        }
                    283:                        break;
                    284:                case 6:
                    285:                        for ( i = 0; i < nd_wpd; i++ ) {
                    286:                                u1 = d1[i]; u2 = d2[i];
                    287:                                t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2;
                    288:                                t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2;
                    289:                                t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2;
                    290:                                t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2;
                    291:                                t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2;
                    292:                                d[i] = u;
                    293:                        }
                    294:                        break;
                    295:                case 8:
                    296:                        for ( i = 0; i < nd_wpd; i++ ) {
                    297:                                u1 = d1[i]; u2 = d2[i];
                    298:                                t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2;
                    299:                                t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2;
                    300:                                t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2;
                    301:                                t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2;
                    302:                                d[i] = u;
                    303:                        }
                    304:                        break;
                    305:                case 16:
                    306:                        for ( i = 0; i < nd_wpd; i++ ) {
                    307:                                u1 = d1[i]; u2 = d2[i];
                    308:                                t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2;
                    309:                                t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2;
                    310:                                d[i] = u;
                    311:                        }
                    312:                        break;
                    313:                case 32:
                    314:                        for ( i = 0; i < nd_wpd; i++ ) {
                    315:                                u1 = d1[i]; u2 = d2[i];
                    316:                                d[i] = u1>u2?u1:u2;
                    317:                        }
                    318:                        break;
                    319:                default:
                    320:                        for ( i = 0; i < nd_wpd; i++ ) {
                    321:                                u1 = d1[i]; u2 = d2[i];
                    322:                                for ( j = 0, u = 0; j < nd_epw; j++ ) {
                    323:                                        t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2;
                    324:                                }
                    325:                                d[i] = u;
                    326:                        }
                    327:                        break;
                    328:        }
                    329: }
                    330:
                    331: int ndl_td(unsigned int *d)
                    332: {
                    333:        unsigned int t,u;
                    334:        int i,j;
                    335:
                    336:        for ( t = 0, i = 0; i < nd_wpd; i++ ) {
                    337:                u = d[i];
                    338:                for ( j = 0; j < nd_epw; j++, u>>=nd_bpe )
                    339:                        t += (u&nd_mask0);
                    340:        }
                    341:        return t;
                    342: }
                    343:
                    344: INLINE int ndl_compare(unsigned int *d1,unsigned int *d2)
                    345: {
                    346:        int i;
                    347:
                    348:        for ( i = 0; i < nd_wpd; i++, d1++, d2++ )
                    349:                if ( *d1 > *d2 )
                    350:                        return is_rlex ? -1 : 1;
                    351:                else if ( *d1 < *d2 )
                    352:                        return is_rlex ? 1 : -1;
                    353:        return 0;
                    354: }
                    355:
                    356: INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)
                    357: {
                    358:        int i;
                    359:
                    360:        for ( i = 0; i < nd_wpd; i++ )
                    361:                if ( d1[i] != d2[i] )
                    362:                        return 0;
                    363:        return 1;
                    364: }
                    365:
                    366: INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d)
                    367: {
                    368:        int i;
                    369:
                    370:        for ( i = 0; i < nd_wpd; i++ ) {
                    371:                d[i] = d1[i]+d2[i];
                    372:        }
                    373: }
                    374:
                    375: void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
                    376: {
                    377:        int i;
                    378:
                    379:        for ( i = 0; i < nd_wpd; i++ )
                    380:                d[i] = d1[i]-d2[i];
                    381: }
                    382:
                    383: int ndl_disjoint(unsigned int *d1,unsigned int *d2)
                    384: {
                    385:        unsigned int t1,t2,u,u1,u2;
                    386:        int i,j;
                    387:
                    388:        switch ( nd_bpe ) {
                    389:                case 4:
                    390:                        for ( i = 0; i < nd_wpd; i++ ) {
                    391:                                u1 = d1[i]; u2 = d2[i];
                    392:                                t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0;
                    393:                                t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0;
                    394:                                t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0;
                    395:                                t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0;
                    396:                                t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0;
                    397:                                t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0;
                    398:                                t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0;
                    399:                                t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0;
                    400:                        }
                    401:                        return 1;
                    402:                        break;
                    403:                case 6:
                    404:                        for ( i = 0; i < nd_wpd; i++ ) {
                    405:                                u1 = d1[i]; u2 = d2[i];
                    406:                                t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0;
                    407:                                t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0;
                    408:                                t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0;
                    409:                                t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0;
                    410:                                t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0;
                    411:                        }
                    412:                        return 1;
                    413:                        break;
                    414:                case 8:
                    415:                        for ( i = 0; i < nd_wpd; i++ ) {
                    416:                                u1 = d1[i]; u2 = d2[i];
                    417:                                t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0;
                    418:                                t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0;
                    419:                                t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0;
                    420:                                t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0;
                    421:                        }
                    422:                        return 1;
                    423:                        break;
                    424:                case 16:
                    425:                        for ( i = 0; i < nd_wpd; i++ ) {
                    426:                                u1 = d1[i]; u2 = d2[i];
                    427:                                t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0;
                    428:                                t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0;
                    429:                        }
                    430:                        return 1;
                    431:                        break;
                    432:                case 32:
                    433:                        for ( i = 0; i < nd_wpd; i++ )
                    434:                                if ( d1[i] && d2[i] ) return 0;
                    435:                        return 1;
                    436:                        break;
                    437:                default:
                    438:                        for ( i = 0; i < nd_wpd; i++ ) {
                    439:                                u1 = d1[i]; u2 = d2[i];
                    440:                                for ( j = 0; j < nd_epw; j++ ) {
                    441:                                        if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0;
                    442:                                        u1 >>= nd_bpe; u2 >>= nd_bpe;
                    443:                                }
                    444:                        }
                    445:                        return 1;
                    446:                        break;
                    447:        }
                    448: }
                    449:
                    450: ND nd_reduce(ND p1,ND p2)
                    451: {
                    452:        int c,c1,c2,t,td,td2,mul;
                    453:        NM m2,prev,head,cur,new;
                    454:        unsigned int *d;
                    455:
                    456:        if ( !p1 )
                    457:                return 0;
                    458:        else {
                    459:                c2 = invm(HC(p2),nd_mod);
                    460:                c1 = nd_mod-HC(p1);
                    461:                DMAR(c1,c2,0,nd_mod,mul);
                    462:                td = HTD(p1)-HTD(p2);
                    463:                d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                    464:                ndl_sub(HDL(p1),HDL(p2),d);
                    465:                prev = 0; head = cur = BDY(p1);
                    466:                NEWNM(new);
                    467:                for ( m2 = BDY(p2); m2; ) {
                    468:                        td2 = new->td = m2->td+td;
                    469:                        ndl_add(m2->dl,d,new->dl);
                    470:                        if ( !cur ) {
                    471:                                c1 = C(m2);
                    472:                                DMAR(c1,mul,0,nd_mod,c2);
                    473:                                C(new) = c2;
                    474:                                if ( !prev ) {
                    475:                                        prev = new;
                    476:                                        NEXT(prev) = 0;
                    477:                                        head = prev;
                    478:                                } else {
                    479:                                        NEXT(prev) = new;
                    480:                                        NEXT(new) = 0;
                    481:                                        prev = new;
                    482:                                }
                    483:                                m2 = NEXT(m2);
                    484:                                NEWNM(new);
                    485:                                continue;
                    486:                        }
                    487:                        if ( cur->td > td2 )
                    488:                                c = 1;
                    489:                        else if ( cur->td < td2 )
                    490:                                c = -1;
                    491:                        else
                    492:                                c = ndl_compare(cur->dl,new->dl);
                    493:                        switch ( c ) {
                    494:                                case 0:
                    495:                                        c2 = C(m2);
                    496:                                        c1 = C(cur);
                    497:                                        DMAR(c2,mul,c1,nd_mod,t);
                    498:                                        if ( t )
                    499:                                                C(cur) = t;
                    500:                                        else if ( !prev ) {
                    501:                                                head = NEXT(cur);
                    502:                                                FREENM(cur);
                    503:                                                cur = head;
                    504:                                        } else {
                    505:                                                NEXT(prev) = NEXT(cur);
                    506:                                                FREENM(cur);
                    507:                                                cur = NEXT(prev);
                    508:                                        }
                    509:                                        m2 = NEXT(m2);
                    510:                                        break;
                    511:                                case 1:
                    512:                                        prev = cur;
                    513:                                        cur = NEXT(cur);
                    514:                                        break;
                    515:                                case -1:
                    516:                                        if ( !prev ) {
                    517:                                                /* cur = head */
                    518:                                                prev = new;
                    519:                                                c2 = C(m2);
                    520:                                                DMAR(c2,mul,0,nd_mod,c1);
                    521:                                                C(prev) = c1;
                    522:                                                NEXT(prev) = head;
                    523:                                                head = prev;
                    524:                                        } else {
                    525:                                                c2 = C(m2);
                    526:                                                DMAR(c2,mul,0,nd_mod,c1);
                    527:                                                C(new) = c1;
                    528:                                                NEXT(prev) = new;
                    529:                                                NEXT(new) = cur;
                    530:                                                prev = new;
                    531:                                        }
                    532:                                        NEWNM(new);
                    533:                                        m2 = NEXT(m2);
                    534:                                        break;
                    535:                        }
                    536:                }
                    537:                FREENM(new);
                    538:                if ( head ) {
                    539:                        BDY(p1) = head;
                    540:                        p1->sugar = MAX(p1->sugar,p2->sugar+td);
                    541:                        return p1;
                    542:                } else {
                    543:                        FREEND(p1);
                    544:                        return 0;
                    545:                }
                    546:
                    547:        }
                    548: }
                    549:
1.3     ! noro      550: ND ndv_add(ND p1,NDV p2)
        !           551: {
        !           552:        int c,c1,c2,t,td,td2,mul,len,i;
        !           553:        NM prev,head,cur,new;
        !           554:        unsigned int *d;
        !           555:        NMV m2;
        !           556:
        !           557:        if ( !p1 )
        !           558:                return 0;
        !           559:        else {
        !           560:                prev = 0; head = cur = BDY(p1);
        !           561:                NEWNM(new);
        !           562:                len = p2->len;
        !           563:                for ( m2 = BDY(p2), i = 0; i < len; ) {
        !           564:                        td2 = new->td = m2->td;
        !           565:                        if ( !cur ) {
        !           566:                                C(new) = C(m2);
        !           567:                                bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
        !           568:                                if ( !prev ) {
        !           569:                                        prev = new;
        !           570:                                        NEXT(prev) = 0;
        !           571:                                        head = prev;
        !           572:                                } else {
        !           573:                                        NEXT(prev) = new;
        !           574:                                        NEXT(new) = 0;
        !           575:                                        prev = new;
        !           576:                                }
        !           577:                                NMV_ADV(m2); i++;
        !           578:                                NEWNM(new);
        !           579:                                continue;
        !           580:                        }
        !           581:                        if ( cur->td > td2 )
        !           582:                                c = 1;
        !           583:                        else if ( cur->td < td2 )
        !           584:                                c = -1;
        !           585:                        else
        !           586:                                c = ndl_compare(cur->dl,m2->dl);
        !           587:                        switch ( c ) {
        !           588:                                case 0:
        !           589:                                        t = C(m2)+C(cur)-nd_mod;
        !           590:                                        if ( t < 0 )
        !           591:                                                t += nd_mod;
        !           592:                                        if ( t )
        !           593:                                                C(cur) = t;
        !           594:                                        else if ( !prev ) {
        !           595:                                                head = NEXT(cur);
        !           596:                                                FREENM(cur);
        !           597:                                                cur = head;
        !           598:                                        } else {
        !           599:                                                NEXT(prev) = NEXT(cur);
        !           600:                                                FREENM(cur);
        !           601:                                                cur = NEXT(prev);
        !           602:                                        }
        !           603:                                        NMV_ADV(m2); i++;
        !           604:                                        break;
        !           605:                                case 1:
        !           606:                                        prev = cur;
        !           607:                                        cur = NEXT(cur);
        !           608:                                        break;
        !           609:                                case -1:
        !           610:                                        if ( !prev ) {
        !           611:                                                /* cur = head */
        !           612:                                                prev = new;
        !           613:                                                bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
        !           614:                                                C(prev) = C(m2);
        !           615:                                                NEXT(prev) = head;
        !           616:                                                head = prev;
        !           617:                                        } else {
        !           618:                                                C(new) = C(m2);
        !           619:                                                bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
        !           620:                                                NEXT(prev) = new;
        !           621:                                                NEXT(new) = cur;
        !           622:                                                prev = new;
        !           623:                                        }
        !           624:                                        NEWNM(new);
        !           625:                                        NMV_ADV(m2); i++;
        !           626:                                        break;
        !           627:                        }
        !           628:                }
        !           629:                FREENM(new);
        !           630:                if ( head ) {
        !           631:                        BDY(p1) = head;
        !           632:                        p1->sugar = MAX(p1->sugar,p2->sugar+td);
        !           633:                        return p1;
        !           634:                } else {
        !           635:                        FREEND(p1);
        !           636:                        return 0;
        !           637:                }
        !           638:
        !           639:        }
        !           640: }
        !           641:
1.1       noro      642: /* HDL(p1) = HDL(p2) */
                    643:
                    644: ND nd_reduce_special(ND p1,ND p2)
                    645: {
                    646:        int c,c1,c2,t,td,td2,mul;
                    647:        NM m2,prev,head,cur,new;
                    648:
                    649:        if ( !p1 )
                    650:                return 0;
                    651:        else {
                    652:                c2 = invm(HC(p2),nd_mod);
                    653:                c1 = nd_mod-HC(p1);
                    654:                DMAR(c1,c2,0,nd_mod,mul);
                    655:                prev = 0; head = cur = BDY(p1);
                    656:                NEWNM(new);
                    657:                for ( m2 = BDY(p2); m2; ) {
                    658:                        td2 = new->td = m2->td;
                    659:                        if ( !cur ) {
                    660:                                c1 = C(m2);
                    661:                                DMAR(c1,mul,0,nd_mod,c2);
                    662:                                C(new) = c2;
                    663:                                bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
                    664:                                if ( !prev ) {
                    665:                                        prev = new;
                    666:                                        NEXT(prev) = 0;
                    667:                                        head = prev;
                    668:                                } else {
                    669:                                        NEXT(prev) = new;
                    670:                                        NEXT(new) = 0;
                    671:                                        prev = new;
                    672:                                }
                    673:                                m2 = NEXT(m2);
                    674:                                NEWNM(new);
                    675:                                continue;
                    676:                        }
                    677:                        if ( cur->td > td2 )
                    678:                                c = 1;
                    679:                        else if ( cur->td < td2 )
                    680:                                c = -1;
                    681:                        else
                    682:                                c = ndl_compare(cur->dl,m2->dl);
                    683:                        switch ( c ) {
                    684:                                case 0:
                    685:                                        c2 = C(m2);
                    686:                                        c1 = C(cur);
                    687:                                        DMAR(c2,mul,c1,nd_mod,t);
                    688:                                        if ( t )
                    689:                                                C(cur) = t;
                    690:                                        else if ( !prev ) {
                    691:                                                head = NEXT(cur);
                    692:                                                FREENM(cur);
                    693:                                                cur = head;
                    694:                                        } else {
                    695:                                                NEXT(prev) = NEXT(cur);
                    696:                                                FREENM(cur);
                    697:                                                cur = NEXT(prev);
                    698:                                        }
                    699:                                        m2 = NEXT(m2);
                    700:                                        break;
                    701:                                case 1:
                    702:                                        prev = cur;
                    703:                                        cur = NEXT(cur);
                    704:                                        break;
                    705:                                case -1:
                    706:                                        bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
                    707:                                        if ( !prev ) {
                    708:                                                /* cur = head */
                    709:                                                prev = new;
                    710:                                                c2 = C(m2);
                    711:                                                DMAR(c2,mul,0,nd_mod,c1);
                    712:                                                C(prev) = c1;
                    713:                                                NEXT(prev) = head;
                    714:                                                head = prev;
                    715:                                        } else {
                    716:                                                c2 = C(m2);
                    717:                                                DMAR(c2,mul,0,nd_mod,c1);
                    718:                                                C(new) = c1;
                    719:                                                NEXT(prev) = new;
                    720:                                                NEXT(new) = cur;
                    721:                                                prev = new;
                    722:                                        }
                    723:                                        NEWNM(new);
                    724:                                        m2 = NEXT(m2);
                    725:                                        break;
                    726:                        }
                    727:                }
                    728:                FREENM(new);
                    729:                if ( head ) {
                    730:                        BDY(p1) = head;
                    731:                        p1->sugar = MAX(p1->sugar,p2->sugar+td);
                    732:                        return p1;
                    733:                } else {
                    734:                        FREEND(p1);
                    735:                        return 0;
                    736:                }
                    737:
                    738:        }
                    739: }
                    740:
                    741: INLINE int ndl_check_bound(unsigned int *d)
                    742: {
                    743:        int i;
                    744:
                    745:        for ( i = 0; i < nd_wpd; i++ )
                    746:                if ( d[i] & nd_mask1 )
                    747:                        return 1;
                    748:        return 0;
                    749: }
                    750:
                    751: int nd_sp(ND_pairs p,ND *rp)
                    752: {
                    753:        NM m;
                    754:        ND p1,p2,t1,t2;
                    755:        unsigned int *lcm,*check;
                    756:        int td;
                    757:
                    758:        check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                    759:        p1 = nd_ps[p->i1];
                    760:        p2 = nd_ps[p->i2];
                    761:        lcm = p->lcm;
                    762:        td = p->td;
                    763:        NEWNM(m);
                    764:        C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl); NEXT(m) = 0;
                    765:        ndl_add(nd_bound[p->i1],m->dl,check);
                    766:        if ( ndl_check_bound(check) )
                    767:                return 0;
                    768:        t1 = nd_mul_nm(p1,m);
                    769:        C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
                    770:        ndl_add(nd_bound[p->i2],m->dl,check);
                    771:        if ( ndl_check_bound(check) ) {
                    772:                nd_free(t1);
                    773:                return 0;
                    774:        }
                    775:        t2 = nd_mul_nm(p2,m);
                    776:        FREENM(m);
                    777:        *rp = nd_add(t1,t2);
                    778:        return 1;
                    779: }
                    780:
                    781: int ndl_hash_value(int td,unsigned int *d)
                    782: {
                    783:        int i;
                    784:        int r;
                    785:
                    786:        r = td;
                    787:        for ( i = 0; i < nd_wpd; i++ )
                    788:                r = ((r<<16)+d[i])%REDTAB_LEN;
                    789:        return r;
                    790: }
                    791:
                    792: int nd_find_reducer(ND g, ND *rp)
                    793: {
                    794:        NM m;
                    795:        ND r,p;
                    796:        int i,c1,c2,c;
                    797:        int d,k,append,index;
                    798:        unsigned int *check;
                    799:        NM t;
                    800:
                    801:        d = ndl_hash_value(HTD(g),HDL(g));
                    802:        for ( m = nd_red[d], k = 0; m; m = NEXT(m), k++ ) {
                    803:                if ( HTD(g) == m->td && ndl_equal(HDL(g),m->dl) ) {
                    804:                        if ( k > 0 ) nd_notfirst++;
                    805:                        index = m->c;
                    806:                        append = 0;
                    807:                        nd_found++;
                    808:                        goto found;
                    809:                }
                    810:        }
                    811:
                    812:        for ( i = 0; i < nd_psn; i++ ) {
                    813:                p = nd_ps[i];
                    814:                if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
                    815:                        index = i;
                    816:                        append = 1;
                    817:                        nd_create++;
                    818:                        goto found;
                    819:                }
                    820:        }
                    821:        return 0;
                    822:
                    823: found:
                    824:        NEWNM(m);
                    825:        p = nd_ps[index];
                    826:        ndl_sub(HDL(g),HDL(p),m->dl);
                    827:
                    828:        check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                    829:        ndl_add(nd_bound[index],m->dl,check);
                    830:        if ( ndl_check_bound(check) ) {
                    831:                FREENM(m);
                    832:                return -1;
                    833:        }
                    834:
                    835:        c1 = invm(HC(p),nd_mod);
                    836:        c2 = nd_mod-HC(g);
                    837:        DMAR(c1,c2,0,nd_mod,c);
                    838:        C(m) = c;
                    839:        m->td = HTD(g)-HTD(p);
                    840:        NEXT(m) = 0;
                    841:        *rp = r = nd_mul_nm(p,m);
                    842:        FREENM(m);
                    843:
                    844:        if ( append ) nd_append_red(HDL(g),HTD(g),i);
                    845:        return 1;
                    846: }
                    847:
1.3     ! noro      848: int ndv_find_reducer(ND g, NDV red)
        !           849: {
        !           850:        NM m;
        !           851:        ND r,p;
        !           852:        int i,c1,c2,c;
        !           853:        int d,k,append,index;
        !           854:        unsigned int *check;
        !           855:        NM t;
        !           856:        NDV pv;
        !           857:
        !           858:        d = ndl_hash_value(HTD(g),HDL(g));
        !           859:        for ( m = nd_red[d], k = 0; m; m = NEXT(m), k++ ) {
        !           860:                if ( HTD(g) == m->td && ndl_equal(HDL(g),m->dl) ) {
        !           861:                        if ( k > 0 ) nd_notfirst++;
        !           862:                        index = m->c;
        !           863:                        append = 0;
        !           864:                        nd_found++;
        !           865:                        goto found;
        !           866:                }
        !           867:        }
        !           868:
        !           869:        for ( i = 0; i < nd_psn; i++ ) {
        !           870:                p = nd_ps[i];
        !           871:                if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
        !           872:                        index = i;
        !           873:                        append = 1;
        !           874:                        nd_create++;
        !           875:                        goto found;
        !           876:                }
        !           877:        }
        !           878:        return 0;
        !           879:
        !           880: found:
        !           881:        NEWNM(m);
        !           882:        p = nd_ps[index];
        !           883:        pv = nd_psv[index];
        !           884:        ndl_sub(HDL(g),HDL(p),m->dl);
        !           885:
        !           886:        check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
        !           887:        ndl_add(nd_bound[index],m->dl,check);
        !           888:        if ( ndl_check_bound(check) ) {
        !           889:                FREENM(m);
        !           890:                return -1;
        !           891:        }
        !           892:
        !           893:        c1 = invm(HC(p),nd_mod);
        !           894:        c2 = nd_mod-HC(g);
        !           895:        DMAR(c1,c2,0,nd_mod,c);
        !           896:        C(m) = c;
        !           897:        m->td = HTD(g)-HTD(p);
        !           898:        NEXT(m) = 0;
        !           899:        ndv_mul_nm(pv,m,red);
        !           900:        FREENM(m);
        !           901:
        !           902:        if ( append ) nd_append_red(HDL(g),HTD(g),i);
        !           903:        return 1;
        !           904: }
        !           905:
1.1       noro      906: ND nd_find_monic_reducer(ND g)
                    907: {
                    908:        int *d;
                    909:        ND p,r;
                    910:        int i;
                    911:
                    912:        for ( i = 0; i < nd_psn; i++ ) {
                    913:                p = nd_ps[i];
                    914:                if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
                    915:                        d = (int *)ALLOCA(nd_wpd*sizeof(int));
                    916:                        ndl_sub(HDL(g),HDL(p),d);
                    917:                        r = nd_mul_term(p,HTD(g)-HTD(p),d);
                    918:                        return r;
                    919:                }
                    920:        }
                    921:        return 0;
                    922: }
                    923:
                    924: ND nd_add(ND p1,ND p2)
                    925: {
                    926:        int n,c;
                    927:        int t;
                    928:        ND r;
                    929:        NM m1,m2,mr0,mr,s;
                    930:
                    931:        if ( !p1 )
                    932:                return p2;
                    933:        else if ( !p2 )
                    934:                return p1;
                    935:        else {
                    936:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
                    937:                        if ( m1->td > m2->td )
                    938:                                c = 1;
                    939:                        else if ( m1->td < m2->td )
                    940:                                c = -1;
                    941:                        else
                    942:                                c = ndl_compare(m1->dl,m2->dl);
                    943:                        switch ( c ) {
                    944:                                case 0:
                    945:                                        t = ((C(m1))+(C(m2))) - nd_mod;
                    946:                                        if ( t < 0 )
                    947:                                                t += nd_mod;
                    948:                                        s = m1; m1 = NEXT(m1);
                    949:                                        if ( t ) {
                    950:                                                NEXTNM2(mr0,mr,s); C(mr) = (t);
                    951:                                        } else {
                    952:                                                FREENM(s);
                    953:                                        }
                    954:                                        s = m2; m2 = NEXT(m2); FREENM(s);
                    955:                                        break;
                    956:                                case 1:
                    957:                                        s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
                    958:                                        break;
                    959:                                case -1:
                    960:                                        s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
                    961:                                        break;
                    962:                        }
                    963:                }
                    964:                if ( !mr0 )
                    965:                        if ( m1 )
                    966:                                mr0 = m1;
                    967:                        else if ( m2 )
                    968:                                mr0 = m2;
                    969:                        else
                    970:                                return 0;
                    971:                else if ( m1 )
                    972:                        NEXT(mr) = m1;
                    973:                else if ( m2 )
                    974:                        NEXT(mr) = m2;
                    975:                else
                    976:                        NEXT(mr) = 0;
                    977:                BDY(p1) = mr0;
                    978:                p1->sugar = MAX(p1->sugar,p2->sugar);
                    979:                FREEND(p2);
                    980:                return p1;
                    981:        }
                    982: }
                    983:
                    984: ND nd_mul_nm(ND p,NM m0)
                    985: {
                    986:        NM m,mr,mr0;
                    987:        unsigned int *d,*dt,*dm;
                    988:        int c,n,td,i,c1,c2;
                    989:        int *pt,*p1,*p2;
                    990:        ND r;
                    991:
                    992:        if ( !p )
                    993:                return 0;
                    994:        else {
                    995:                n = NV(p); m = BDY(p);
                    996:                d = m0->dl; td = m0->td; c = C(m0);
                    997:                mr0 = 0;
                    998:                for ( ; m; m = NEXT(m) ) {
                    999:                        NEXTNM(mr0,mr);
                   1000:                        c1 = C(m);
                   1001:                        DMAR(c1,c,0,nd_mod,c2);
                   1002:                        C(mr) = c2;
                   1003:                        mr->td = m->td+td;
                   1004:                        ndl_add(m->dl,d,mr->dl);
                   1005:                }
                   1006:                NEXT(mr) = 0;
                   1007:                MKND(NV(p),mr0,r);
                   1008:                r->sugar = p->sugar + td;
                   1009:                return r;
                   1010:        }
                   1011: }
                   1012:
1.3     ! noro     1013: void ndv_mul_nm(NDV p,NM m0,NDV r)
        !          1014: {
        !          1015:        NMV m,mr,mr0;
        !          1016:        unsigned int *d,*dt,*dm;
        !          1017:        int c,n,td,i,c1,c2,len;
        !          1018:
        !          1019:        if ( !p )
        !          1020:                /* XXX */
        !          1021:                r->len = 0;
        !          1022:        else {
        !          1023:                n = NV(p); m = BDY(p); len = p->len;
        !          1024:                d = m0->dl; td = m0->td; c = C(m0);
        !          1025:                mr = BDY(r);
        !          1026:                for ( i = 0; i < len; i++, NMV_ADV(m), NMV_ADV(mr) ) {
        !          1027:                        c1 = C(m);
        !          1028:                        DMAR(c1,c,0,nd_mod,c2);
        !          1029:                        C(mr) = c2;
        !          1030:                        mr->td = m->td+td;
        !          1031:                        ndl_add(m->dl,d,mr->dl);
        !          1032:                }
        !          1033:                NV(r) = NV(p);
        !          1034:                r->len = p->len;
        !          1035:                r->sugar = p->sugar + td;
        !          1036:        }
        !          1037: }
        !          1038:
1.1       noro     1039: ND nd_mul_term(ND p,int td,unsigned int *d)
                   1040: {
                   1041:        NM m,mr,mr0;
                   1042:        int c,n;
                   1043:        ND r;
                   1044:
                   1045:        if ( !p )
                   1046:                return 0;
                   1047:        else {
                   1048:                n = NV(p); m = BDY(p);
                   1049:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   1050:                        NEXTNM(mr0,mr);
                   1051:                        C(mr) = C(m);
                   1052:                        mr->td = m->td+td;
                   1053:                        ndl_add(m->dl,d,mr->dl);
                   1054:                }
                   1055:                NEXT(mr) = 0;
                   1056:                MKND(NV(p),mr0,r);
                   1057:                r->sugar = p->sugar + td;
                   1058:                return r;
                   1059:        }
                   1060: }
                   1061:
                   1062: #if 1
                   1063: /* ret=1 : success, ret=0 : overflow */
                   1064: int nd_nf(ND g,int full,ND *rp)
                   1065: {
1.3     ! noro     1066:        ND p,d;
1.1       noro     1067:        NM m,mrd,tail;
                   1068:        int n,sugar,psugar,stat;
1.3     ! noro     1069:        NDV red;
1.1       noro     1070:
                   1071:        if ( !g ) {
                   1072:                *rp = 0;
                   1073:                return 1;
                   1074:        }
                   1075:        sugar = g->sugar;
                   1076:        n = NV(g);
                   1077:        for ( d = 0; g; ) {
                   1078:                /* stat=1 : found, stat=0 : not found, stat=-1 : overflow */
1.3     ! noro     1079: #if 0
1.1       noro     1080:                stat = nd_find_reducer(g,&red);
1.3     ! noro     1081: #else
        !          1082:                stat = ndv_find_reducer(g,ndv_red);
        !          1083: #endif
1.1       noro     1084:                if ( stat == 1 ) {
1.3     ! noro     1085: #if 0
1.1       noro     1086:                        g = nd_add(g,red);
                   1087: #else
1.3     ! noro     1088:                        g = ndv_add(g,ndv_red);
1.1       noro     1089: #endif
1.3     ! noro     1090:                sugar = MAX(sugar,ndv_red->sugar);
1.1       noro     1091:                } else if ( stat == -1 ) {
                   1092:                        nd_free(g);
                   1093:                        nd_free(d);
                   1094:                        return 0;
                   1095:                } else if ( !full ) {
                   1096:                        *rp = g;
                   1097:                        return 1;
                   1098:                } else {
                   1099:                        m = BDY(g);
                   1100:                        if ( NEXT(m) ) {
                   1101:                                BDY(g) = NEXT(m); NEXT(m) = 0;
                   1102:                        } else {
                   1103:                                FREEND(g); g = 0;
                   1104:                        }
                   1105:                        if ( d ) {
                   1106:                                NEXT(tail)=m;
                   1107:                                tail=m;
                   1108:                        } else {
                   1109:                                MKND(n,m,d);
                   1110:                                tail = BDY(d);
                   1111:                        }
                   1112:                }
                   1113:        }
                   1114:        if ( d )
                   1115:                d->sugar = sugar;
                   1116:        *rp = d;
                   1117:        return 1;
                   1118: }
                   1119: #else
                   1120:
                   1121: ND nd_remove_head(ND p)
                   1122: {
                   1123:        NM m;
                   1124:
                   1125:        m = BDY(p);
                   1126:        if ( !NEXT(m) ) {
                   1127:                FREEND(p);
                   1128:                p = 0;
                   1129:        } else
                   1130:                BDY(p) = NEXT(m);
                   1131:        FREENM(m);
                   1132:        return p;
                   1133: }
                   1134:
                   1135: PGeoBucket create_pbucket()
                   1136: {
                   1137:     PGeoBucket g;
                   1138:
                   1139:        g = CALLOC(1,sizeof(struct oPGeoBucket));
                   1140:        g->m = -1;
                   1141:        return g;
                   1142: }
                   1143:
                   1144: void add_pbucket(PGeoBucket g,ND d)
                   1145: {
                   1146:        int l,k,m;
                   1147:
                   1148:        l = nd_length(d);
                   1149:        for ( k = 0, m = 1; l > m; k++, m <<= 2 );
                   1150:        /* 4^(k-1) < l <= 4^k */
                   1151:        d = nd_add(g->body[k],d);
                   1152:        for ( ; d && nd_length(d) > 1<<(2*k); k++ ) {
                   1153:                g->body[k] = 0;
                   1154:                d = nd_add(g->body[k+1],d);
                   1155:        }
                   1156:        g->body[k] = d;
                   1157:        g->m = MAX(g->m,k);
                   1158: }
                   1159:
                   1160: int head_pbucket(PGeoBucket g)
                   1161: {
                   1162:        int j,i,c,k,nv,sum;
                   1163:        unsigned int *di,*dj;
                   1164:        ND gi,gj;
                   1165:
                   1166:        k = g->m;
                   1167:        while ( 1 ) {
                   1168:                j = -1;
                   1169:                for ( i = 0; i <= k; i++ ) {
                   1170:                        if ( !(gi = g->body[i]) )
                   1171:                                continue;
                   1172:                        if ( j < 0 ) {
                   1173:                                j = i;
                   1174:                                gj = g->body[j];
                   1175:                                dj = HDL(gj);
                   1176:                                sum = HC(gj);
                   1177:                        } else {
                   1178:                                di = HDL(gi);
                   1179:                                nv = NV(gi);
                   1180:                                if ( HTD(gi) > HTD(gj) )
                   1181:                                        c = 1;
                   1182:                                else if ( HTD(gi) < HTD(gj) )
                   1183:                                        c = -1;
                   1184:                                else
                   1185:                                        c = ndl_compare(di,dj);
                   1186:                                if ( c > 0 ) {
                   1187:                                        if ( sum )
                   1188:                                                HC(gj) = sum;
                   1189:                                        else
                   1190:                                                g->body[j] = nd_remove_head(gj);
                   1191:                                        j = i;
                   1192:                                        gj = g->body[j];
                   1193:                                        dj = HDL(gj);
                   1194:                                        sum = HC(gj);
                   1195:                                } else if ( c == 0 ) {
                   1196:                                        sum = sum+HC(gi)-nd_mod;
                   1197:                                        if ( sum < 0 )
                   1198:                                                sum += nd_mod;
                   1199:                                        g->body[i] = nd_remove_head(gi);
                   1200:                                }
                   1201:                        }
                   1202:                }
                   1203:                if ( j < 0 )
                   1204:                        return -1;
                   1205:                else if ( sum ) {
                   1206:                        HC(gj) = sum;
                   1207:                        return j;
                   1208:                } else
                   1209:                        g->body[j] = nd_remove_head(gj);
                   1210:        }
                   1211: }
                   1212:
                   1213: ND normalize_pbucket(PGeoBucket g)
                   1214: {
                   1215:        int i;
                   1216:        ND r,t;
                   1217:
                   1218:        r = 0;
                   1219:        for ( i = 0; i <= g->m; i++ )
                   1220:                r = nd_add(r,g->body[i]);
                   1221:        return r;
                   1222: }
                   1223:
                   1224: ND nd_nf(ND g,int full)
                   1225: {
                   1226:        ND u,p,d,red;
                   1227:        NODE l;
                   1228:        NM m,mrd;
                   1229:        int sugar,psugar,n,h_reducible,h;
                   1230:        PGeoBucket bucket;
                   1231:
                   1232:        if ( !g ) {
                   1233:                return 0;
                   1234:        }
                   1235:        sugar = g->sugar;
                   1236:        n = g->nv;
                   1237:        bucket = create_pbucket();
                   1238:        add_pbucket(bucket,g);
                   1239:        d = 0;
                   1240:        while ( 1 ) {
                   1241:                h = head_pbucket(bucket);
                   1242:                if ( h < 0 ) {
                   1243:                        if ( d )
                   1244:                                d->sugar = sugar;
                   1245:                        return d;
                   1246:                }
                   1247:                g = bucket->body[h];
                   1248:                red = nd_find_reducer(g);
                   1249:                if ( red ) {
                   1250:                        bucket->body[h] = nd_remove_head(g);
                   1251:                        red = nd_remove_head(red);
                   1252:                        add_pbucket(bucket,red);
                   1253:                        sugar = MAX(sugar,red->sugar);
                   1254:                } else if ( !full ) {
                   1255:                        g = normalize_pbucket(bucket);
                   1256:                        if ( g )
                   1257:                                g->sugar = sugar;
                   1258:                        return g;
                   1259:                } else {
                   1260:                        m = BDY(g);
                   1261:                        if ( NEXT(m) ) {
                   1262:                                BDY(g) = NEXT(m); NEXT(m) = 0;
                   1263:                        } else {
                   1264:                                FREEND(g); g = 0;
                   1265:                        }
                   1266:                        bucket->body[h] = g;
                   1267:                        NEXT(m) = 0;
                   1268:                        if ( d ) {
                   1269:                                for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
                   1270:                                NEXT(mrd) = m;
                   1271:                        } else {
                   1272:                                MKND(n,m,d);
                   1273:                        }
                   1274:                }
                   1275:        }
                   1276: }
                   1277: #endif
                   1278:
                   1279: NODE nd_gb(NODE f)
                   1280: {
                   1281:        int i,nh,sugar,stat;
                   1282:        NODE r,g,gall;
                   1283:        ND_pairs d;
                   1284:        ND_pairs l;
                   1285:        ND h,nf;
                   1286:
                   1287:        for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                   1288:                i = (int)BDY(r);
                   1289:                d = update_pairs(d,g,i);
                   1290:                g = update_base(g,i);
                   1291:                gall = append_one(gall,i);
                   1292:        }
                   1293:        sugar = 0;
                   1294:        while ( d ) {
                   1295: again:
                   1296:                l = nd_minp(d,&d);
                   1297:                if ( l->sugar != sugar ) {
                   1298:                        sugar = l->sugar;
                   1299:                        fprintf(asir_out,"%d",sugar);
                   1300:                }
                   1301:                stat = nd_sp(l,&h);
                   1302:                if ( !stat ) {
                   1303:                        NEXT(l) = d; d = l;
                   1304:                        d = nd_reconstruct(d);
                   1305:                        goto again;
                   1306:                }
                   1307:                stat = nd_nf(h,!Top,&nf);
                   1308:                if ( !stat ) {
                   1309:                        NEXT(l) = d; d = l;
                   1310:                        d = nd_reconstruct(d);
                   1311:                        goto again;
                   1312:                } else if ( nf ) {
                   1313:                        printf("+"); fflush(stdout);
                   1314:                        nh = nd_newps(nf);
                   1315:                        d = update_pairs(d,g,nh);
                   1316:                        g = update_base(g,nh);
                   1317:                        gall = append_one(gall,nh);
                   1318:                        FREENDP(l);
                   1319:                } else {
                   1320:                        printf("."); fflush(stdout);
                   1321:                        FREENDP(l);
                   1322:                }
                   1323:        }
                   1324:        return g;
                   1325: }
                   1326:
                   1327: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t)
                   1328: {
                   1329:        ND_pairs d1,nd,cur,head,prev,remove;
                   1330:
                   1331:        if ( !g ) return d;
                   1332:        d = crit_B(d,t);
                   1333:        d1 = nd_newpairs(g,t);
                   1334:        d1 = crit_M(d1);
                   1335:        d1 = crit_F(d1);
                   1336:        prev = 0; cur = head = d1;
                   1337:        while ( cur ) {
                   1338:                if ( crit_2( cur->i1,cur->i2 ) ) {
                   1339:                        remove = cur;
                   1340:                        if ( !prev ) {
                   1341:                                head = cur = NEXT(cur);
                   1342:                        } else {
                   1343:                                cur = NEXT(prev) = NEXT(cur);
                   1344:                        }
                   1345:                        FREENDP(remove);
                   1346:                } else {
                   1347:                        prev = cur;
                   1348:                        cur = NEXT(cur);
                   1349:                }
                   1350:        }
                   1351:        if ( !d )
                   1352:                return head;
                   1353:        else {
                   1354:                nd = d;
                   1355:                while ( NEXT(nd) )
                   1356:                        nd = NEXT(nd);
                   1357:                NEXT(nd) = head;
                   1358:                return d;
                   1359:        }
                   1360: }
                   1361:
                   1362: ND_pairs nd_newpairs( NODE g, int t )
                   1363: {
                   1364:        NODE h;
                   1365:        unsigned int *dl;
                   1366:        int td,ts,s;
                   1367:        ND_pairs r,r0;
                   1368:
                   1369:        dl = HDL(nd_ps[t]);
                   1370:        td = HTD(nd_ps[t]);
                   1371:        ts = nd_ps[t]->sugar - td;
                   1372:        for ( r0 = 0, h = g; h; h = NEXT(h) ) {
                   1373:                NEXTND_pairs(r0,r);
                   1374:                r->i1 = (int)BDY(h);
                   1375:                r->i2 = t;
                   1376:                ndl_lcm(HDL(nd_ps[r->i1]),dl,r->lcm);
                   1377:                r->td = ndl_td(r->lcm);
                   1378:                s = nd_ps[r->i1]->sugar-HTD(nd_ps[r->i1]);
                   1379:                r->sugar = MAX(s,ts) + r->td;
                   1380:        }
                   1381:        NEXT(r) = 0;
                   1382:        return r0;
                   1383: }
                   1384:
                   1385: ND_pairs crit_B( ND_pairs d, int s )
                   1386: {
                   1387:        ND_pairs cur,head,prev,remove;
                   1388:        unsigned int *t,*tl,*lcm;
                   1389:        int td,tdl;
                   1390:
                   1391:        if ( !d ) return 0;
                   1392:        t = HDL(nd_ps[s]);
                   1393:        prev = 0;
                   1394:        head = cur = d;
                   1395:        lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                   1396:        while ( cur ) {
                   1397:                tl = cur->lcm;
                   1398:                if ( ndl_reducible(tl,t)
                   1399:                        && (ndl_lcm(HDL(nd_ps[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
                   1400:                        && (ndl_lcm(HDL(nd_ps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
                   1401:                        remove = cur;
                   1402:                        if ( !prev ) {
                   1403:                                head = cur = NEXT(cur);
                   1404:                        } else {
                   1405:                                cur = NEXT(prev) = NEXT(cur);
                   1406:                        }
                   1407:                        FREENDP(remove);
                   1408:                } else {
                   1409:                        prev = cur;
                   1410:                        cur = NEXT(cur);
                   1411:                }
                   1412:        }
                   1413:        return head;
                   1414: }
                   1415:
                   1416: ND_pairs crit_M( ND_pairs d1 )
                   1417: {
                   1418:        ND_pairs e,d2,d3,dd,p;
                   1419:        unsigned int *id,*jd;
                   1420:        int itd,jtd;
                   1421:
                   1422:        for ( dd = 0, e = d1; e; e = d3 ) {
                   1423:                if ( !(d2 = NEXT(e)) ) {
                   1424:                        NEXT(e) = dd;
                   1425:                        return e;
                   1426:                }
                   1427:                id = e->lcm;
                   1428:                itd = e->td;
                   1429:                for ( d3 = 0; d2; d2 = p ) {
                   1430:                        p = NEXT(d2),
                   1431:                        jd = d2->lcm;
                   1432:                        jtd = d2->td;
                   1433:                        if ( jtd == itd  )
                   1434:                                if ( id == jd );
                   1435:                                else if ( ndl_reducible(jd,id) ) continue;
                   1436:                                else if ( ndl_reducible(id,jd) ) goto delit;
                   1437:                                else ;
                   1438:                        else if ( jtd > itd )
                   1439:                                if ( ndl_reducible(jd,id) ) continue;
                   1440:                                else ;
                   1441:                        else if ( ndl_reducible(id,jd ) ) goto delit;
                   1442:                        NEXT(d2) = d3;
                   1443:                        d3 = d2;
                   1444:                }
                   1445:                NEXT(e) = dd;
                   1446:                dd = e;
                   1447:                continue;
                   1448:                /**/
                   1449:        delit:  NEXT(d2) = d3;
                   1450:                d3 = d2;
                   1451:                for ( ; p; p = d2 ) {
                   1452:                        d2 = NEXT(p);
                   1453:                        NEXT(p) = d3;
                   1454:                        d3 = p;
                   1455:                }
                   1456:                FREENDP(e);
                   1457:        }
                   1458:        return dd;
                   1459: }
                   1460:
                   1461: ND_pairs crit_F( ND_pairs d1 )
                   1462: {
                   1463:        ND_pairs rest, head,remove;
                   1464:        ND_pairs last, p, r, w;
                   1465:        int s;
                   1466:
                   1467:        for ( head = last = 0, p = d1; NEXT(p); ) {
                   1468:                r = w = equivalent_pairs(p,&rest);
                   1469:                s = r->sugar;
                   1470:                w = NEXT(w);
                   1471:                while ( w ) {
                   1472:                        if ( crit_2(w->i1,w->i2) ) {
                   1473:                                r = w;
                   1474:                                w = NEXT(w);
                   1475:                                while ( w ) {
                   1476:                                        remove = w;
                   1477:                                        w = NEXT(w);
                   1478:                                        FREENDP(remove);
                   1479:                                }
                   1480:                                break;
                   1481:                        } else if ( w->sugar < s ) {
                   1482:                                FREENDP(r);
                   1483:                                r = w;
                   1484:                                s = r->sugar;
                   1485:                                w = NEXT(w);
                   1486:                        } else {
                   1487:                                remove = w;
                   1488:                                w = NEXT(w);
                   1489:                                FREENDP(remove);
                   1490:                        }
                   1491:                }
                   1492:                if ( last ) NEXT(last) = r;
                   1493:                else head = r;
                   1494:                NEXT(last = r) = 0;
                   1495:                p = rest;
                   1496:                if ( !p ) return head;
                   1497:        }
                   1498:        if ( !last ) return p;
                   1499:        NEXT(last) = p;
                   1500:        return head;
                   1501: }
                   1502:
                   1503: int crit_2( int dp1, int dp2 )
                   1504: {
                   1505:        return ndl_disjoint(HDL(nd_ps[dp1]),HDL(nd_ps[dp2]));
                   1506: }
                   1507:
                   1508: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest )
                   1509: {
                   1510:        ND_pairs w,p,r,s;
                   1511:        unsigned int *d;
                   1512:        int td;
                   1513:
                   1514:        w = d1;
                   1515:        d = w->lcm;
                   1516:        td = w->td;
                   1517:        s = NEXT(w);
                   1518:        NEXT(w) = 0;
                   1519:        for ( r = 0; s; s = p ) {
                   1520:                p = NEXT(s);
                   1521:                if ( td == s->td && ndl_equal(d,s->lcm) ) {
                   1522:                        NEXT(s) = w;
                   1523:                        w = s;
                   1524:                } else {
                   1525:                        NEXT(s) = r;
                   1526:                        r = s;
                   1527:                }
                   1528:        }
                   1529:        *prest = r;
                   1530:        return w;
                   1531: }
                   1532:
                   1533: NODE update_base(NODE nd,int ndp)
                   1534: {
                   1535:        unsigned int *dl, *dln;
                   1536:        NODE last, p, head;
                   1537:        int td,tdn;
                   1538:
                   1539:        dl = HDL(nd_ps[ndp]);
                   1540:        td = HTD(nd_ps[ndp]);
                   1541:        for ( head = last = 0, p = nd; p; ) {
                   1542:                dln = HDL(nd_ps[(int)BDY(p)]);
                   1543:                tdn = HTD(nd_ps[(int)BDY(p)]);
                   1544:                if ( tdn >= td && ndl_reducible( dln, dl ) ) {
                   1545:                        p = NEXT(p);
                   1546:                        if ( last ) NEXT(last) = p;
                   1547:                } else {
                   1548:                        if ( !last ) head = p;
                   1549:                        p = NEXT(last = p);
                   1550:                }
                   1551:        }
                   1552:        head = append_one(head,ndp);
                   1553:        return head;
                   1554: }
                   1555:
                   1556: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
                   1557: {
                   1558:        ND_pairs m,ml,p,l;
                   1559:        unsigned int *lcm;
                   1560:        int s,td,len,tlen,c;
                   1561:
                   1562:        if ( !(p = NEXT(m = d)) ) {
                   1563:                *prest = p;
                   1564:                NEXT(m) = 0;
                   1565:                return m;
                   1566:        }
                   1567:        lcm = m->lcm;
                   1568:        s = m->sugar;
                   1569:        td = m->td;
                   1570:        len = nd_length(nd_ps[m->i1])+nd_length(nd_ps[m->i2]);
                   1571:        for ( ml = 0, l = m; p; p = NEXT(l = p) ) {
                   1572:                if (p->sugar < s)
                   1573:                        goto find;
                   1574:                else if ( p->sugar == s ) {
                   1575:                        if ( p->td < td )
                   1576:                                goto find;
                   1577:                        else if ( p->td == td ) {
                   1578:                                c = ndl_compare(p->lcm,lcm);
                   1579:                                if ( c < 0 )
                   1580:                                        goto find;
                   1581:                                else if ( c == 0 ) {
                   1582:                                        tlen = nd_length(nd_ps[p->i1])+nd_length(nd_ps[p->i2]);
                   1583:                                        if ( tlen < len )
                   1584:                                                goto find;
                   1585:                                }
                   1586:                        }
                   1587:                }
                   1588:                continue;
                   1589: find:
                   1590:                ml = l;
                   1591:                m = p;
                   1592:                lcm = m->lcm;
                   1593:                s = m->sugar;
                   1594:                td = m->td;
                   1595:                len = tlen;
                   1596:        }
                   1597:        if ( !ml ) *prest = NEXT(m);
                   1598:        else {
                   1599:                NEXT(ml) = NEXT(m);
                   1600:                *prest = d;
                   1601:        }
                   1602:        NEXT(m) = 0;
                   1603:        return m;
                   1604: }
                   1605:
                   1606: int nd_newps(ND a)
                   1607: {
1.3     ! noro     1608:        int len;
        !          1609:
1.1       noro     1610:        if ( nd_psn == nd_pslen ) {
                   1611:                nd_pslen *= 2;
                   1612:                nd_ps = (ND *)REALLOC((char *)nd_ps,nd_pslen*sizeof(ND));
1.3     ! noro     1613:                nd_psv = (NDV *)REALLOC((char *)nd_psv,nd_pslen*sizeof(NDV));
1.1       noro     1614:                nd_bound = (unsigned int **)
                   1615:                        REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
                   1616:        }
                   1617:        nd_monic(a);
                   1618:        nd_ps[nd_psn] = a;
1.3     ! noro     1619:        nd_psv[nd_psn] = ndtondv(a);
1.1       noro     1620:        nd_bound[nd_psn] = nd_compute_bound(a);
1.3     ! noro     1621:        len = nd_psv[nd_psn]->len;
        !          1622:        if ( len > nmv_len ) {
        !          1623:                nmv_len = 2*len;
        !          1624:                BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
        !          1625:        }
1.1       noro     1626:        return nd_psn++;
                   1627: }
                   1628:
                   1629: NODE NODE_sortb(NODE f,int);
                   1630: ND dptond(DP);
                   1631: DP ndtodp(ND);
                   1632:
                   1633: NODE nd_setup(NODE f)
                   1634: {
1.3     ! noro     1635:        int i,td,len;
1.1       noro     1636:        NODE s,s0,f0;
                   1637:
                   1638:        nd_found = 0;
                   1639:        nd_notfirst = 0;
                   1640:        nd_create = 0;
                   1641: #if 0
                   1642:        f0 = f = NODE_sortb(f,1);
                   1643: #endif
                   1644:        nd_psn = length(f); nd_pslen = 2*nd_psn;
                   1645:        nd_ps = (ND *)MALLOC(nd_pslen*sizeof(ND));
1.3     ! noro     1646:        nd_psv = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
1.1       noro     1647:        nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
                   1648:        nd_bpe = 4;
                   1649:        nd_setup_parameters();
                   1650:        nd_free_private_storage();
1.3     ! noro     1651:        len = 0;
1.1       noro     1652:        for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
                   1653:                nd_ps[i] = dptond((DP)BDY(f));
1.3     ! noro     1654:                nd_psv[i] = ndtondv(nd_ps[i]);
1.1       noro     1655:                nd_monic(nd_ps[i]);
                   1656:                nd_bound[i] = nd_compute_bound(nd_ps[i]);
1.3     ! noro     1657:                len = MAX(len,nd_psv[i]->len);
1.1       noro     1658:        }
                   1659:        nd_red = (NM *)MALLOC(REDTAB_LEN*sizeof(NM));
1.3     ! noro     1660:        nmv_len = 16*len;
        !          1661:        NEWNDV(ndv_red);
        !          1662:        BDY(ndv_red) = (NMV)MALLOC_ATOMIC(nmv_len*nmv_adv);
1.1       noro     1663:        for ( s0 = 0, i = 0; i < nd_psn; i++ ) {
                   1664:                NEXTNODE(s0,s); BDY(s) = (pointer)i;
                   1665:        }
                   1666:        if ( s0 ) NEXT(s) = 0;
                   1667:        return s0;
                   1668: }
                   1669:
                   1670: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
                   1671: {
                   1672:        struct order_spec ord1;
                   1673:        VL fv,vv,vc;
                   1674:        NODE fd,fd0,r,r0,t,x,s,xx;
                   1675:        DP a,b,c;
                   1676:
                   1677:        get_vars((Obj)f,&fv); pltovl(v,&vv);
                   1678:        nd_nvar = length(vv);
                   1679:        if ( ord->id )
                   1680:                error("nd_gr : unsupported order");
                   1681:        switch ( ord->ord.simple ) {
                   1682:                case 0:
                   1683:                        is_rlex = 1;
                   1684:                        break;
                   1685:                case 1:
                   1686:                        is_rlex = 0;
                   1687:                        break;
                   1688:                default:
                   1689:                        error("nd_gr : unsupported order");
                   1690:        }
                   1691:        initd(ord);
                   1692:        nd_mod = m;
                   1693:        for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                   1694:                ptod(CO,vv,(P)BDY(t),&b);
                   1695:                _dp_mod(b,m,0,&c);
                   1696:                if ( c ) {
                   1697:                        NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
                   1698:                }
                   1699:        }
                   1700:        if ( fd0 ) NEXT(fd) = 0;
                   1701:        s = nd_setup(fd0);
                   1702:        x = nd_gb(s);
                   1703: #if 0
                   1704:        x = nd_reduceall(x,m);
                   1705: #endif
                   1706:        for ( r0 = 0; x; x = NEXT(x) ) {
                   1707:                NEXTNODE(r0,r);
                   1708:                a = ndtodp(nd_ps[(int)BDY(x)]);
                   1709:                _dtop_mod(CO,vv,a,(P *)&BDY(r));
                   1710:        }
                   1711:        if ( r0 ) NEXT(r) = 0;
                   1712:        MKLIST(*rp,r0);
                   1713:        fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
                   1714:                nd_found,nd_notfirst,nd_create);
                   1715: }
                   1716:
                   1717: void dltondl(int n,DL dl,unsigned int *r)
                   1718: {
                   1719:        unsigned int *d;
                   1720:        int i;
                   1721:
                   1722:        d = dl->d;
                   1723:        bzero(r,nd_wpd*sizeof(unsigned int));
                   1724:        if ( is_rlex )
                   1725:                for ( i = 0; i < n; i++ )
                   1726:                        r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe));
                   1727:        else
                   1728:                for ( i = 0; i < n; i++ )
                   1729:                        r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe);
                   1730: }
                   1731:
                   1732: DL ndltodl(int n,int td,unsigned int *ndl)
                   1733: {
                   1734:        DL dl;
                   1735:        int *d;
                   1736:        int i;
                   1737:
                   1738:        NEWDL(dl,n);
                   1739:        dl->td = td;
                   1740:        d = dl->d;
                   1741:        if ( is_rlex )
                   1742:                for ( i = 0; i < n; i++ )
                   1743:                        d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
                   1744:                                &((1<<nd_bpe)-1);
                   1745:        else
                   1746:                for ( i = 0; i < n; i++ )
                   1747:                        d[i] = (ndl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
                   1748:                                &((1<<nd_bpe)-1);
                   1749:        return dl;
                   1750: }
                   1751:
                   1752: ND dptond(DP p)
                   1753: {
                   1754:        ND d;
                   1755:        NM m0,m;
                   1756:        MP t;
                   1757:        int n;
                   1758:
                   1759:        if ( !p )
                   1760:                return 0;
                   1761:        n = NV(p);
                   1762:        m0 = 0;
                   1763:        for ( t = BDY(p); t; t = NEXT(t) ) {
                   1764:                NEXTNM(m0,m);
                   1765:                m->c = ITOS(t->c);
                   1766:                m->td = t->dl->td;
                   1767:                dltondl(n,t->dl,m->dl);
                   1768:        }
                   1769:        NEXT(m) = 0;
                   1770:        MKND(n,m0,d);
                   1771:        d->nv = n;
                   1772:        d->sugar = p->sugar;
                   1773:        return d;
                   1774: }
                   1775:
                   1776: DP ndtodp(ND p)
                   1777: {
                   1778:        DP d;
                   1779:        MP m0,m;
                   1780:        NM t;
                   1781:        int n;
                   1782:
                   1783:        if ( !p )
                   1784:                return 0;
                   1785:        n = NV(p);
                   1786:        m0 = 0;
                   1787:        for ( t = BDY(p); t; t = NEXT(t) ) {
                   1788:                NEXTMP(m0,m);
                   1789:                m->c = STOI(t->c);
                   1790:                m->dl = ndltodl(n,t->td,t->dl);
                   1791:        }
                   1792:        NEXT(m) = 0;
                   1793:        MKDP(n,m0,d);
                   1794:        d->sugar = p->sugar;
                   1795:        return d;
                   1796: }
                   1797:
                   1798: void ndl_print(unsigned int *dl)
                   1799: {
                   1800:        int n;
                   1801:        int i;
                   1802:
                   1803:        n = nd_nvar;
                   1804:        printf("<<");
                   1805:        if ( is_rlex )
                   1806:                for ( i = 0; i < n; i++ )
                   1807:                        printf(i==n-1?"%d":"%d,",
                   1808:                                (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
                   1809:                                        &((1<<nd_bpe)-1));
                   1810:        else
                   1811:                for ( i = 0; i < n; i++ )
                   1812:                        printf(i==n-1?"%d":"%d,",
                   1813:                                (dl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
                   1814:                                        &((1<<nd_bpe)-1));
                   1815:        printf(">>");
                   1816: }
                   1817:
                   1818: void nd_print(ND p)
                   1819: {
                   1820:        NM m;
                   1821:
                   1822:        if ( !p )
                   1823:                printf("0\n");
                   1824:        else {
                   1825:                for ( m = BDY(p); m; m = NEXT(m) ) {
                   1826:                        printf("+%d*",m->c);
                   1827:                        ndl_print(m->dl);
                   1828:                }
                   1829:                printf("\n");
                   1830:        }
                   1831: }
                   1832:
                   1833: void ndp_print(ND_pairs d)
                   1834: {
                   1835:        ND_pairs t;
                   1836:
                   1837:        for ( t = d; t; t = NEXT(t) ) {
                   1838:                printf("%d,%d ",t->i1,t->i2);
                   1839:        }
                   1840:        printf("\n");
                   1841: }
                   1842:
                   1843: void nd_monic(ND p)
                   1844: {
                   1845:        if ( !p )
                   1846:                return;
                   1847:        else
                   1848:                nd_mul_c(p,invm(HC(p),nd_mod));
                   1849: }
                   1850:
                   1851: void nd_mul_c(ND p,int mul)
                   1852: {
                   1853:        NM m;
                   1854:        int c,c1;
                   1855:
                   1856:        if ( !p )
                   1857:                return;
                   1858:        for ( m = BDY(p); m; m = NEXT(m) ) {
                   1859:                c1 = C(m);
                   1860:                DMAR(c1,mul,0,nd_mod,c);
                   1861:                C(m) = c;
                   1862:        }
                   1863: }
                   1864:
                   1865: void nd_free(ND p)
                   1866: {
                   1867:        NM t,s;
                   1868:
                   1869:        if ( !p )
                   1870:                return;
                   1871:        t = BDY(p);
                   1872:        while ( t ) {
                   1873:                s = NEXT(t);
                   1874:                FREENM(t);
                   1875:                t = s;
                   1876:        }
                   1877:        FREEND(p);
                   1878: }
                   1879:
                   1880: void nd_append_red(unsigned int *d,int td,int i)
                   1881: {
                   1882:        NM m,m0;
                   1883:        int h;
                   1884:
                   1885:        NEWNM(m);
                   1886:        h = ndl_hash_value(td,d);
                   1887:        m->c = i;
                   1888:        m->td = td;
                   1889:        bcopy(d,m->dl,nd_wpd*sizeof(unsigned int));
                   1890:        NEXT(m) = nd_red[h];
                   1891:        nd_red[h] = m;
                   1892: }
                   1893:
                   1894: unsigned int *nd_compute_bound(ND p)
                   1895: {
                   1896:        unsigned int *d1,*d2,*t;
                   1897:        NM m;
                   1898:
                   1899:        if ( !p )
                   1900:                return 0;
                   1901:        d1 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                   1902:        d2 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
                   1903:        bcopy(HDL(p),d1,nd_wpd*sizeof(unsigned int));
                   1904:        for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
                   1905:                ndl_lcm(m->dl,d1,d2);
                   1906:                t = d1; d1 = d2; d2 = t;
                   1907:        }
                   1908:        t = (unsigned int *)MALLOC_ATOMIC(nd_wpd*sizeof(unsigned int));
                   1909:        bcopy(d1,t,nd_wpd*sizeof(unsigned int));
                   1910:        return t;
                   1911: }
                   1912:
                   1913: void nd_setup_parameters() {
                   1914:        int i;
                   1915:
                   1916:        nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
                   1917:        nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
                   1918:        if ( nd_bpe < 32 ) {
                   1919:                nd_mask0 = (1<<nd_bpe)-1;
                   1920:        } else {
                   1921:                nd_mask0 = 0xffffffff;
                   1922:        }
                   1923:        bzero(nd_mask,sizeof(nd_mask));
                   1924:        nd_mask1 = 0;
                   1925:        for ( i = 0; i < nd_epw; i++ ) {
                   1926:                nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
                   1927:                nd_mask1 |= (1<<(nd_bpe-1))<<(i*nd_bpe);
                   1928:        }
1.3     ! noro     1929:        nmv_adv = sizeof(struct oNMV)+(nd_wpd-1)*sizeof(unsigned int);
1.1       noro     1930: }
                   1931:
                   1932: ND_pairs nd_reconstruct(ND_pairs d)
                   1933: {
                   1934:        int i,obpe;
                   1935:        NM prev_nm_free_list;
                   1936:        ND_pairs s0,s,t,prev_ndp_free_list;
                   1937:
                   1938:        obpe = nd_bpe;
                   1939:        switch ( nd_bpe ) {
                   1940:                case 4: nd_bpe = 6; break;
                   1941:                case 6: nd_bpe = 8; break;
                   1942:                case 8: nd_bpe = 16; break;
                   1943:                case 16: nd_bpe = 32; break;
                   1944:        }
                   1945:        nd_setup_parameters();
                   1946:        prev_nm_free_list = _nm_free_list;
                   1947:        prev_ndp_free_list = _ndp_free_list;
                   1948:        _nm_free_list = 0;
                   1949:        _ndp_free_list = 0;
                   1950:        for ( i = 0; i < nd_psn; i++ ) {
                   1951:                nd_ps[i] = nd_dup(nd_ps[i],obpe);
1.3     ! noro     1952:                nd_psv[i] = ndtondv(nd_ps[i]);
1.1       noro     1953:                nd_bound[i] = nd_compute_bound(nd_ps[i]);
                   1954:        }
                   1955:        s0 = 0;
                   1956:        for ( t = d; t; t = NEXT(t) ) {
                   1957:                NEXTND_pairs(s0,s);
                   1958:                s->i1 = t->i1;
                   1959:                s->i2 = t->i2;
                   1960:                s->td = t->td;
                   1961:                s->sugar = t->sugar;
                   1962:                ndl_dup(obpe,t->lcm,s->lcm);
                   1963:        }
                   1964:        if ( s0 ) NEXT(s) = 0;
                   1965:        prev_nm_free_list = 0;
                   1966:        prev_ndp_free_list = 0;
1.3     ! noro     1967:        BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
1.1       noro     1968:        GC_gcollect();
                   1969:        return s0;
                   1970: }
                   1971:
                   1972: void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
                   1973: {
                   1974:        int n,i,ei,oepw,cepw,cbpe;
                   1975:
                   1976:        n = nd_nvar;
                   1977:        oepw = (sizeof(unsigned int)*8)/obpe;
                   1978:        cepw = nd_epw;
                   1979:        cbpe = nd_bpe;
                   1980:        if ( is_rlex )
                   1981:                for ( i = 0; i < n; i++ ) {
                   1982:                        ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe))
                   1983:                                &((1<<obpe)-1);
                   1984:                        r[(n-1-i)/cepw] |= (ei<<((cepw-((n-1-i)%cepw)-1)*cbpe));
                   1985:                }
                   1986:        else
                   1987:                for ( i = 0; i < n; i++ ) {
                   1988:                        ei = (d[i/oepw]>>((oepw-(i%oepw)-1)*obpe))
                   1989:                                &((1<<obpe)-1);
                   1990:                        r[i/cepw] |= (ei<<((cepw-(i%cepw)-1)*cbpe));
                   1991:                }
                   1992: }
                   1993:
                   1994: ND nd_dup(ND p,int obpe)
                   1995: {
                   1996:        NM m,mr,mr0;
                   1997:        int c,n;
                   1998:        ND r;
                   1999:
                   2000:        if ( !p )
                   2001:                return 0;
                   2002:        else {
                   2003:                n = NV(p); m = BDY(p);
                   2004:                for ( mr0 = 0; m; m = NEXT(m) ) {
                   2005:                        NEXTNM(mr0,mr);
                   2006:                        C(mr) = C(m);
                   2007:                        mr->td = m->td;
                   2008:                        ndl_dup(obpe,m->dl,mr->dl);
                   2009:                }
                   2010:                NEXT(mr) = 0;
                   2011:                MKND(NV(p),mr0,r);
                   2012:                r->sugar = p->sugar;
                   2013:                return r;
                   2014:        }
                   2015: }
1.3     ! noro     2016:
        !          2017: NDV ndtondv(ND p)
        !          2018: {
        !          2019:        NDV d;
        !          2020:        NMV m,m0;
        !          2021:        NM t;
        !          2022:        int i,len;
        !          2023:
        !          2024:        if ( !p )
        !          2025:                return 0;
        !          2026:        len = nd_length(p);
        !          2027:        m0 = m = (NMV)MALLOC_ATOMIC(len*nmv_adv);
        !          2028:        for ( t = BDY(p), i = 0; t; t = NEXT(t), i++, NMV_ADV(m) ) {
        !          2029:                m->td = t->td;
        !          2030:                bcopy(t->dl,m->dl,nd_wpd*sizeof(unsigned int));
        !          2031:                m->c = t->c;
        !          2032:        }
        !          2033:        MKNDV(NV(p),m0,len,d);
        !          2034:        d->sugar = p->sugar;
        !          2035:        return d;
        !          2036: }
        !          2037:
        !          2038: void ndv_print(NDV p)
        !          2039: {
        !          2040:        NMV m;
        !          2041:        int i,len;
        !          2042:
        !          2043:        if ( !p )
        !          2044:                printf("0\n");
        !          2045:        else {
        !          2046:                len = p->len;
        !          2047:                for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
        !          2048:                        printf("+%d*",m->c);
        !          2049:                        ndl_print(m->dl);
        !          2050:                }
        !          2051:                printf("\n");
        !          2052:        }
        !          2053: }
        !          2054:

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