Annotation of OpenXM_contrib2/asir2000/engine/nd.c, Revision 1.11
1.11 ! noro 1: /* $OpenXM$ */
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.9 noro 1792: l = ((nd_nvar+(sizeof(unsigned int)-1))
1793: /sizeof(unsigned int))*sizeof(unsigned int);
1.7 noro 1794: t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
1795: bzero(t,l*sizeof(unsigned int));
1.5 noro 1796: bcopy(d1,t,nd_nvar*sizeof(unsigned int));
1797: return t;
1798: }
1799:
1.1 noro 1800: unsigned int *nd_compute_bound(ND p)
1801: {
1802: unsigned int *d1,*d2,*t;
1.9 noro 1803: int i,l;
1.1 noro 1804: NM m;
1805:
1806: if ( !p )
1807: return 0;
1808: d1 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1809: d2 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1810: bcopy(HDL(p),d1,nd_wpd*sizeof(unsigned int));
1811: for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
1812: ndl_lcm(m->dl,d1,d2);
1813: t = d1; d1 = d2; d2 = t;
1814: }
1.9 noro 1815: l = ((nd_nvar+(sizeof(unsigned int)-1))
1816: /sizeof(unsigned int))*sizeof(unsigned int);
1817: t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
1818: bzero(t,l*sizeof(unsigned int));
1.5 noro 1819: for ( i = 0; i < nd_nvar; i++ )
1820: t[i] = (d1[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))&nd_mask0;
1.1 noro 1821: return t;
1822: }
1823:
1824: void nd_setup_parameters() {
1825: int i;
1826:
1827: nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
1828: nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
1829: if ( nd_bpe < 32 ) {
1830: nd_mask0 = (1<<nd_bpe)-1;
1831: } else {
1832: nd_mask0 = 0xffffffff;
1833: }
1834: bzero(nd_mask,sizeof(nd_mask));
1835: nd_mask1 = 0;
1836: for ( i = 0; i < nd_epw; i++ ) {
1837: nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
1838: nd_mask1 |= (1<<(nd_bpe-1))<<(i*nd_bpe);
1839: }
1.4 noro 1840: #if USE_NDV
1.3 noro 1841: nmv_adv = sizeof(struct oNMV)+(nd_wpd-1)*sizeof(unsigned int);
1.4 noro 1842: #endif
1.1 noro 1843: }
1844:
1845: ND_pairs nd_reconstruct(ND_pairs d)
1846: {
1.11 ! noro 1847: int i,obpe,oadv;
1.6 noro 1848: NM prev_nm_free_list,mr0,mr,m;
1.1 noro 1849: ND_pairs s0,s,t,prev_ndp_free_list;
1850:
1851: obpe = nd_bpe;
1.11 ! noro 1852: #if USE_NDV
! 1853: oadv = nmv_adv;
! 1854: #endif
1.5 noro 1855: if ( obpe < 4 )
1856: nd_bpe = 4;
1857: else if ( obpe < 6 )
1858: nd_bpe = 6;
1859: else if ( obpe < 8 )
1860: nd_bpe = 8;
1861: else if ( obpe < 16 )
1862: nd_bpe = 16;
1863: else if ( obpe < 32 )
1864: nd_bpe = 32;
1865: else
1866: error("nd_reconstruct : exponent too large");
1867:
1.1 noro 1868: nd_setup_parameters();
1869: prev_nm_free_list = _nm_free_list;
1870: prev_ndp_free_list = _ndp_free_list;
1871: _nm_free_list = 0;
1872: _ndp_free_list = 0;
1873: for ( i = 0; i < nd_psn; i++ ) {
1.4 noro 1874: #if USE_NDV
1.11 ! noro 1875: ndv_realloc(nd_ps[i],obpe,oadv);
! 1876: #else
! 1877: nd_realloc(nd_ps[i],obpe);
1.4 noro 1878: #endif
1.1 noro 1879: }
1880: s0 = 0;
1881: for ( t = d; t; t = NEXT(t) ) {
1882: NEXTND_pairs(s0,s);
1883: s->i1 = t->i1;
1884: s->i2 = t->i2;
1885: s->td = t->td;
1886: s->sugar = t->sugar;
1887: ndl_dup(obpe,t->lcm,s->lcm);
1888: }
1.6 noro 1889: for ( i = 0; i < REDTAB_LEN; i++ ) {
1890: for ( mr0 = 0, m = nd_red[i]; m; m = NEXT(m) ) {
1891: NEXTNM(mr0,mr);
1892: mr->c = m->c;
1893: mr->td = m->td;
1894: ndl_dup(obpe,m->dl,mr->dl);
1895: }
1896: if ( mr0 ) NEXT(mr) = 0;
1897: nd_red[i] = mr0;
1898: }
1.11 ! noro 1899: for ( i = 0; i < nd_psn; i++ ) {
! 1900: NEWNM(m); m->td = nd_psh[i]->td; ndl_dup(obpe,nd_psh[i]->dl,m->dl);
! 1901: nd_psh[i] = m;
! 1902: }
1.1 noro 1903: if ( s0 ) NEXT(s) = 0;
1904: prev_nm_free_list = 0;
1905: prev_ndp_free_list = 0;
1.4 noro 1906: #if USE_NDV
1.3 noro 1907: BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
1.4 noro 1908: #endif
1.1 noro 1909: GC_gcollect();
1910: return s0;
1911: }
1912:
1913: void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
1914: {
1915: int n,i,ei,oepw,cepw,cbpe;
1916:
1917: n = nd_nvar;
1918: oepw = (sizeof(unsigned int)*8)/obpe;
1919: cepw = nd_epw;
1920: cbpe = nd_bpe;
1921: if ( is_rlex )
1922: for ( i = 0; i < n; i++ ) {
1923: ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe))
1924: &((1<<obpe)-1);
1925: r[(n-1-i)/cepw] |= (ei<<((cepw-((n-1-i)%cepw)-1)*cbpe));
1926: }
1927: else
1928: for ( i = 0; i < n; i++ ) {
1929: ei = (d[i/oepw]>>((oepw-(i%oepw)-1)*obpe))
1930: &((1<<obpe)-1);
1931: r[i/cepw] |= (ei<<((cepw-(i%cepw)-1)*cbpe));
1932: }
1933: }
1934:
1.11 ! noro 1935: void nd_realloc(ND p,int obpe)
1.1 noro 1936: {
1937: NM m,mr,mr0;
1938:
1.11 ! noro 1939: if ( p ) {
! 1940: m = BDY(p);
1.1 noro 1941: for ( mr0 = 0; m; m = NEXT(m) ) {
1942: NEXTNM(mr0,mr);
1943: C(mr) = C(m);
1944: mr->td = m->td;
1945: ndl_dup(obpe,m->dl,mr->dl);
1946: }
1947: NEXT(mr) = 0;
1.11 ! noro 1948: BDY(p) = mr0;
1.1 noro 1949: }
1950: }
1.3 noro 1951:
1.6 noro 1952: ND nd_copy(ND p)
1953: {
1954: NM m,mr,mr0;
1955: int c,n,s;
1956: ND r;
1957:
1958: if ( !p )
1959: return 0;
1960: else {
1961: s = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);
1962: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1963: NEXTNM(mr0,mr);
1964: C(mr) = C(m);
1965: mr->td = m->td;
1966: ndl_copy(m->dl,mr->dl);
1967: }
1968: NEXT(mr) = 0;
1969: MKND(NV(p),mr0,r);
1970: r->sugar = p->sugar;
1971: return r;
1972: }
1973: }
1974:
1.4 noro 1975: #if USE_NDV
1.11 ! noro 1976: int nd_sp(ND_pairs p,ND *rp)
! 1977: {
! 1978: NM m;
! 1979: NDV p1,p2;
! 1980: ND t1,t2;
! 1981: unsigned int *lcm;
! 1982: int td;
! 1983:
! 1984: p1 = nd_ps[p->i1];
! 1985: p2 = nd_ps[p->i2];
! 1986: lcm = p->lcm;
! 1987: td = p->td;
! 1988: NEWNM(m);
! 1989: C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl);
! 1990: if ( ndl_check_bound2(p->i1,m->dl) )
! 1991: return 0;
! 1992: t1 = ndv_mul_nm_create(p1,m);
! 1993: C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
! 1994: if ( ndl_check_bound2(p->i2,m->dl) ) {
! 1995: nd_free(t1);
! 1996: return 0;
! 1997: }
! 1998: ndv_mul_nm(p2,m,ndv_red);
! 1999: FREENM(m);
! 2000: *rp = ndv_add(t1,ndv_red);
! 2001: return 1;
! 2002: }
! 2003:
! 2004: void ndv_mul_c(NDV p,int mul)
! 2005: {
! 2006: NMV m;
! 2007: int c,c1,len,i;
! 2008:
! 2009: if ( !p )
! 2010: return;
! 2011: len = p->len;
! 2012: for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
! 2013: c1 = C(m);
! 2014: DMAR(c1,mul,0,nd_mod,c);
! 2015: C(m) = c;
! 2016: }
! 2017: }
! 2018:
1.4 noro 2019: void ndv_mul_nm(NDV p,NM m0,NDV r)
2020: {
2021: NMV m,mr,mr0;
2022: unsigned int *d,*dt,*dm;
2023: int c,n,td,i,c1,c2,len;
2024:
2025: if ( !p )
2026: /* XXX */
2027: r->len = 0;
2028: else {
2029: n = NV(p); m = BDY(p); len = p->len;
2030: d = m0->dl; td = m0->td; c = C(m0);
2031: mr = BDY(r);
1.9 noro 2032: for ( ; len > 0; len--, NMV_ADV(m), NMV_ADV(mr) ) {
2033: c1 = C(m); DMAR(c1,c,0,nd_mod,c2); C(mr) = c2;
2034: mr->td = m->td+td; ndl_add(m->dl,d,mr->dl);
2035: }
2036: NV(r) = NV(p);
2037: r->len = p->len;
2038: r->sugar = p->sugar + td;
2039: }
2040: }
2041:
2042: ND ndv_mul_nm_create(NDV p,NM m0)
2043: {
2044: NM mr,mr0;
2045: NMV m;
2046: unsigned int *d,*dt,*dm;
2047: int c,n,td,i,c1,c2,len;
2048: ND r;
2049:
2050: if ( !p )
2051: return 0;
2052: else {
2053: n = NV(p); m = BDY(p);
2054: d = m0->dl; td = m0->td; c = C(m0);
2055: len = p->len;
2056: mr0 = 0;
2057: for ( i = 0; i < len; i++, NMV_ADV(m) ) {
2058: NEXTNM(mr0,mr);
1.4 noro 2059: c1 = C(m);
2060: DMAR(c1,c,0,nd_mod,c2);
2061: C(mr) = c2;
2062: mr->td = m->td+td;
2063: ndl_add(m->dl,d,mr->dl);
2064: }
1.9 noro 2065: NEXT(mr) = 0;
2066: MKND(NV(p),mr0,r);
1.4 noro 2067: r->sugar = p->sugar + td;
1.9 noro 2068: return r;
1.4 noro 2069: }
2070: }
2071:
2072: ND ndv_add(ND p1,NDV p2)
2073: {
1.9 noro 2074: register NM prev,cur,new;
1.4 noro 2075: int c,c1,c2,t,td,td2,mul,len,i;
1.9 noro 2076: NM head;
1.4 noro 2077: unsigned int *d;
2078: NMV m2;
2079:
2080: if ( !p1 )
2081: return 0;
2082: else {
2083: prev = 0; head = cur = BDY(p1);
1.9 noro 2084: NEWNM(new); len = p2->len;
2085: for ( m2 = BDY(p2), i = 0; cur && i < len; ) {
1.4 noro 2086: td2 = new->td = m2->td;
1.9 noro 2087: if ( cur->td > td2 ) c = 1;
2088: else if ( cur->td < td2 ) c = -1;
2089: else c = ndl_compare(cur->dl,m2->dl);
1.4 noro 2090: switch ( c ) {
2091: case 0:
2092: t = C(m2)+C(cur)-nd_mod;
1.9 noro 2093: if ( t < 0 ) t += nd_mod;
2094: if ( t ) C(cur) = t;
1.4 noro 2095: else if ( !prev ) {
1.9 noro 2096: head = NEXT(cur); FREENM(cur); cur = head;
1.4 noro 2097: } else {
1.9 noro 2098: NEXT(prev) = NEXT(cur); FREENM(cur); cur = NEXT(prev);
1.4 noro 2099: }
2100: NMV_ADV(m2); i++;
2101: break;
2102: case 1:
1.9 noro 2103: prev = cur; cur = NEXT(cur);
1.4 noro 2104: break;
2105: case -1:
1.6 noro 2106: ndl_copy(m2->dl,new->dl);
1.9 noro 2107: C(new) = C(m2);
1.4 noro 2108: if ( !prev ) {
2109: /* cur = head */
1.9 noro 2110: prev = new; NEXT(prev) = head; head = prev;
1.4 noro 2111: } else {
1.9 noro 2112: NEXT(prev) = new; NEXT(new) = cur; prev = new;
1.4 noro 2113: }
1.9 noro 2114: NEWNM(new); NMV_ADV(m2); i++;
1.4 noro 2115: break;
2116: }
2117: }
1.9 noro 2118: for ( ; i < len; i++, NMV_ADV(m2) ) {
2119: td2 = new->td = m2->td; C(new) = C(m2); ndl_copy(m2->dl,new->dl);
2120: if ( !prev ) {
2121: prev = new; NEXT(prev) = 0; head = prev;
2122: } else {
2123: NEXT(prev) = new; NEXT(new) = 0; prev = new;
2124: }
2125: NEWNM(new);
2126: }
1.4 noro 2127: FREENM(new);
2128: if ( head ) {
1.9 noro 2129: BDY(p1) = head; p1->sugar = MAX(p1->sugar,p2->sugar);
1.4 noro 2130: return p1;
2131: } else {
2132: FREEND(p1);
2133: return 0;
2134: }
2135:
2136: }
2137: }
2138:
1.11 ! noro 2139: void ndv_realloc(NDV p,int obpe,int oadv)
! 2140: {
! 2141: NMV m,mr,mr0;
! 2142: int len,i;
! 2143:
! 2144: #define NMV_OADV(m) (m = (NMV)(((char *)m)+oadv))
! 2145:
! 2146: if ( p ) {
! 2147: m = BDY(p); len = p->len;
! 2148: mr0 = mr = (NMV)MALLOC_ATOMIC(len*nmv_adv);
! 2149: bzero(mr0,len*nmv_adv);
! 2150: for ( i = 0; i < len; i++, NMV_OADV(m), NMV_ADV(mr) ) {
! 2151: C(mr) = C(m);
! 2152: mr->td = m->td;
! 2153: ndl_dup(obpe,m->dl,mr->dl);
! 2154: }
! 2155: BDY(p) = mr0;
! 2156: }
! 2157: }
! 2158:
1.3 noro 2159: NDV ndtondv(ND p)
2160: {
2161: NDV d;
2162: NMV m,m0;
2163: NM t;
2164: int i,len;
2165:
2166: if ( !p )
2167: return 0;
2168: len = nd_length(p);
2169: m0 = m = (NMV)MALLOC_ATOMIC(len*nmv_adv);
2170: for ( t = BDY(p), i = 0; t; t = NEXT(t), i++, NMV_ADV(m) ) {
2171: m->td = t->td;
1.7 noro 2172: ndl_copy(t->dl,m->dl);
1.3 noro 2173: m->c = t->c;
2174: }
2175: MKNDV(NV(p),m0,len,d);
2176: d->sugar = p->sugar;
2177: return d;
2178: }
2179:
1.11 ! noro 2180: NDV dptondv(DP p)
! 2181: {
! 2182: NDV d;
! 2183: NMV m0,m;
! 2184: MP t;
! 2185: int l,i,n;
! 2186:
! 2187: if ( !p )
! 2188: return 0;
! 2189: for ( t = BDY(p), l = 0; t; t = NEXT(t), l++ );
! 2190: m0 = m = (NMV)MALLOC_ATOMIC(l*nmv_adv);
! 2191: n = NV(p);
! 2192: for ( t = BDY(p), i = 0; i < l; i++, t = NEXT(t), NMV_ADV(m) ) {
! 2193: m->c = ITOS(t->c);
! 2194: m->td = t->dl->td;
! 2195: dltondl(n,t->dl,m->dl);
! 2196: }
! 2197: MKNDV(n,m0,l,d);
! 2198: d->sugar = p->sugar;
! 2199: return d;
! 2200: }
! 2201:
! 2202: DP ndvtodp(NDV p)
! 2203: {
! 2204: DP d;
! 2205: MP m0,m;
! 2206: NMV t;
! 2207: int len,i,n;
! 2208:
! 2209: if ( !p )
! 2210: return 0;
! 2211: m0 = 0;
! 2212: len = p->len;
! 2213: n = NV(p);
! 2214: for ( t = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) {
! 2215: NEXTMP(m0,m);
! 2216: m->c = STOI(t->c);
! 2217: m->dl = ndltodl(n,t->td,t->dl);
! 2218: }
! 2219: NEXT(m) = 0;
! 2220: MKDP(NV(p),m0,d);
! 2221: d->sugar = p->sugar;
! 2222: return d;
! 2223: }
! 2224:
1.3 noro 2225: void ndv_print(NDV p)
2226: {
2227: NMV m;
2228: int i,len;
2229:
2230: if ( !p )
2231: printf("0\n");
2232: else {
2233: len = p->len;
2234: for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
2235: printf("+%d*",m->c);
2236: ndl_print(m->dl);
2237: }
2238: printf("\n");
2239: }
1.11 ! noro 2240: }
! 2241: #else
! 2242: int nd_sp(ND_pairs p,ND *rp)
! 2243: {
! 2244: NM m;
! 2245: ND p1,p2;
! 2246: ND t1,t2;
! 2247: unsigned int *lcm;
! 2248: int td;
! 2249:
! 2250: p1 = nd_ps[p->i1];
! 2251: p2 = nd_ps[p->i2];
! 2252: lcm = p->lcm;
! 2253: td = p->td;
! 2254: NEWNM(m);
! 2255: C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl);
! 2256: if ( ndl_check_bound2(p->i1,m->dl) )
! 2257: return 0;
! 2258: t1 = nd_mul_ind_nm(p->i1,m);
! 2259: C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
! 2260: if ( ndl_check_bound2(p->i2,m->dl) ) {
! 2261: nd_free(t1);
! 2262: return 0;
! 2263: }
! 2264: t2 = nd_mul_ind_nm(p->i2,m);
! 2265: FREENM(m);
! 2266: *rp = nd_add(t1,t2);
! 2267: return 1;
! 2268: }
! 2269:
! 2270: ND nd_mul_nm(ND p,NM m0)
! 2271: {
! 2272: NM m,mr,mr0;
! 2273: unsigned int *d,*dt,*dm;
! 2274: int c,n,td,i,c1,c2;
! 2275: int *pt,*p1,*p2;
! 2276: ND r;
! 2277:
! 2278: if ( !p )
! 2279: return 0;
! 2280: else {
! 2281: n = NV(p); m = BDY(p);
! 2282: d = m0->dl; td = m0->td; c = C(m0);
! 2283: mr0 = 0;
! 2284: for ( ; m; m = NEXT(m) ) {
! 2285: NEXTNM(mr0,mr);
! 2286: c1 = C(m);
! 2287: DMAR(c1,c,0,nd_mod,c2);
! 2288: C(mr) = c2;
! 2289: mr->td = m->td+td;
! 2290: ndl_add(m->dl,d,mr->dl);
! 2291: }
! 2292: NEXT(mr) = 0;
! 2293: MKND(NV(p),mr0,r);
! 2294: r->sugar = p->sugar + td;
! 2295: return r;
! 2296: }
! 2297: }
! 2298:
! 2299: ND nd_mul_ind_nm(int index,NM m0)
! 2300: {
! 2301: register int c1,c2,c;
! 2302: register NM m,new,prev;
! 2303: NM mr0;
! 2304: unsigned int *d;
! 2305: int n,td,i,len,d0,d1;
! 2306: ND p,r;
! 2307:
! 2308: p = nd_ps[index];
! 2309: len = nd_psl[index];
! 2310: n = NV(p); m = BDY(p);
! 2311: d = m0->dl; td = m0->td; c = C(m0);
! 2312:
! 2313: NEWNM(mr0);
! 2314: c1 = C(m); DMAR(c1,c,0,nd_mod,c2); C(mr0) = c2;
! 2315: mr0->td = m->td+td; ndl_add(m->dl,d,mr0->dl);
! 2316: prev = mr0; m = NEXT(m);
! 2317: len--;
! 2318:
! 2319: switch ( nd_wpd ) {
! 2320: case 1:
! 2321: d0 = d[0];
! 2322: while ( len-- ) {
! 2323: c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
! 2324: NEWNM(new); C(new) = c2;
! 2325: new->td = m->td+td; new->dl[0] = m->dl[0]+d0;
! 2326: m = NEXT(m); NEXT(prev) = new; prev = new;
! 2327: }
! 2328: break;
! 2329: case 2:
! 2330: d0 = d[0]; d1 = d[1];
! 2331: while ( len-- ) {
! 2332: c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
! 2333: NEWNM(new); C(new) = c2;
! 2334: new->td = m->td+td;
! 2335: new->dl[0] = m->dl[0]+d0;
! 2336: new->dl[1] = m->dl[1]+d1;
! 2337: m = NEXT(m); NEXT(prev) = new; prev = new;
! 2338: }
! 2339: break;
! 2340: default:
! 2341: while ( len-- ) {
! 2342: c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
! 2343: NEWNM(new); C(new) = c2;
! 2344: new->td = m->td+td; ndl_add(m->dl,d,new->dl);
! 2345: m = NEXT(m); NEXT(prev) = new; prev = new;
! 2346: }
! 2347: break;
! 2348: }
! 2349:
! 2350: NEXT(prev) = 0;
! 2351: MKND(NV(p),mr0,r);
! 2352: r->sugar = p->sugar + td;
! 2353: return r;
1.3 noro 2354: }
1.4 noro 2355: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>