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