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

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

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