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

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

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