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