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