![[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.12 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.11 2003/07/26 01:00:32 noro Exp $ */
1.2 noro 2:
1.1 noro 3: #include "ca.h"
4: #include "inline.h"
5:
1.6 noro 6: #define USE_NDV 1
1.4 noro 7:
1.1 noro 8: #if defined(__GNUC__)
9: #define INLINE inline
10: #elif defined(VISUAL)
11: #define INLINE __inline
12: #else
13: #define INLINE
14: #endif
15:
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;
26: int sugar;
27: } *ND;
28:
1.3 noro 29: typedef struct oNDV {
30: struct oNMV *body;
31: int nv;
32: int sugar;
33: int len;
34: } *NDV;
35:
1.1 noro 36: typedef struct oNM {
37: struct oNM *next;
38: int td;
39: int c;
40: unsigned int dl[1];
41: } *NM;
42:
1.3 noro 43: typedef struct oNMV {
44: int td;
45: int c;
46: unsigned int dl[1];
47: } *NMV;
48:
1.1 noro 49: typedef struct oND_pairs {
50: struct oND_pairs *next;
51: int i1,i2;
52: int td,sugar;
53: unsigned int lcm[1];
54: } *ND_pairs;
55:
56: static unsigned int **nd_bound;
57: int nd_mod,nd_nvar;
58: int is_rlex;
59: int nd_epw,nd_bpe,nd_wpd;
60: unsigned int nd_mask[32];
61: unsigned int nd_mask0,nd_mask1;
62: NM _nm_free_list;
63: ND _nd_free_list;
64: ND_pairs _ndp_free_list;
65: NM *nd_red;
66: int nd_red_len;
67:
68: extern int Top,Reverse;
69: int nd_psn,nd_pslen;
70: int nd_found,nd_create,nd_notfirst;
1.6 noro 71: int *nd_psl;
1.11 noro 72: NM *nd_psh;
1.1 noro 73:
74: void GC_gcollect();
75: NODE append_one(NODE,int);
76:
77: #define HTD(d) ((d)->body->td)
78: #define HDL(d) ((d)->body->dl)
79: #define HC(d) ((d)->body->c)
80:
81: #define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list)
82: #define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list)
83: #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)
84:
85: #define NEXTNM(r,c) \
86: if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);}
87: #define NEXTNM2(r,c,s) \
88: if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
89: #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)
90: #define FREENDP(m) NEXT(m)=_ndp_free_list; _ndp_free_list=(m)
91: #define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m)
92:
93: #define NEXTND_pairs(r,c) \
94: if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
95:
96: ND_pairs crit_B( ND_pairs d, int s );
97: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp);
98: NODE nd_setup(NODE f);
99: int nd_newps(ND a);
100: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest );
101: NODE update_base(NODE nd,int ndp);
102: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest );
103: int crit_2( int dp1, int dp2 );
104: ND_pairs crit_F( ND_pairs d1 );
105: ND_pairs crit_M( ND_pairs d1 );
106: ND_pairs nd_newpairs( NODE g, int t );
107: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t);
108: NODE nd_gb(NODE f);
109: void nd_free_private_storage();
110: void _NM_alloc();
111: void _ND_alloc();
112: int ndl_td(unsigned int *d);
113: ND nd_add(ND p1,ND p2);
114: ND nd_mul_nm(ND p,NM m0);
1.6 noro 115: ND nd_mul_ind_nm(int index,NM m0);
1.1 noro 116: ND nd_mul_term(ND p,int td,unsigned int *d);
117: int nd_sp(ND_pairs p,ND *nf);
1.6 noro 118: int nd_find_reducer(ND g);
1.1 noro 119: int nd_nf(ND g,int full,ND *nf);
120: ND nd_reduce(ND p1,ND p2);
121: ND nd_reduce_special(ND p1,ND p2);
122: void nd_free(ND p);
123: void ndl_print(unsigned int *dl);
124: void nd_print(ND p);
125: void ndp_print(ND_pairs d);
126: int nd_length(ND p);
127: void nd_monic(ND p);
128: void nd_mul_c(ND p,int mul);
129: void nd_free_redlist();
130: void nd_append_red(unsigned int *d,int td,int i);
131: unsigned int *nd_compute_bound(ND p);
1.5 noro 132: unsigned int *dp_compute_bound(DP p);
1.1 noro 133: ND_pairs nd_reconstruct(ND_pairs);
134: void nd_setup_parameters();
1.11 noro 135: void nd_realloc(ND p,int obpe);
1.6 noro 136: ND nd_copy(ND p);
1.1 noro 137: void ndl_dup(int obpe,unsigned int *d,unsigned int *r);
1.4 noro 138:
139: #if USE_NDV
1.11 noro 140: static NDV *nd_ps;
141:
1.4 noro 142: #define NMV_ADV(m) (m = (NMV)(((char *)m)+nmv_adv))
143: #define NEWNDV(d) ((d)=(NDV)MALLOC(sizeof(struct oNDV)))
144: #define MKNDV(n,m,l,d) NEWNDV(d); (d)->nv=(n); BDY(d)=(m); (d)->len = l;
145:
146: int nmv_adv;
147: int nmv_len;
148: NDV ndv_red;
149:
1.11 noro 150: void ndv_mul_c(NDV p,int mul);
1.4 noro 151: ND ndv_add(ND p1,NDV p2);
1.3 noro 152: NDV ndtondv(ND p);
153: void ndv_mul_nm(NDV pv,NM m,NDV r);
1.9 noro 154: ND ndv_mul_nm_create(NDV p,NM m0);
1.11 noro 155: void ndv_realloc(NDV p,int obpe,int oadv);
156: NDV dptondv(DP);
157: DP ndvtodp(NDV);
158: #else
159: static ND *nd_ps;
160: ND dptond(DP);
161: DP ndtodp(ND);
1.4 noro 162: #endif
1.1 noro 163:
164: void nd_free_private_storage()
165: {
166: _nd_free_list = 0;
167: _nm_free_list = 0;
1.5 noro 168: _ndp_free_list = 0;
1.11 noro 169: bzero(nd_red,sizeof(REDTAB_LEN*sizeof(NM)));
1.1 noro 170: GC_gcollect();
171: }
172:
173: void _NM_alloc()
174: {
175: NM p;
176: int i;
177:
1.11 noro 178: for ( i = 0; i < 1024; i++ ) {
1.1 noro 179: p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
180: p->next = _nm_free_list; _nm_free_list = p;
181: }
182: }
183:
184: void _ND_alloc()
185: {
186: ND p;
187: int i;
188:
189: for ( i = 0; i < 1024; i++ ) {
190: p = (ND)GC_malloc(sizeof(struct oND));
191: p->body = (NM)_nd_free_list; _nd_free_list = p;
192: }
193: }
194:
195: void _NDP_alloc()
196: {
197: ND_pairs p;
198: int i;
199:
1.11 noro 200: for ( i = 0; i < 1024; i++ ) {
1.1 noro 201: p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs)
202: +(nd_wpd-1)*sizeof(unsigned int));
203: p->next = _ndp_free_list; _ndp_free_list = p;
204: }
205: }
206:
207: INLINE nd_length(ND p)
208: {
209: NM m;
210: int i;
211:
212: if ( !p )
213: return 0;
214: else {
215: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
216: return i;
217: }
218: }
219:
220: INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2)
221: {
222: unsigned int u1,u2;
223: int i,j;
224:
225: switch ( nd_bpe ) {
226: case 4:
227: for ( i = 0; i < nd_wpd; i++ ) {
228: u1 = d1[i]; u2 = d2[i];
229: if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0;
230: if ( (u1&0xf000000) < (u2&0xf000000) ) return 0;
231: if ( (u1&0xf00000) < (u2&0xf00000) ) return 0;
232: if ( (u1&0xf0000) < (u2&0xf0000) ) return 0;
233: if ( (u1&0xf000) < (u2&0xf000) ) return 0;
234: if ( (u1&0xf00) < (u2&0xf00) ) return 0;
235: if ( (u1&0xf0) < (u2&0xf0) ) return 0;
236: if ( (u1&0xf) < (u2&0xf) ) return 0;
237: }
238: return 1;
239: break;
240: case 6:
241: for ( i = 0; i < nd_wpd; i++ ) {
242: u1 = d1[i]; u2 = d2[i];
243: if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0;
244: if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0;
245: if ( (u1&0x3f000) < (u2&0x3f000) ) return 0;
246: if ( (u1&0xfc0) < (u2&0xfc0) ) return 0;
247: if ( (u1&0x3f) < (u2&0x3f) ) return 0;
248: }
249: return 1;
250: break;
251: case 8:
252: for ( i = 0; i < nd_wpd; i++ ) {
253: u1 = d1[i]; u2 = d2[i];
254: if ( (u1&0xff000000) < (u2&0xff000000) ) return 0;
255: if ( (u1&0xff0000) < (u2&0xff0000) ) return 0;
256: if ( (u1&0xff00) < (u2&0xff00) ) return 0;
257: if ( (u1&0xff) < (u2&0xff) ) return 0;
258: }
259: return 1;
260: break;
261: case 16:
262: for ( i = 0; i < nd_wpd; i++ ) {
263: u1 = d1[i]; u2 = d2[i];
264: if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0;
265: if ( (u1&0xffff) < (u2&0xffff) ) return 0;
266: }
267: return 1;
268: break;
269: case 32:
270: for ( i = 0; i < nd_wpd; i++ )
271: if ( d1[i] < d2[i] ) return 0;
272: return 1;
273: break;
274: default:
275: for ( i = 0; i < nd_wpd; i++ ) {
276: u1 = d1[i]; u2 = d2[i];
277: for ( j = 0; j < nd_epw; j++ )
278: if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0;
279: }
280: return 1;
281: }
282: }
283:
284: void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
285: {
286: unsigned int t1,t2,u,u1,u2;
287: int i,j;
288:
289: switch ( nd_bpe ) {
290: case 4:
291: for ( i = 0; i < nd_wpd; i++ ) {
292: u1 = d1[i]; u2 = d2[i];
293: t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2;
294: t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2;
295: t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2;
296: t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2;
297: t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2;
298: t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2;
299: t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2;
300: t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2;
301: d[i] = u;
302: }
303: break;
304: case 6:
305: for ( i = 0; i < nd_wpd; i++ ) {
306: u1 = d1[i]; u2 = d2[i];
307: t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2;
308: t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2;
309: t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2;
310: t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2;
311: t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2;
312: d[i] = u;
313: }
314: break;
315: case 8:
316: for ( i = 0; i < nd_wpd; i++ ) {
317: u1 = d1[i]; u2 = d2[i];
318: t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2;
319: t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2;
320: t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2;
321: t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2;
322: d[i] = u;
323: }
324: break;
325: case 16:
326: for ( i = 0; i < nd_wpd; i++ ) {
327: u1 = d1[i]; u2 = d2[i];
328: t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2;
329: t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2;
330: d[i] = u;
331: }
332: break;
333: case 32:
334: for ( i = 0; i < nd_wpd; i++ ) {
335: u1 = d1[i]; u2 = d2[i];
336: d[i] = u1>u2?u1:u2;
337: }
338: break;
339: default:
340: for ( i = 0; i < nd_wpd; i++ ) {
341: u1 = d1[i]; u2 = d2[i];
342: for ( j = 0, u = 0; j < nd_epw; j++ ) {
343: t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2;
344: }
345: d[i] = u;
346: }
347: break;
348: }
349: }
350:
351: int ndl_td(unsigned int *d)
352: {
353: unsigned int t,u;
354: int i,j;
355:
356: for ( t = 0, i = 0; i < nd_wpd; i++ ) {
357: u = d[i];
358: for ( j = 0; j < nd_epw; j++, u>>=nd_bpe )
359: t += (u&nd_mask0);
360: }
361: return t;
362: }
363:
364: INLINE int ndl_compare(unsigned int *d1,unsigned int *d2)
365: {
366: int i;
367:
368: for ( i = 0; i < nd_wpd; i++, d1++, d2++ )
369: if ( *d1 > *d2 )
370: return is_rlex ? -1 : 1;
371: else if ( *d1 < *d2 )
372: return is_rlex ? 1 : -1;
373: return 0;
374: }
375:
376: INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)
377: {
378: int i;
379:
380: for ( i = 0; i < nd_wpd; i++ )
381: if ( d1[i] != d2[i] )
382: return 0;
383: return 1;
384: }
385:
1.6 noro 386: INLINE void ndl_copy(unsigned int *d1,unsigned int *d2)
387: {
388: int i;
389:
390: switch ( nd_wpd ) {
391: case 1:
392: d2[0] = d1[0];
393: break;
394: case 2:
395: d2[0] = d1[0];
396: d2[1] = d1[1];
397: break;
398: default:
399: for ( i = 0; i < nd_wpd; i++ )
400: d2[i] = d1[i];
401: break;
402: }
403: }
404:
1.1 noro 405: INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d)
406: {
407: int i;
408:
1.6 noro 409: switch ( nd_wpd ) {
410: case 1:
411: d[0] = d1[0]+d2[0];
412: break;
413: case 2:
414: d[0] = d1[0]+d2[0];
415: d[1] = d1[1]+d2[1];
416: break;
417: default:
418: for ( i = 0; i < nd_wpd; i++ )
419: d[i] = d1[i]+d2[i];
420: break;
421: }
422: }
423:
424: INLINE void ndl_add2(unsigned int *d1,unsigned int *d2)
425: {
426: int i;
427:
428: switch ( nd_wpd ) {
429: case 1:
430: d2[0] += d1[0];
431: break;
432: case 2:
433: d2[0] += d1[0];
434: d2[1] += d1[1];
435: break;
436: default:
437: for ( i = 0; i < nd_wpd; i++ )
438: d2[i] += d1[i];
439: break;
1.1 noro 440: }
441: }
442:
443: void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
444: {
445: int i;
446:
447: for ( i = 0; i < nd_wpd; i++ )
448: d[i] = d1[i]-d2[i];
449: }
450:
451: int ndl_disjoint(unsigned int *d1,unsigned int *d2)
452: {
453: unsigned int t1,t2,u,u1,u2;
454: int i,j;
455:
456: switch ( nd_bpe ) {
457: case 4:
458: for ( i = 0; i < nd_wpd; i++ ) {
459: u1 = d1[i]; u2 = d2[i];
460: t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0;
461: t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0;
462: t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0;
463: t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0;
464: t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0;
465: t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0;
466: t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0;
467: t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0;
468: }
469: return 1;
470: break;
471: case 6:
472: for ( i = 0; i < nd_wpd; i++ ) {
473: u1 = d1[i]; u2 = d2[i];
474: t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0;
475: t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0;
476: t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0;
477: t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0;
478: t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0;
479: }
480: return 1;
481: break;
482: case 8:
483: for ( i = 0; i < nd_wpd; i++ ) {
484: u1 = d1[i]; u2 = d2[i];
485: t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0;
486: t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0;
487: t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0;
488: t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0;
489: }
490: return 1;
491: break;
492: case 16:
493: for ( i = 0; i < nd_wpd; i++ ) {
494: u1 = d1[i]; u2 = d2[i];
495: t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0;
496: t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0;
497: }
498: return 1;
499: break;
500: case 32:
501: for ( i = 0; i < nd_wpd; i++ )
502: if ( d1[i] && d2[i] ) return 0;
503: return 1;
504: break;
505: default:
506: for ( i = 0; i < nd_wpd; i++ ) {
507: u1 = d1[i]; u2 = d2[i];
508: for ( j = 0; j < nd_epw; j++ ) {
509: if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0;
510: u1 >>= nd_bpe; u2 >>= nd_bpe;
511: }
512: }
513: return 1;
514: break;
515: }
516: }
517:
518: ND nd_reduce(ND p1,ND p2)
519: {
520: int c,c1,c2,t,td,td2,mul;
521: NM m2,prev,head,cur,new;
522: unsigned int *d;
523:
524: if ( !p1 )
525: return 0;
526: else {
527: c2 = invm(HC(p2),nd_mod);
528: c1 = nd_mod-HC(p1);
529: DMAR(c1,c2,0,nd_mod,mul);
530: td = HTD(p1)-HTD(p2);
531: d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
532: ndl_sub(HDL(p1),HDL(p2),d);
533: prev = 0; head = cur = BDY(p1);
534: NEWNM(new);
535: for ( m2 = BDY(p2); m2; ) {
536: td2 = new->td = m2->td+td;
537: ndl_add(m2->dl,d,new->dl);
538: if ( !cur ) {
539: c1 = C(m2);
540: DMAR(c1,mul,0,nd_mod,c2);
541: C(new) = c2;
542: if ( !prev ) {
543: prev = new;
544: NEXT(prev) = 0;
545: head = prev;
546: } else {
547: NEXT(prev) = new;
548: NEXT(new) = 0;
549: prev = new;
550: }
551: m2 = NEXT(m2);
552: NEWNM(new);
553: continue;
554: }
555: if ( cur->td > td2 )
556: c = 1;
557: else if ( cur->td < td2 )
558: c = -1;
559: else
560: c = ndl_compare(cur->dl,new->dl);
561: switch ( c ) {
562: case 0:
563: c2 = C(m2);
564: c1 = C(cur);
565: DMAR(c2,mul,c1,nd_mod,t);
566: if ( t )
567: C(cur) = t;
568: else if ( !prev ) {
569: head = NEXT(cur);
570: FREENM(cur);
571: cur = head;
572: } else {
573: NEXT(prev) = NEXT(cur);
574: FREENM(cur);
575: cur = NEXT(prev);
576: }
577: m2 = NEXT(m2);
578: break;
579: case 1:
580: prev = cur;
581: cur = NEXT(cur);
582: break;
583: case -1:
584: if ( !prev ) {
585: /* cur = head */
586: prev = new;
587: c2 = C(m2);
588: DMAR(c2,mul,0,nd_mod,c1);
589: C(prev) = c1;
590: NEXT(prev) = head;
591: head = prev;
592: } else {
593: c2 = C(m2);
594: DMAR(c2,mul,0,nd_mod,c1);
595: C(new) = c1;
596: NEXT(prev) = new;
597: NEXT(new) = cur;
598: prev = new;
599: }
600: NEWNM(new);
601: m2 = NEXT(m2);
602: break;
603: }
604: }
605: FREENM(new);
606: if ( head ) {
607: BDY(p1) = head;
608: p1->sugar = MAX(p1->sugar,p2->sugar+td);
609: return p1;
610: } else {
611: FREEND(p1);
612: return 0;
613: }
614:
615: }
616: }
617:
618: /* HDL(p1) = HDL(p2) */
619:
620: ND nd_reduce_special(ND p1,ND p2)
621: {
622: int c,c1,c2,t,td,td2,mul;
623: NM m2,prev,head,cur,new;
624:
625: if ( !p1 )
626: return 0;
627: else {
628: c2 = invm(HC(p2),nd_mod);
629: c1 = nd_mod-HC(p1);
630: DMAR(c1,c2,0,nd_mod,mul);
631: prev = 0; head = cur = BDY(p1);
632: NEWNM(new);
633: for ( m2 = BDY(p2); m2; ) {
634: td2 = new->td = m2->td;
635: if ( !cur ) {
636: c1 = C(m2);
637: DMAR(c1,mul,0,nd_mod,c2);
638: C(new) = c2;
639: bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
640: if ( !prev ) {
641: prev = new;
642: NEXT(prev) = 0;
643: head = prev;
644: } else {
645: NEXT(prev) = new;
646: NEXT(new) = 0;
647: prev = new;
648: }
649: m2 = NEXT(m2);
650: NEWNM(new);
651: continue;
652: }
653: if ( cur->td > td2 )
654: c = 1;
655: else if ( cur->td < td2 )
656: c = -1;
657: else
658: c = ndl_compare(cur->dl,m2->dl);
659: switch ( c ) {
660: case 0:
661: c2 = C(m2);
662: c1 = C(cur);
663: DMAR(c2,mul,c1,nd_mod,t);
664: if ( t )
665: C(cur) = t;
666: else if ( !prev ) {
667: head = NEXT(cur);
668: FREENM(cur);
669: cur = head;
670: } else {
671: NEXT(prev) = NEXT(cur);
672: FREENM(cur);
673: cur = NEXT(prev);
674: }
675: m2 = NEXT(m2);
676: break;
677: case 1:
678: prev = cur;
679: cur = NEXT(cur);
680: break;
681: case -1:
682: bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
683: if ( !prev ) {
684: /* cur = head */
685: prev = new;
686: c2 = C(m2);
687: DMAR(c2,mul,0,nd_mod,c1);
688: C(prev) = c1;
689: NEXT(prev) = head;
690: head = prev;
691: } else {
692: c2 = C(m2);
693: DMAR(c2,mul,0,nd_mod,c1);
694: C(new) = c1;
695: NEXT(prev) = new;
696: NEXT(new) = cur;
697: prev = new;
698: }
699: NEWNM(new);
700: m2 = NEXT(m2);
701: break;
702: }
703: }
704: FREENM(new);
705: if ( head ) {
706: BDY(p1) = head;
707: p1->sugar = MAX(p1->sugar,p2->sugar+td);
708: return p1;
709: } else {
710: FREEND(p1);
711: return 0;
712: }
713:
714: }
715: }
716:
1.5 noro 717: int ndl_check_bound2(int index,unsigned int *d2)
1.1 noro 718: {
1.5 noro 719: unsigned int u2;
720: unsigned int *d1;
721: int i,j,ind,k;
1.1 noro 722:
1.5 noro 723: d1 = nd_bound[index];
724: ind = 0;
725: switch ( nd_bpe ) {
726: case 4:
727: for ( i = 0; i < nd_wpd; i++ ) {
728: u2 = d2[i];
729: if ( d1[ind++]+((u2>>28)&0xf) >= 0x10 ) return 1;
730: if ( d1[ind++]+((u2>>24)&0xf) >= 0x10 ) return 1;
731: if ( d1[ind++]+((u2>>20)&0xf) >= 0x10 ) return 1;
732: if ( d1[ind++]+((u2>>16)&0xf) >= 0x10 ) return 1;
733: if ( d1[ind++]+((u2>>12)&0xf) >= 0x10 ) return 1;
734: if ( d1[ind++]+((u2>>8)&0xf) >= 0x10 ) return 1;
735: if ( d1[ind++]+((u2>>4)&0xf) >= 0x10 ) return 1;
736: if ( d1[ind++]+(u2&0xf) >= 0x10 ) return 1;
737: }
738: return 0;
739: break;
740: case 6:
741: for ( i = 0; i < nd_wpd; i++ ) {
742: u2 = d2[i];
743: if ( d1[ind++]+((u2>>24)&0x3f) >= 0x40 ) return 1;
744: if ( d1[ind++]+((u2>>18)&0x3f) >= 0x40 ) return 1;
745: if ( d1[ind++]+((u2>>12)&0x3f) >= 0x40 ) return 1;
746: if ( d1[ind++]+((u2>>6)&0x3f) >= 0x40 ) return 1;
747: if ( d1[ind++]+(u2&0x3f) >= 0x40 ) return 1;
748: }
749: return 0;
750: break;
751: case 8:
752: for ( i = 0; i < nd_wpd; i++ ) {
753: u2 = d2[i];
754: if ( d1[ind++]+((u2>>24)&0xff) >= 0x100 ) return 1;
755: if ( d1[ind++]+((u2>>16)&0xff) >= 0x100 ) return 1;
756: if ( d1[ind++]+((u2>>8)&0xff) >= 0x100 ) return 1;
757: if ( d1[ind++]+(u2&0xff) >= 0x100 ) return 1;
758: }
759: return 0;
760: break;
761: case 16:
762: for ( i = 0; i < nd_wpd; i++ ) {
763: u2 = d2[i];
764: if ( d1[ind++]+((u2>>16)&0xffff) > 0x10000 ) return 1;
765: if ( d1[ind++]+(u2&0xffff) > 0x10000 ) return 1;
766: }
767: return 0;
768: break;
769: case 32:
770: for ( i = 0; i < nd_wpd; i++ )
771: if ( d1[i]+d2[i]<d1[i] ) return 1;
772: return 0;
773: break;
774: default:
775: for ( i = 0; i < nd_wpd; i++ ) {
776: u2 = d2[i];
777: k = (nd_epw-1)*nd_bpe;
778: for ( j = 0; j < nd_epw; j++, k -= nd_bpe )
779: if ( d1[ind++]+((u2>>k)&nd_mask0) > nd_mask0 ) return 1;
780: }
781: return 0;
782: break;
783: }
1.1 noro 784: }
785:
1.6 noro 786: INLINE int ndl_hash_value(int td,unsigned int *d)
1.1 noro 787: {
788: int i;
789: int r;
790:
791: r = td;
792: for ( i = 0; i < nd_wpd; i++ )
793: r = ((r<<16)+d[i])%REDTAB_LEN;
794: return r;
795: }
796:
1.9 noro 797: INLINE int nd_find_reducer(ND g)
1.1 noro 798: {
799: NM m;
1.6 noro 800: int d,k,i;
1.1 noro 801:
802: d = ndl_hash_value(HTD(g),HDL(g));
803: for ( m = nd_red[d], k = 0; m; m = NEXT(m), k++ ) {
804: if ( HTD(g) == m->td && ndl_equal(HDL(g),m->dl) ) {
805: if ( k > 0 ) nd_notfirst++;
806: nd_found++;
1.6 noro 807: return m->c;
1.1 noro 808: }
809: }
810:
811: for ( i = 0; i < nd_psn; i++ ) {
1.11 noro 812: m = nd_psh[i];
813: if ( HTD(g) >= m->td && ndl_reducible(HDL(g),m->dl) ) {
1.1 noro 814: nd_create++;
1.6 noro 815: nd_append_red(HDL(g),HTD(g),i);
816: return i;
1.1 noro 817: }
818: }
1.6 noro 819: return -1;
1.1 noro 820: }
821:
822: ND nd_add(ND p1,ND p2)
823: {
824: int n,c;
825: int t;
826: ND r;
827: NM m1,m2,mr0,mr,s;
828:
829: if ( !p1 )
830: return p2;
831: else if ( !p2 )
832: return p1;
833: else {
834: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
835: if ( m1->td > m2->td )
836: c = 1;
837: else if ( m1->td < m2->td )
838: c = -1;
839: else
840: c = ndl_compare(m1->dl,m2->dl);
841: switch ( c ) {
842: case 0:
843: t = ((C(m1))+(C(m2))) - nd_mod;
844: if ( t < 0 )
845: t += nd_mod;
846: s = m1; m1 = NEXT(m1);
847: if ( t ) {
848: NEXTNM2(mr0,mr,s); C(mr) = (t);
849: } else {
850: FREENM(s);
851: }
852: s = m2; m2 = NEXT(m2); FREENM(s);
853: break;
854: case 1:
855: s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
856: break;
857: case -1:
858: s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
859: break;
860: }
861: }
862: if ( !mr0 )
863: if ( m1 )
864: mr0 = m1;
865: else if ( m2 )
866: mr0 = m2;
867: else
868: return 0;
869: else if ( m1 )
870: NEXT(mr) = m1;
871: else if ( m2 )
872: NEXT(mr) = m2;
873: else
874: NEXT(mr) = 0;
875: BDY(p1) = mr0;
876: p1->sugar = MAX(p1->sugar,p2->sugar);
877: FREEND(p2);
878: return p1;
879: }
880: }
881:
882: #if 1
883: /* ret=1 : success, ret=0 : overflow */
884: int nd_nf(ND g,int full,ND *rp)
885: {
1.11 noro 886: ND d;
1.1 noro 887: NM m,mrd,tail;
1.7 noro 888: NM mul;
1.10 noro 889: int n,sugar,psugar,sugar0,stat,index;
1.6 noro 890: int c,c1,c2;
1.4 noro 891: #if USE_NDV
1.11 noro 892: NDV p,red;
1.4 noro 893: #else
1.11 noro 894: ND p,red;
1.4 noro 895: #endif
1.1 noro 896:
897: if ( !g ) {
898: *rp = 0;
899: return 1;
900: }
1.10 noro 901: sugar0 = sugar = g->sugar;
1.1 noro 902: n = NV(g);
1.7 noro 903: mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
1.1 noro 904: for ( d = 0; g; ) {
1.6 noro 905: index = nd_find_reducer(g);
906: if ( index >= 0 ) {
907: p = nd_ps[index];
1.7 noro 908: ndl_sub(HDL(g),HDL(p),mul->dl);
909: mul->td = HTD(g)-HTD(p);
1.10 noro 910: #if 0
911: if ( d && (p->sugar+mul->td) > sugar ) {
912: goto afo;
913: }
914: #endif
1.7 noro 915: if ( ndl_check_bound2(index,mul->dl) ) {
1.6 noro 916: nd_free(g); nd_free(d);
917: return 0;
918: }
919: c1 = invm(HC(p),nd_mod); c2 = nd_mod-HC(g);
1.7 noro 920: DMAR(c1,c2,0,nd_mod,c); C(mul) = c;
1.4 noro 921: #if USE_NDV
1.11 noro 922: ndv_mul_nm(nd_ps[index],mul,ndv_red);
1.4 noro 923: g = ndv_add(g,ndv_red);
924: sugar = MAX(sugar,ndv_red->sugar);
925: #else
1.7 noro 926: red = nd_mul_ind_nm(index,mul);
1.1 noro 927: g = nd_add(g,red);
1.4 noro 928: sugar = MAX(sugar,red->sugar);
1.1 noro 929: #endif
930: } else if ( !full ) {
931: *rp = g;
932: return 1;
933: } else {
1.10 noro 934: afo:
1.1 noro 935: m = BDY(g);
936: if ( NEXT(m) ) {
937: BDY(g) = NEXT(m); NEXT(m) = 0;
938: } else {
939: FREEND(g); g = 0;
940: }
941: if ( d ) {
942: NEXT(tail)=m;
943: tail=m;
944: } else {
945: MKND(n,m,d);
946: tail = BDY(d);
947: }
948: }
949: }
950: if ( d )
951: d->sugar = sugar;
952: *rp = d;
953: return 1;
954: }
955: #else
956:
957: ND nd_remove_head(ND p)
958: {
959: NM m;
960:
961: m = BDY(p);
962: if ( !NEXT(m) ) {
963: FREEND(p);
964: p = 0;
965: } else
966: BDY(p) = NEXT(m);
967: FREENM(m);
968: return p;
969: }
970:
971: PGeoBucket create_pbucket()
972: {
973: PGeoBucket g;
974:
975: g = CALLOC(1,sizeof(struct oPGeoBucket));
976: g->m = -1;
977: return g;
978: }
979:
980: void add_pbucket(PGeoBucket g,ND d)
981: {
982: int l,k,m;
983:
984: l = nd_length(d);
985: for ( k = 0, m = 1; l > m; k++, m <<= 2 );
986: /* 4^(k-1) < l <= 4^k */
987: d = nd_add(g->body[k],d);
988: for ( ; d && nd_length(d) > 1<<(2*k); k++ ) {
989: g->body[k] = 0;
990: d = nd_add(g->body[k+1],d);
991: }
992: g->body[k] = d;
993: g->m = MAX(g->m,k);
994: }
995:
996: int head_pbucket(PGeoBucket g)
997: {
998: int j,i,c,k,nv,sum;
999: unsigned int *di,*dj;
1000: ND gi,gj;
1001:
1002: k = g->m;
1003: while ( 1 ) {
1004: j = -1;
1005: for ( i = 0; i <= k; i++ ) {
1006: if ( !(gi = g->body[i]) )
1007: continue;
1008: if ( j < 0 ) {
1009: j = i;
1010: gj = g->body[j];
1011: dj = HDL(gj);
1012: sum = HC(gj);
1013: } else {
1014: di = HDL(gi);
1015: nv = NV(gi);
1016: if ( HTD(gi) > HTD(gj) )
1017: c = 1;
1018: else if ( HTD(gi) < HTD(gj) )
1019: c = -1;
1020: else
1021: c = ndl_compare(di,dj);
1022: if ( c > 0 ) {
1023: if ( sum )
1024: HC(gj) = sum;
1025: else
1026: g->body[j] = nd_remove_head(gj);
1027: j = i;
1028: gj = g->body[j];
1029: dj = HDL(gj);
1030: sum = HC(gj);
1031: } else if ( c == 0 ) {
1032: sum = sum+HC(gi)-nd_mod;
1033: if ( sum < 0 )
1034: sum += nd_mod;
1035: g->body[i] = nd_remove_head(gi);
1036: }
1037: }
1038: }
1039: if ( j < 0 )
1040: return -1;
1041: else if ( sum ) {
1042: HC(gj) = sum;
1043: return j;
1044: } else
1045: g->body[j] = nd_remove_head(gj);
1046: }
1047: }
1048:
1049: ND normalize_pbucket(PGeoBucket g)
1050: {
1051: int i;
1052: ND r,t;
1053:
1054: r = 0;
1055: for ( i = 0; i <= g->m; i++ )
1056: r = nd_add(r,g->body[i]);
1057: return r;
1058: }
1059:
1060: ND nd_nf(ND g,int full)
1061: {
1062: ND u,p,d,red;
1063: NODE l;
1064: NM m,mrd;
1065: int sugar,psugar,n,h_reducible,h;
1066: PGeoBucket bucket;
1067:
1068: if ( !g ) {
1069: return 0;
1070: }
1071: sugar = g->sugar;
1072: n = g->nv;
1073: bucket = create_pbucket();
1074: add_pbucket(bucket,g);
1075: d = 0;
1076: while ( 1 ) {
1077: h = head_pbucket(bucket);
1078: if ( h < 0 ) {
1079: if ( d )
1080: d->sugar = sugar;
1081: return d;
1082: }
1083: g = bucket->body[h];
1084: red = nd_find_reducer(g);
1085: if ( red ) {
1086: bucket->body[h] = nd_remove_head(g);
1087: red = nd_remove_head(red);
1088: add_pbucket(bucket,red);
1089: sugar = MAX(sugar,red->sugar);
1090: } else if ( !full ) {
1091: g = normalize_pbucket(bucket);
1092: if ( g )
1093: g->sugar = sugar;
1094: return g;
1095: } else {
1096: m = BDY(g);
1097: if ( NEXT(m) ) {
1098: BDY(g) = NEXT(m); NEXT(m) = 0;
1099: } else {
1100: FREEND(g); g = 0;
1101: }
1102: bucket->body[h] = g;
1103: NEXT(m) = 0;
1104: if ( d ) {
1105: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1106: NEXT(mrd) = m;
1107: } else {
1108: MKND(n,m,d);
1109: }
1110: }
1111: }
1112: }
1113: #endif
1114:
1115: NODE nd_gb(NODE f)
1116: {
1117: int i,nh,sugar,stat;
1118: NODE r,g,gall;
1119: ND_pairs d;
1120: ND_pairs l;
1121: ND h,nf;
1122:
1123: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1124: i = (int)BDY(r);
1125: d = update_pairs(d,g,i);
1126: g = update_base(g,i);
1127: gall = append_one(gall,i);
1128: }
1129: sugar = 0;
1130: while ( d ) {
1131: again:
1132: l = nd_minp(d,&d);
1133: if ( l->sugar != sugar ) {
1134: sugar = l->sugar;
1135: fprintf(asir_out,"%d",sugar);
1136: }
1137: stat = nd_sp(l,&h);
1138: if ( !stat ) {
1139: NEXT(l) = d; d = l;
1140: d = nd_reconstruct(d);
1141: goto again;
1142: }
1143: stat = nd_nf(h,!Top,&nf);
1144: if ( !stat ) {
1145: NEXT(l) = d; d = l;
1146: d = nd_reconstruct(d);
1147: goto again;
1148: } else if ( nf ) {
1149: printf("+"); fflush(stdout);
1150: nh = nd_newps(nf);
1151: d = update_pairs(d,g,nh);
1152: g = update_base(g,nh);
1153: gall = append_one(gall,nh);
1154: FREENDP(l);
1155: } else {
1156: printf("."); fflush(stdout);
1157: FREENDP(l);
1158: }
1159: }
1160: return g;
1161: }
1162:
1163: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t)
1164: {
1165: ND_pairs d1,nd,cur,head,prev,remove;
1166:
1167: if ( !g ) return d;
1168: d = crit_B(d,t);
1169: d1 = nd_newpairs(g,t);
1170: d1 = crit_M(d1);
1171: d1 = crit_F(d1);
1172: prev = 0; cur = head = d1;
1173: while ( cur ) {
1174: if ( crit_2( cur->i1,cur->i2 ) ) {
1175: remove = cur;
1176: if ( !prev ) {
1177: head = cur = NEXT(cur);
1178: } else {
1179: cur = NEXT(prev) = NEXT(cur);
1180: }
1181: FREENDP(remove);
1182: } else {
1183: prev = cur;
1184: cur = NEXT(cur);
1185: }
1186: }
1187: if ( !d )
1188: return head;
1189: else {
1190: nd = d;
1191: while ( NEXT(nd) )
1192: nd = NEXT(nd);
1193: NEXT(nd) = head;
1194: return d;
1195: }
1196: }
1197:
1198: ND_pairs nd_newpairs( NODE g, int t )
1199: {
1200: NODE h;
1201: unsigned int *dl;
1202: int td,ts,s;
1203: ND_pairs r,r0;
1204:
1205: dl = HDL(nd_ps[t]);
1206: td = HTD(nd_ps[t]);
1207: ts = nd_ps[t]->sugar - td;
1208: for ( r0 = 0, h = g; h; h = NEXT(h) ) {
1209: NEXTND_pairs(r0,r);
1210: r->i1 = (int)BDY(h);
1211: r->i2 = t;
1212: ndl_lcm(HDL(nd_ps[r->i1]),dl,r->lcm);
1213: r->td = ndl_td(r->lcm);
1214: s = nd_ps[r->i1]->sugar-HTD(nd_ps[r->i1]);
1215: r->sugar = MAX(s,ts) + r->td;
1216: }
1217: NEXT(r) = 0;
1218: return r0;
1219: }
1220:
1221: ND_pairs crit_B( ND_pairs d, int s )
1222: {
1223: ND_pairs cur,head,prev,remove;
1224: unsigned int *t,*tl,*lcm;
1225: int td,tdl;
1226:
1227: if ( !d ) return 0;
1228: t = HDL(nd_ps[s]);
1229: prev = 0;
1230: head = cur = d;
1231: lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1232: while ( cur ) {
1233: tl = cur->lcm;
1234: if ( ndl_reducible(tl,t)
1235: && (ndl_lcm(HDL(nd_ps[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
1236: && (ndl_lcm(HDL(nd_ps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
1237: remove = cur;
1238: if ( !prev ) {
1239: head = cur = NEXT(cur);
1240: } else {
1241: cur = NEXT(prev) = NEXT(cur);
1242: }
1243: FREENDP(remove);
1244: } else {
1245: prev = cur;
1246: cur = NEXT(cur);
1247: }
1248: }
1249: return head;
1250: }
1251:
1252: ND_pairs crit_M( ND_pairs d1 )
1253: {
1254: ND_pairs e,d2,d3,dd,p;
1255: unsigned int *id,*jd;
1256: int itd,jtd;
1257:
1258: for ( dd = 0, e = d1; e; e = d3 ) {
1259: if ( !(d2 = NEXT(e)) ) {
1260: NEXT(e) = dd;
1261: return e;
1262: }
1263: id = e->lcm;
1264: itd = e->td;
1265: for ( d3 = 0; d2; d2 = p ) {
1266: p = NEXT(d2),
1267: jd = d2->lcm;
1268: jtd = d2->td;
1269: if ( jtd == itd )
1270: if ( id == jd );
1271: else if ( ndl_reducible(jd,id) ) continue;
1272: else if ( ndl_reducible(id,jd) ) goto delit;
1273: else ;
1274: else if ( jtd > itd )
1275: if ( ndl_reducible(jd,id) ) continue;
1276: else ;
1277: else if ( ndl_reducible(id,jd ) ) goto delit;
1278: NEXT(d2) = d3;
1279: d3 = d2;
1280: }
1281: NEXT(e) = dd;
1282: dd = e;
1283: continue;
1284: /**/
1285: delit: NEXT(d2) = d3;
1286: d3 = d2;
1287: for ( ; p; p = d2 ) {
1288: d2 = NEXT(p);
1289: NEXT(p) = d3;
1290: d3 = p;
1291: }
1292: FREENDP(e);
1293: }
1294: return dd;
1295: }
1296:
1297: ND_pairs crit_F( ND_pairs d1 )
1298: {
1299: ND_pairs rest, head,remove;
1300: ND_pairs last, p, r, w;
1301: int s;
1302:
1303: for ( head = last = 0, p = d1; NEXT(p); ) {
1304: r = w = equivalent_pairs(p,&rest);
1305: s = r->sugar;
1306: w = NEXT(w);
1307: while ( w ) {
1308: if ( crit_2(w->i1,w->i2) ) {
1309: r = w;
1310: w = NEXT(w);
1311: while ( w ) {
1312: remove = w;
1313: w = NEXT(w);
1314: FREENDP(remove);
1315: }
1316: break;
1317: } else if ( w->sugar < s ) {
1318: FREENDP(r);
1319: r = w;
1320: s = r->sugar;
1321: w = NEXT(w);
1322: } else {
1323: remove = w;
1324: w = NEXT(w);
1325: FREENDP(remove);
1326: }
1327: }
1328: if ( last ) NEXT(last) = r;
1329: else head = r;
1330: NEXT(last = r) = 0;
1331: p = rest;
1332: if ( !p ) return head;
1333: }
1334: if ( !last ) return p;
1335: NEXT(last) = p;
1336: return head;
1337: }
1338:
1339: int crit_2( int dp1, int dp2 )
1340: {
1341: return ndl_disjoint(HDL(nd_ps[dp1]),HDL(nd_ps[dp2]));
1342: }
1343:
1344: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest )
1345: {
1346: ND_pairs w,p,r,s;
1347: unsigned int *d;
1348: int td;
1349:
1350: w = d1;
1351: d = w->lcm;
1352: td = w->td;
1353: s = NEXT(w);
1354: NEXT(w) = 0;
1355: for ( r = 0; s; s = p ) {
1356: p = NEXT(s);
1357: if ( td == s->td && ndl_equal(d,s->lcm) ) {
1358: NEXT(s) = w;
1359: w = s;
1360: } else {
1361: NEXT(s) = r;
1362: r = s;
1363: }
1364: }
1365: *prest = r;
1366: return w;
1367: }
1368:
1369: NODE update_base(NODE nd,int ndp)
1370: {
1371: unsigned int *dl, *dln;
1372: NODE last, p, head;
1373: int td,tdn;
1374:
1375: dl = HDL(nd_ps[ndp]);
1376: td = HTD(nd_ps[ndp]);
1377: for ( head = last = 0, p = nd; p; ) {
1378: dln = HDL(nd_ps[(int)BDY(p)]);
1379: tdn = HTD(nd_ps[(int)BDY(p)]);
1380: if ( tdn >= td && ndl_reducible( dln, dl ) ) {
1381: p = NEXT(p);
1382: if ( last ) NEXT(last) = p;
1383: } else {
1384: if ( !last ) head = p;
1385: p = NEXT(last = p);
1386: }
1387: }
1388: head = append_one(head,ndp);
1389: return head;
1390: }
1391:
1392: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
1393: {
1394: ND_pairs m,ml,p,l;
1395: unsigned int *lcm;
1396: int s,td,len,tlen,c;
1397:
1398: if ( !(p = NEXT(m = d)) ) {
1399: *prest = p;
1400: NEXT(m) = 0;
1401: return m;
1402: }
1403: lcm = m->lcm;
1404: s = m->sugar;
1405: td = m->td;
1.6 noro 1406: len = nd_psl[m->i1]+nd_psl[m->i2];
1.1 noro 1407: for ( ml = 0, l = m; p; p = NEXT(l = p) ) {
1408: if (p->sugar < s)
1409: goto find;
1410: else if ( p->sugar == s ) {
1411: if ( p->td < td )
1412: goto find;
1413: else if ( p->td == td ) {
1414: c = ndl_compare(p->lcm,lcm);
1415: if ( c < 0 )
1416: goto find;
1.10 noro 1417: #if 0
1.1 noro 1418: else if ( c == 0 ) {
1.6 noro 1419: tlen = nd_psl[p->i1]+nd_psl[p->i2];
1.1 noro 1420: if ( tlen < len )
1421: goto find;
1422: }
1.10 noro 1423: #endif
1.1 noro 1424: }
1425: }
1426: continue;
1427: find:
1428: ml = l;
1429: m = p;
1430: lcm = m->lcm;
1431: s = m->sugar;
1432: td = m->td;
1433: len = tlen;
1434: }
1435: if ( !ml ) *prest = NEXT(m);
1436: else {
1437: NEXT(ml) = NEXT(m);
1438: *prest = d;
1439: }
1440: NEXT(m) = 0;
1441: return m;
1442: }
1443:
1444: int nd_newps(ND a)
1445: {
1.3 noro 1446: int len;
1.11 noro 1447: NM m;
1.3 noro 1448:
1.1 noro 1449: if ( nd_psn == nd_pslen ) {
1450: nd_pslen *= 2;
1.6 noro 1451: nd_psl = (int *)REALLOC((char *)nd_psl,nd_pslen*sizeof(int));
1.4 noro 1452: #if USE_NDV
1.11 noro 1453: nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV));
1454: #else
1455: nd_ps = (ND *)REALLOC((char *)nd_ps,nd_pslen*sizeof(ND));
1.4 noro 1456: #endif
1.11 noro 1457: nd_psh = (NM *)REALLOC((char *)nd_psh,nd_pslen*sizeof(NM));
1.1 noro 1458: nd_bound = (unsigned int **)
1459: REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
1460: }
1461: nd_monic(a);
1462: nd_bound[nd_psn] = nd_compute_bound(a);
1.11 noro 1463: NEWNM(m); m->td = HTD(a); ndl_copy(HDL(a),m->dl); nd_psh[nd_psn] = m;
1.5 noro 1464: #if USE_NDV
1.11 noro 1465: nd_ps[nd_psn]= ndtondv(a);
1466: nd_psl[nd_psn] = len = nd_ps[nd_psn]->len;
1.3 noro 1467: if ( len > nmv_len ) {
1468: nmv_len = 2*len;
1469: BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
1470: }
1.11 noro 1471: #else
1472: nd_ps[nd_psn] = a;
1473: nd_psl[nd_psn] = nd_length(a);
1.4 noro 1474: #endif
1.1 noro 1475: return nd_psn++;
1476: }
1477:
1478: NODE nd_setup(NODE f)
1479: {
1.5 noro 1480: int i,j,td,len,max;
1.1 noro 1481: NODE s,s0,f0;
1.5 noro 1482: unsigned int *d;
1.11 noro 1483: NM m;
1484:
1485: nd_found = 0; nd_notfirst = 0; nd_create = 0;
1.1 noro 1486:
1487: nd_psn = length(f); nd_pslen = 2*nd_psn;
1.6 noro 1488: nd_psl = (int *)MALLOC(nd_pslen*sizeof(int));
1.4 noro 1489: #if USE_NDV
1.11 noro 1490: nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
1491: #else
1492: nd_ps = (ND *)MALLOC(nd_pslen*sizeof(ND));
1.4 noro 1493: #endif
1.11 noro 1494: nd_psh = (NM *)MALLOC(nd_pslen*sizeof(NM));
1.1 noro 1495: nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
1.5 noro 1496: for ( max = 0, i = 0, s = f; i < nd_psn; i++, s = NEXT(s) ) {
1497: nd_bound[i] = d = dp_compute_bound((DP)BDY(s));
1498: for ( j = 0; j < nd_nvar; j++ )
1499: max = MAX(d[j],max);
1500: }
1.11 noro 1501: if ( !nd_red )
1502: nd_red = (NM *)MALLOC(REDTAB_LEN*sizeof(NM));
1503: bzero(nd_red,REDTAB_LEN*sizeof(NM));
1.5 noro 1504:
1505: if ( max < 2 )
1506: nd_bpe = 2;
1507: else if ( max < 4 )
1508: nd_bpe = 4;
1509: else if ( max < 64 )
1510: nd_bpe = 6;
1511: else if ( max < 256 )
1512: nd_bpe = 8;
1513: else if ( max < 65536 )
1514: nd_bpe = 16;
1515: else
1516: nd_bpe = 32;
1517:
1.1 noro 1518: nd_setup_parameters();
1519: nd_free_private_storage();
1.3 noro 1520: len = 0;
1.1 noro 1521: for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
1.11 noro 1522: #if USE_NDV
1523: nd_ps[i] = dptondv((DP)BDY(f));
1524: ndv_mul_c(nd_ps[i],invm(HC(nd_ps[i]),nd_mod));
1525: len = MAX(len,nd_ps[i]->len);
1526: #else
1.1 noro 1527: nd_ps[i] = dptond((DP)BDY(f));
1.11 noro 1528: nd_mul_c(nd_ps[i],1);
1.6 noro 1529: nd_psl[i] = nd_length(nd_ps[i]);
1.4 noro 1530: #endif
1.11 noro 1531: NEWNM(m); m->td = HTD(nd_ps[i]); ndl_copy(HDL(nd_ps[i]),m->dl);
1532: nd_psh[i] = m;
1.1 noro 1533: }
1.4 noro 1534: #if USE_NDV
1.3 noro 1535: nmv_len = 16*len;
1536: NEWNDV(ndv_red);
1537: BDY(ndv_red) = (NMV)MALLOC_ATOMIC(nmv_len*nmv_adv);
1.4 noro 1538: #endif
1.1 noro 1539: for ( s0 = 0, i = 0; i < nd_psn; i++ ) {
1540: NEXTNODE(s0,s); BDY(s) = (pointer)i;
1541: }
1542: if ( s0 ) NEXT(s) = 0;
1543: return s0;
1544: }
1545:
1546: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
1547: {
1548: struct order_spec ord1;
1549: VL fv,vv,vc;
1550: NODE fd,fd0,r,r0,t,x,s,xx;
1551: DP a,b,c;
1552:
1553: get_vars((Obj)f,&fv); pltovl(v,&vv);
1554: nd_nvar = length(vv);
1555: if ( ord->id )
1556: error("nd_gr : unsupported order");
1557: switch ( ord->ord.simple ) {
1558: case 0:
1559: is_rlex = 1;
1560: break;
1561: case 1:
1562: is_rlex = 0;
1563: break;
1564: default:
1565: error("nd_gr : unsupported order");
1566: }
1567: initd(ord);
1568: nd_mod = m;
1569: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
1570: ptod(CO,vv,(P)BDY(t),&b);
1571: _dp_mod(b,m,0,&c);
1572: if ( c ) {
1573: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
1574: }
1575: }
1576: if ( fd0 ) NEXT(fd) = 0;
1577: s = nd_setup(fd0);
1578: x = nd_gb(s);
1579: #if 0
1580: x = nd_reduceall(x,m);
1581: #endif
1582: for ( r0 = 0; x; x = NEXT(x) ) {
1583: NEXTNODE(r0,r);
1.11 noro 1584: #if USE_NDV
1585: a = ndvtodp(nd_ps[(int)BDY(x)]);
1586: #else
1.1 noro 1587: a = ndtodp(nd_ps[(int)BDY(x)]);
1.11 noro 1588: #endif
1.1 noro 1589: _dtop_mod(CO,vv,a,(P *)&BDY(r));
1590: }
1591: if ( r0 ) NEXT(r) = 0;
1592: MKLIST(*rp,r0);
1593: fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
1594: nd_found,nd_notfirst,nd_create);
1595: }
1596:
1597: void dltondl(int n,DL dl,unsigned int *r)
1598: {
1599: unsigned int *d;
1600: int i;
1601:
1602: d = dl->d;
1603: bzero(r,nd_wpd*sizeof(unsigned int));
1604: if ( is_rlex )
1605: for ( i = 0; i < n; i++ )
1606: r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe));
1607: else
1608: for ( i = 0; i < n; i++ )
1609: r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe);
1610: }
1611:
1612: DL ndltodl(int n,int td,unsigned int *ndl)
1613: {
1614: DL dl;
1615: int *d;
1616: int i;
1617:
1618: NEWDL(dl,n);
1619: dl->td = td;
1620: d = dl->d;
1621: if ( is_rlex )
1622: for ( i = 0; i < n; i++ )
1623: d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
1624: &((1<<nd_bpe)-1);
1625: else
1626: for ( i = 0; i < n; i++ )
1627: d[i] = (ndl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
1628: &((1<<nd_bpe)-1);
1629: return dl;
1630: }
1631:
1632: ND dptond(DP p)
1633: {
1634: ND d;
1635: NM m0,m;
1636: MP t;
1637: int n;
1638:
1639: if ( !p )
1640: return 0;
1641: n = NV(p);
1642: m0 = 0;
1643: for ( t = BDY(p); t; t = NEXT(t) ) {
1644: NEXTNM(m0,m);
1645: m->c = ITOS(t->c);
1646: m->td = t->dl->td;
1647: dltondl(n,t->dl,m->dl);
1648: }
1649: NEXT(m) = 0;
1650: MKND(n,m0,d);
1651: d->nv = n;
1652: d->sugar = p->sugar;
1653: return d;
1654: }
1655:
1656: DP ndtodp(ND p)
1657: {
1658: DP d;
1659: MP m0,m;
1660: NM t;
1661: int n;
1662:
1663: if ( !p )
1664: return 0;
1665: n = NV(p);
1666: m0 = 0;
1667: for ( t = BDY(p); t; t = NEXT(t) ) {
1668: NEXTMP(m0,m);
1669: m->c = STOI(t->c);
1670: m->dl = ndltodl(n,t->td,t->dl);
1671: }
1672: NEXT(m) = 0;
1673: MKDP(n,m0,d);
1674: d->sugar = p->sugar;
1675: return d;
1676: }
1677:
1678: void ndl_print(unsigned int *dl)
1679: {
1680: int n;
1681: int i;
1682:
1683: n = nd_nvar;
1684: printf("<<");
1685: if ( is_rlex )
1686: for ( i = 0; i < n; i++ )
1687: printf(i==n-1?"%d":"%d,",
1688: (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
1689: &((1<<nd_bpe)-1));
1690: else
1691: for ( i = 0; i < n; i++ )
1692: printf(i==n-1?"%d":"%d,",
1693: (dl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
1694: &((1<<nd_bpe)-1));
1695: printf(">>");
1696: }
1697:
1698: void nd_print(ND p)
1699: {
1700: NM m;
1701:
1702: if ( !p )
1703: printf("0\n");
1704: else {
1705: for ( m = BDY(p); m; m = NEXT(m) ) {
1706: printf("+%d*",m->c);
1707: ndl_print(m->dl);
1708: }
1709: printf("\n");
1710: }
1711: }
1712:
1713: void ndp_print(ND_pairs d)
1714: {
1715: ND_pairs t;
1716:
1717: for ( t = d; t; t = NEXT(t) ) {
1718: printf("%d,%d ",t->i1,t->i2);
1719: }
1720: printf("\n");
1721: }
1722:
1723: void nd_monic(ND p)
1724: {
1725: if ( !p )
1726: return;
1727: else
1728: nd_mul_c(p,invm(HC(p),nd_mod));
1729: }
1730:
1731: void nd_mul_c(ND p,int mul)
1732: {
1733: NM m;
1734: int c,c1;
1735:
1736: if ( !p )
1737: return;
1738: for ( m = BDY(p); m; m = NEXT(m) ) {
1739: c1 = C(m);
1740: DMAR(c1,mul,0,nd_mod,c);
1741: C(m) = c;
1742: }
1743: }
1744:
1745: void nd_free(ND p)
1746: {
1747: NM t,s;
1748:
1749: if ( !p )
1750: return;
1751: t = BDY(p);
1752: while ( t ) {
1753: s = NEXT(t);
1754: FREENM(t);
1755: t = s;
1756: }
1757: FREEND(p);
1758: }
1759:
1760: void nd_append_red(unsigned int *d,int td,int i)
1761: {
1762: NM m,m0;
1763: int h;
1764:
1765: NEWNM(m);
1766: h = ndl_hash_value(td,d);
1767: m->c = i;
1768: m->td = td;
1.7 noro 1769: ndl_copy(d,m->dl);
1.1 noro 1770: NEXT(m) = nd_red[h];
1771: nd_red[h] = m;
1772: }
1773:
1.5 noro 1774: unsigned int *dp_compute_bound(DP p)
1775: {
1776: unsigned int *d,*d1,*d2,*t;
1777: MP m;
1.7 noro 1778: int i,l;
1.5 noro 1779:
1780: if ( !p )
1781: return 0;
1782: d1 = (unsigned int *)ALLOCA(nd_nvar*sizeof(unsigned int));
1783: d2 = (unsigned int *)ALLOCA(nd_nvar*sizeof(unsigned int));
1784: m = BDY(p);
1785: bcopy(m->dl->d,d1,nd_nvar*sizeof(unsigned int));
1786: for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
1787: d = m->dl->d;
1788: for ( i = 0; i < nd_nvar; i++ )
1789: d2[i] = d[i] > d1[i] ? d[i] : d1[i];
1790: t = d1; d1 = d2; d2 = t;
1791: }
1.12 ! noro 1792: l = nd_nvar+31;
1.7 noro 1793: t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
1794: bzero(t,l*sizeof(unsigned int));
1.5 noro 1795: bcopy(d1,t,nd_nvar*sizeof(unsigned int));
1796: return t;
1797: }
1798:
1.1 noro 1799: unsigned int *nd_compute_bound(ND p)
1800: {
1801: unsigned int *d1,*d2,*t;
1.9 noro 1802: int i,l;
1.1 noro 1803: NM m;
1804:
1805: if ( !p )
1806: return 0;
1807: d1 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1808: d2 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1809: bcopy(HDL(p),d1,nd_wpd*sizeof(unsigned int));
1810: for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
1811: ndl_lcm(m->dl,d1,d2);
1812: t = d1; d1 = d2; d2 = t;
1813: }
1.12 ! noro 1814: l = nd_nvar+31;
1.9 noro 1815: t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
1816: bzero(t,l*sizeof(unsigned int));
1.5 noro 1817: for ( i = 0; i < nd_nvar; i++ )
1818: t[i] = (d1[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))&nd_mask0;
1.1 noro 1819: return t;
1820: }
1821:
1822: void nd_setup_parameters() {
1823: int i;
1824:
1825: nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
1826: nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
1827: if ( nd_bpe < 32 ) {
1828: nd_mask0 = (1<<nd_bpe)-1;
1829: } else {
1830: nd_mask0 = 0xffffffff;
1831: }
1832: bzero(nd_mask,sizeof(nd_mask));
1833: nd_mask1 = 0;
1834: for ( i = 0; i < nd_epw; i++ ) {
1835: nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
1836: nd_mask1 |= (1<<(nd_bpe-1))<<(i*nd_bpe);
1837: }
1.4 noro 1838: #if USE_NDV
1.3 noro 1839: nmv_adv = sizeof(struct oNMV)+(nd_wpd-1)*sizeof(unsigned int);
1.4 noro 1840: #endif
1.1 noro 1841: }
1842:
1843: ND_pairs nd_reconstruct(ND_pairs d)
1844: {
1.11 noro 1845: int i,obpe,oadv;
1.6 noro 1846: NM prev_nm_free_list,mr0,mr,m;
1.1 noro 1847: ND_pairs s0,s,t,prev_ndp_free_list;
1848:
1849: obpe = nd_bpe;
1.11 noro 1850: #if USE_NDV
1851: oadv = nmv_adv;
1852: #endif
1.5 noro 1853: if ( obpe < 4 )
1854: nd_bpe = 4;
1855: else if ( obpe < 6 )
1856: nd_bpe = 6;
1857: else if ( obpe < 8 )
1858: nd_bpe = 8;
1859: else if ( obpe < 16 )
1860: nd_bpe = 16;
1861: else if ( obpe < 32 )
1862: nd_bpe = 32;
1863: else
1864: error("nd_reconstruct : exponent too large");
1865:
1.1 noro 1866: nd_setup_parameters();
1867: prev_nm_free_list = _nm_free_list;
1868: prev_ndp_free_list = _ndp_free_list;
1869: _nm_free_list = 0;
1870: _ndp_free_list = 0;
1871: for ( i = 0; i < nd_psn; i++ ) {
1.4 noro 1872: #if USE_NDV
1.11 noro 1873: ndv_realloc(nd_ps[i],obpe,oadv);
1874: #else
1875: nd_realloc(nd_ps[i],obpe);
1.4 noro 1876: #endif
1.1 noro 1877: }
1878: s0 = 0;
1879: for ( t = d; t; t = NEXT(t) ) {
1880: NEXTND_pairs(s0,s);
1881: s->i1 = t->i1;
1882: s->i2 = t->i2;
1883: s->td = t->td;
1884: s->sugar = t->sugar;
1885: ndl_dup(obpe,t->lcm,s->lcm);
1886: }
1.6 noro 1887: for ( i = 0; i < REDTAB_LEN; i++ ) {
1888: for ( mr0 = 0, m = nd_red[i]; m; m = NEXT(m) ) {
1889: NEXTNM(mr0,mr);
1890: mr->c = m->c;
1891: mr->td = m->td;
1892: ndl_dup(obpe,m->dl,mr->dl);
1893: }
1894: if ( mr0 ) NEXT(mr) = 0;
1895: nd_red[i] = mr0;
1896: }
1.11 noro 1897: for ( i = 0; i < nd_psn; i++ ) {
1898: NEWNM(m); m->td = nd_psh[i]->td; ndl_dup(obpe,nd_psh[i]->dl,m->dl);
1899: nd_psh[i] = m;
1900: }
1.1 noro 1901: if ( s0 ) NEXT(s) = 0;
1902: prev_nm_free_list = 0;
1903: prev_ndp_free_list = 0;
1.4 noro 1904: #if USE_NDV
1.3 noro 1905: BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
1.4 noro 1906: #endif
1.1 noro 1907: GC_gcollect();
1908: return s0;
1909: }
1910:
1911: void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
1912: {
1913: int n,i,ei,oepw,cepw,cbpe;
1914:
1915: n = nd_nvar;
1916: oepw = (sizeof(unsigned int)*8)/obpe;
1917: cepw = nd_epw;
1918: cbpe = nd_bpe;
1919: if ( is_rlex )
1920: for ( i = 0; i < n; i++ ) {
1921: ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe))
1922: &((1<<obpe)-1);
1923: r[(n-1-i)/cepw] |= (ei<<((cepw-((n-1-i)%cepw)-1)*cbpe));
1924: }
1925: else
1926: for ( i = 0; i < n; i++ ) {
1927: ei = (d[i/oepw]>>((oepw-(i%oepw)-1)*obpe))
1928: &((1<<obpe)-1);
1929: r[i/cepw] |= (ei<<((cepw-(i%cepw)-1)*cbpe));
1930: }
1931: }
1932:
1.11 noro 1933: void nd_realloc(ND p,int obpe)
1.1 noro 1934: {
1935: NM m,mr,mr0;
1936:
1.11 noro 1937: if ( p ) {
1938: m = BDY(p);
1.1 noro 1939: for ( mr0 = 0; m; m = NEXT(m) ) {
1940: NEXTNM(mr0,mr);
1941: C(mr) = C(m);
1942: mr->td = m->td;
1943: ndl_dup(obpe,m->dl,mr->dl);
1944: }
1945: NEXT(mr) = 0;
1.11 noro 1946: BDY(p) = mr0;
1.1 noro 1947: }
1948: }
1.3 noro 1949:
1.6 noro 1950: ND nd_copy(ND p)
1951: {
1952: NM m,mr,mr0;
1953: int c,n,s;
1954: ND r;
1955:
1956: if ( !p )
1957: return 0;
1958: else {
1959: s = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);
1960: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1961: NEXTNM(mr0,mr);
1962: C(mr) = C(m);
1963: mr->td = m->td;
1964: ndl_copy(m->dl,mr->dl);
1965: }
1966: NEXT(mr) = 0;
1967: MKND(NV(p),mr0,r);
1968: r->sugar = p->sugar;
1969: return r;
1970: }
1971: }
1972:
1.4 noro 1973: #if USE_NDV
1.11 noro 1974: int nd_sp(ND_pairs p,ND *rp)
1975: {
1976: NM m;
1977: NDV p1,p2;
1978: ND t1,t2;
1979: unsigned int *lcm;
1980: int td;
1981:
1982: p1 = nd_ps[p->i1];
1983: p2 = nd_ps[p->i2];
1984: lcm = p->lcm;
1985: td = p->td;
1986: NEWNM(m);
1987: C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl);
1988: if ( ndl_check_bound2(p->i1,m->dl) )
1989: return 0;
1990: t1 = ndv_mul_nm_create(p1,m);
1991: C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
1992: if ( ndl_check_bound2(p->i2,m->dl) ) {
1993: nd_free(t1);
1994: return 0;
1995: }
1996: ndv_mul_nm(p2,m,ndv_red);
1997: FREENM(m);
1998: *rp = ndv_add(t1,ndv_red);
1999: return 1;
2000: }
2001:
2002: void ndv_mul_c(NDV p,int mul)
2003: {
2004: NMV m;
2005: int c,c1,len,i;
2006:
2007: if ( !p )
2008: return;
2009: len = p->len;
2010: for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
2011: c1 = C(m);
2012: DMAR(c1,mul,0,nd_mod,c);
2013: C(m) = c;
2014: }
2015: }
2016:
1.4 noro 2017: void ndv_mul_nm(NDV p,NM m0,NDV r)
2018: {
2019: NMV m,mr,mr0;
2020: unsigned int *d,*dt,*dm;
2021: int c,n,td,i,c1,c2,len;
2022:
2023: if ( !p )
2024: /* XXX */
2025: r->len = 0;
2026: else {
2027: n = NV(p); m = BDY(p); len = p->len;
2028: d = m0->dl; td = m0->td; c = C(m0);
2029: mr = BDY(r);
1.9 noro 2030: for ( ; len > 0; len--, NMV_ADV(m), NMV_ADV(mr) ) {
2031: c1 = C(m); DMAR(c1,c,0,nd_mod,c2); C(mr) = c2;
2032: mr->td = m->td+td; ndl_add(m->dl,d,mr->dl);
2033: }
2034: NV(r) = NV(p);
2035: r->len = p->len;
2036: r->sugar = p->sugar + td;
2037: }
2038: }
2039:
2040: ND ndv_mul_nm_create(NDV p,NM m0)
2041: {
2042: NM mr,mr0;
2043: NMV m;
2044: unsigned int *d,*dt,*dm;
2045: int c,n,td,i,c1,c2,len;
2046: ND r;
2047:
2048: if ( !p )
2049: return 0;
2050: else {
2051: n = NV(p); m = BDY(p);
2052: d = m0->dl; td = m0->td; c = C(m0);
2053: len = p->len;
2054: mr0 = 0;
2055: for ( i = 0; i < len; i++, NMV_ADV(m) ) {
2056: NEXTNM(mr0,mr);
1.4 noro 2057: c1 = C(m);
2058: DMAR(c1,c,0,nd_mod,c2);
2059: C(mr) = c2;
2060: mr->td = m->td+td;
2061: ndl_add(m->dl,d,mr->dl);
2062: }
1.9 noro 2063: NEXT(mr) = 0;
2064: MKND(NV(p),mr0,r);
1.4 noro 2065: r->sugar = p->sugar + td;
1.9 noro 2066: return r;
1.4 noro 2067: }
2068: }
2069:
2070: ND ndv_add(ND p1,NDV p2)
2071: {
1.9 noro 2072: register NM prev,cur,new;
1.4 noro 2073: int c,c1,c2,t,td,td2,mul,len,i;
1.9 noro 2074: NM head;
1.4 noro 2075: unsigned int *d;
2076: NMV m2;
2077:
2078: if ( !p1 )
2079: return 0;
2080: else {
2081: prev = 0; head = cur = BDY(p1);
1.9 noro 2082: NEWNM(new); len = p2->len;
2083: for ( m2 = BDY(p2), i = 0; cur && i < len; ) {
1.4 noro 2084: td2 = new->td = m2->td;
1.9 noro 2085: if ( cur->td > td2 ) c = 1;
2086: else if ( cur->td < td2 ) c = -1;
2087: else c = ndl_compare(cur->dl,m2->dl);
1.4 noro 2088: switch ( c ) {
2089: case 0:
2090: t = C(m2)+C(cur)-nd_mod;
1.9 noro 2091: if ( t < 0 ) t += nd_mod;
2092: if ( t ) C(cur) = t;
1.4 noro 2093: else if ( !prev ) {
1.9 noro 2094: head = NEXT(cur); FREENM(cur); cur = head;
1.4 noro 2095: } else {
1.9 noro 2096: NEXT(prev) = NEXT(cur); FREENM(cur); cur = NEXT(prev);
1.4 noro 2097: }
2098: NMV_ADV(m2); i++;
2099: break;
2100: case 1:
1.9 noro 2101: prev = cur; cur = NEXT(cur);
1.4 noro 2102: break;
2103: case -1:
1.6 noro 2104: ndl_copy(m2->dl,new->dl);
1.9 noro 2105: C(new) = C(m2);
1.4 noro 2106: if ( !prev ) {
2107: /* cur = head */
1.9 noro 2108: prev = new; NEXT(prev) = head; head = prev;
1.4 noro 2109: } else {
1.9 noro 2110: NEXT(prev) = new; NEXT(new) = cur; prev = new;
1.4 noro 2111: }
1.9 noro 2112: NEWNM(new); NMV_ADV(m2); i++;
1.4 noro 2113: break;
2114: }
2115: }
1.9 noro 2116: for ( ; i < len; i++, NMV_ADV(m2) ) {
2117: td2 = new->td = m2->td; C(new) = C(m2); ndl_copy(m2->dl,new->dl);
2118: if ( !prev ) {
2119: prev = new; NEXT(prev) = 0; head = prev;
2120: } else {
2121: NEXT(prev) = new; NEXT(new) = 0; prev = new;
2122: }
2123: NEWNM(new);
2124: }
1.4 noro 2125: FREENM(new);
2126: if ( head ) {
1.9 noro 2127: BDY(p1) = head; p1->sugar = MAX(p1->sugar,p2->sugar);
1.4 noro 2128: return p1;
2129: } else {
2130: FREEND(p1);
2131: return 0;
2132: }
2133:
2134: }
2135: }
2136:
1.11 noro 2137: void ndv_realloc(NDV p,int obpe,int oadv)
2138: {
2139: NMV m,mr,mr0;
2140: int len,i;
2141:
2142: #define NMV_OADV(m) (m = (NMV)(((char *)m)+oadv))
2143:
2144: if ( p ) {
2145: m = BDY(p); len = p->len;
2146: mr0 = mr = (NMV)MALLOC_ATOMIC(len*nmv_adv);
2147: bzero(mr0,len*nmv_adv);
2148: for ( i = 0; i < len; i++, NMV_OADV(m), NMV_ADV(mr) ) {
2149: C(mr) = C(m);
2150: mr->td = m->td;
2151: ndl_dup(obpe,m->dl,mr->dl);
2152: }
2153: BDY(p) = mr0;
2154: }
2155: }
2156:
1.3 noro 2157: NDV ndtondv(ND p)
2158: {
2159: NDV d;
2160: NMV m,m0;
2161: NM t;
2162: int i,len;
2163:
2164: if ( !p )
2165: return 0;
2166: len = nd_length(p);
2167: m0 = m = (NMV)MALLOC_ATOMIC(len*nmv_adv);
2168: for ( t = BDY(p), i = 0; t; t = NEXT(t), i++, NMV_ADV(m) ) {
2169: m->td = t->td;
1.7 noro 2170: ndl_copy(t->dl,m->dl);
1.3 noro 2171: m->c = t->c;
2172: }
2173: MKNDV(NV(p),m0,len,d);
2174: d->sugar = p->sugar;
2175: return d;
2176: }
2177:
1.11 noro 2178: NDV dptondv(DP p)
2179: {
2180: NDV d;
2181: NMV m0,m;
2182: MP t;
2183: int l,i,n;
2184:
2185: if ( !p )
2186: return 0;
2187: for ( t = BDY(p), l = 0; t; t = NEXT(t), l++ );
2188: m0 = m = (NMV)MALLOC_ATOMIC(l*nmv_adv);
2189: n = NV(p);
2190: for ( t = BDY(p), i = 0; i < l; i++, t = NEXT(t), NMV_ADV(m) ) {
2191: m->c = ITOS(t->c);
2192: m->td = t->dl->td;
2193: dltondl(n,t->dl,m->dl);
2194: }
2195: MKNDV(n,m0,l,d);
2196: d->sugar = p->sugar;
2197: return d;
2198: }
2199:
2200: DP ndvtodp(NDV p)
2201: {
2202: DP d;
2203: MP m0,m;
2204: NMV t;
2205: int len,i,n;
2206:
2207: if ( !p )
2208: return 0;
2209: m0 = 0;
2210: len = p->len;
2211: n = NV(p);
2212: for ( t = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) {
2213: NEXTMP(m0,m);
2214: m->c = STOI(t->c);
2215: m->dl = ndltodl(n,t->td,t->dl);
2216: }
2217: NEXT(m) = 0;
2218: MKDP(NV(p),m0,d);
2219: d->sugar = p->sugar;
2220: return d;
2221: }
2222:
1.3 noro 2223: void ndv_print(NDV p)
2224: {
2225: NMV m;
2226: int i,len;
2227:
2228: if ( !p )
2229: printf("0\n");
2230: else {
2231: len = p->len;
2232: for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
2233: printf("+%d*",m->c);
2234: ndl_print(m->dl);
2235: }
2236: printf("\n");
2237: }
1.11 noro 2238: }
2239: #else
2240: int nd_sp(ND_pairs p,ND *rp)
2241: {
2242: NM m;
2243: ND p1,p2;
2244: ND t1,t2;
2245: unsigned int *lcm;
2246: int td;
2247:
2248: p1 = nd_ps[p->i1];
2249: p2 = nd_ps[p->i2];
2250: lcm = p->lcm;
2251: td = p->td;
2252: NEWNM(m);
2253: C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl);
2254: if ( ndl_check_bound2(p->i1,m->dl) )
2255: return 0;
2256: t1 = nd_mul_ind_nm(p->i1,m);
2257: C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
2258: if ( ndl_check_bound2(p->i2,m->dl) ) {
2259: nd_free(t1);
2260: return 0;
2261: }
2262: t2 = nd_mul_ind_nm(p->i2,m);
2263: FREENM(m);
2264: *rp = nd_add(t1,t2);
2265: return 1;
2266: }
2267:
2268: ND nd_mul_nm(ND p,NM m0)
2269: {
2270: NM m,mr,mr0;
2271: unsigned int *d,*dt,*dm;
2272: int c,n,td,i,c1,c2;
2273: int *pt,*p1,*p2;
2274: ND r;
2275:
2276: if ( !p )
2277: return 0;
2278: else {
2279: n = NV(p); m = BDY(p);
2280: d = m0->dl; td = m0->td; c = C(m0);
2281: mr0 = 0;
2282: for ( ; m; m = NEXT(m) ) {
2283: NEXTNM(mr0,mr);
2284: c1 = C(m);
2285: DMAR(c1,c,0,nd_mod,c2);
2286: C(mr) = c2;
2287: mr->td = m->td+td;
2288: ndl_add(m->dl,d,mr->dl);
2289: }
2290: NEXT(mr) = 0;
2291: MKND(NV(p),mr0,r);
2292: r->sugar = p->sugar + td;
2293: return r;
2294: }
2295: }
2296:
2297: ND nd_mul_ind_nm(int index,NM m0)
2298: {
2299: register int c1,c2,c;
2300: register NM m,new,prev;
2301: NM mr0;
2302: unsigned int *d;
2303: int n,td,i,len,d0,d1;
2304: ND p,r;
2305:
2306: p = nd_ps[index];
2307: len = nd_psl[index];
2308: n = NV(p); m = BDY(p);
2309: d = m0->dl; td = m0->td; c = C(m0);
2310:
2311: NEWNM(mr0);
2312: c1 = C(m); DMAR(c1,c,0,nd_mod,c2); C(mr0) = c2;
2313: mr0->td = m->td+td; ndl_add(m->dl,d,mr0->dl);
2314: prev = mr0; m = NEXT(m);
2315: len--;
2316:
2317: switch ( nd_wpd ) {
2318: case 1:
2319: d0 = d[0];
2320: while ( len-- ) {
2321: c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
2322: NEWNM(new); C(new) = c2;
2323: new->td = m->td+td; new->dl[0] = m->dl[0]+d0;
2324: m = NEXT(m); NEXT(prev) = new; prev = new;
2325: }
2326: break;
2327: case 2:
2328: d0 = d[0]; d1 = d[1];
2329: while ( len-- ) {
2330: c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
2331: NEWNM(new); C(new) = c2;
2332: new->td = m->td+td;
2333: new->dl[0] = m->dl[0]+d0;
2334: new->dl[1] = m->dl[1]+d1;
2335: m = NEXT(m); NEXT(prev) = new; prev = new;
2336: }
2337: break;
2338: default:
2339: while ( len-- ) {
2340: c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
2341: NEWNM(new); C(new) = c2;
2342: new->td = m->td+td; ndl_add(m->dl,d,new->dl);
2343: m = NEXT(m); NEXT(prev) = new; prev = new;
2344: }
2345: break;
2346: }
2347:
2348: NEXT(prev) = 0;
2349: MKND(NV(p),mr0,r);
2350: r->sugar = p->sugar + td;
2351: return r;
1.3 noro 2352: }
1.4 noro 2353: #endif