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