![[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.2 ! noro 1: /* $OpenXM$ */
! 2:
1.1 noro 3: #include "ca.h"
4: #include "inline.h"
5:
6: #if defined(__GNUC__)
7: #define INLINE inline
8: #elif defined(VISUAL)
9: #define INLINE __inline
10: #else
11: #define INLINE
12: #endif
13:
14: #define REDTAB_LEN 32003
15:
16: typedef struct oPGeoBucket {
17: int m;
18: struct oND *body[32];
19: } *PGeoBucket;
20:
21: typedef struct oND {
22: struct oNM *body;
23: int nv;
24: int sugar;
25: } *ND;
26:
27: typedef struct oNM {
28: struct oNM *next;
29: int td;
30: int c;
31: unsigned int dl[1];
32: } *NM;
33:
34: typedef struct oND_pairs {
35: struct oND_pairs *next;
36: int i1,i2;
37: int td,sugar;
38: unsigned int lcm[1];
39: } *ND_pairs;
40:
41: static ND *nd_ps;
42: static unsigned int **nd_bound;
43: int nd_mod,nd_nvar;
44: int is_rlex;
45: int nd_epw,nd_bpe,nd_wpd;
46: unsigned int nd_mask[32];
47: unsigned int nd_mask0,nd_mask1;
48: NM _nm_free_list;
49: ND _nd_free_list;
50: ND_pairs _ndp_free_list;
51: NM *nd_red;
52: int nd_red_len;
53:
54: extern int Top,Reverse;
55: int nd_psn,nd_pslen;
56: int nd_found,nd_create,nd_notfirst;
57:
58: void GC_gcollect();
59: NODE append_one(NODE,int);
60:
61: #define HTD(d) ((d)->body->td)
62: #define HDL(d) ((d)->body->dl)
63: #define HC(d) ((d)->body->c)
64:
65: #define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list)
66: #define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list)
67: #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)
68:
69: #define NEXTNM(r,c) \
70: if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);}
71: #define NEXTNM2(r,c,s) \
72: if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
73: #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)
74: #define FREENDP(m) NEXT(m)=_ndp_free_list; _ndp_free_list=(m)
75: #define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m)
76:
77: #define NEXTND_pairs(r,c) \
78: if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
79:
80: ND_pairs crit_B( ND_pairs d, int s );
81: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp);
82: NODE nd_setup(NODE f);
83: int nd_newps(ND a);
84: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest );
85: NODE update_base(NODE nd,int ndp);
86: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest );
87: int crit_2( int dp1, int dp2 );
88: ND_pairs crit_F( ND_pairs d1 );
89: ND_pairs crit_M( ND_pairs d1 );
90: ND_pairs nd_newpairs( NODE g, int t );
91: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t);
92: NODE nd_gb(NODE f);
93: void nd_free_private_storage();
94: void _NM_alloc();
95: void _ND_alloc();
96: int ndl_td(unsigned int *d);
97: ND nd_add(ND p1,ND p2);
98: ND nd_mul_nm(ND p,NM m0);
99: ND nd_mul_term(ND p,int td,unsigned int *d);
100: int nd_sp(ND_pairs p,ND *nf);
101: int nd_find_reducer(ND g,ND *red);
102: int nd_nf(ND g,int full,ND *nf);
103: ND nd_reduce(ND p1,ND p2);
104: ND nd_reduce_special(ND p1,ND p2);
105: void nd_free(ND p);
106: void ndl_print(unsigned int *dl);
107: void nd_print(ND p);
108: void ndp_print(ND_pairs d);
109: int nd_length(ND p);
110: void nd_monic(ND p);
111: void nd_mul_c(ND p,int mul);
112: void nd_free_redlist();
113: void nd_append_red(unsigned int *d,int td,int i);
114: unsigned int *nd_compute_bound(ND p);
115: ND_pairs nd_reconstruct(ND_pairs);
116: void nd_setup_parameters();
117: ND nd_dup(ND p,int obpe);
118: void ndl_dup(int obpe,unsigned int *d,unsigned int *r);
119:
120: void nd_free_private_storage()
121: {
122: _nd_free_list = 0;
123: _nm_free_list = 0;
124: nd_red = 0;
125: GC_gcollect();
126: }
127:
128: void _NM_alloc()
129: {
130: NM p;
131: int i;
132:
133: for ( i = 0; i < 16; i++ ) {
134: p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
135: p->next = _nm_free_list; _nm_free_list = p;
136: }
137: }
138:
139: void _ND_alloc()
140: {
141: ND p;
142: int i;
143:
144: for ( i = 0; i < 1024; i++ ) {
145: p = (ND)GC_malloc(sizeof(struct oND));
146: p->body = (NM)_nd_free_list; _nd_free_list = p;
147: }
148: }
149:
150: void _NDP_alloc()
151: {
152: ND_pairs p;
153: int i;
154:
155: for ( i = 0; i < 10240; i++ ) {
156: p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs)
157: +(nd_wpd-1)*sizeof(unsigned int));
158: p->next = _ndp_free_list; _ndp_free_list = p;
159: }
160: }
161:
162: INLINE nd_length(ND p)
163: {
164: NM m;
165: int i;
166:
167: if ( !p )
168: return 0;
169: else {
170: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
171: return i;
172: }
173: }
174:
175: INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2)
176: {
177: unsigned int u1,u2;
178: int i,j;
179:
180: switch ( nd_bpe ) {
181: case 4:
182: for ( i = 0; i < nd_wpd; i++ ) {
183: u1 = d1[i]; u2 = d2[i];
184: if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0;
185: if ( (u1&0xf000000) < (u2&0xf000000) ) return 0;
186: if ( (u1&0xf00000) < (u2&0xf00000) ) return 0;
187: if ( (u1&0xf0000) < (u2&0xf0000) ) return 0;
188: if ( (u1&0xf000) < (u2&0xf000) ) return 0;
189: if ( (u1&0xf00) < (u2&0xf00) ) return 0;
190: if ( (u1&0xf0) < (u2&0xf0) ) return 0;
191: if ( (u1&0xf) < (u2&0xf) ) return 0;
192: }
193: return 1;
194: break;
195: case 6:
196: for ( i = 0; i < nd_wpd; i++ ) {
197: u1 = d1[i]; u2 = d2[i];
198: if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0;
199: if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0;
200: if ( (u1&0x3f000) < (u2&0x3f000) ) return 0;
201: if ( (u1&0xfc0) < (u2&0xfc0) ) return 0;
202: if ( (u1&0x3f) < (u2&0x3f) ) return 0;
203: }
204: return 1;
205: break;
206: case 8:
207: for ( i = 0; i < nd_wpd; i++ ) {
208: u1 = d1[i]; u2 = d2[i];
209: if ( (u1&0xff000000) < (u2&0xff000000) ) return 0;
210: if ( (u1&0xff0000) < (u2&0xff0000) ) return 0;
211: if ( (u1&0xff00) < (u2&0xff00) ) return 0;
212: if ( (u1&0xff) < (u2&0xff) ) return 0;
213: }
214: return 1;
215: break;
216: case 16:
217: for ( i = 0; i < nd_wpd; i++ ) {
218: u1 = d1[i]; u2 = d2[i];
219: if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0;
220: if ( (u1&0xffff) < (u2&0xffff) ) return 0;
221: }
222: return 1;
223: break;
224: case 32:
225: for ( i = 0; i < nd_wpd; i++ )
226: if ( d1[i] < d2[i] ) return 0;
227: return 1;
228: break;
229: default:
230: for ( i = 0; i < nd_wpd; i++ ) {
231: u1 = d1[i]; u2 = d2[i];
232: for ( j = 0; j < nd_epw; j++ )
233: if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0;
234: }
235: return 1;
236: }
237: }
238:
239: void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
240: {
241: unsigned int t1,t2,u,u1,u2;
242: int i,j;
243:
244: switch ( nd_bpe ) {
245: case 4:
246: for ( i = 0; i < nd_wpd; i++ ) {
247: u1 = d1[i]; u2 = d2[i];
248: t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2;
249: t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2;
250: t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2;
251: t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2;
252: t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2;
253: t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2;
254: t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2;
255: t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2;
256: d[i] = u;
257: }
258: break;
259: case 6:
260: for ( i = 0; i < nd_wpd; i++ ) {
261: u1 = d1[i]; u2 = d2[i];
262: t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2;
263: t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2;
264: t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2;
265: t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2;
266: t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2;
267: d[i] = u;
268: }
269: break;
270: case 8:
271: for ( i = 0; i < nd_wpd; i++ ) {
272: u1 = d1[i]; u2 = d2[i];
273: t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2;
274: t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2;
275: t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2;
276: t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2;
277: d[i] = u;
278: }
279: break;
280: case 16:
281: for ( i = 0; i < nd_wpd; i++ ) {
282: u1 = d1[i]; u2 = d2[i];
283: t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2;
284: t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2;
285: d[i] = u;
286: }
287: break;
288: case 32:
289: for ( i = 0; i < nd_wpd; i++ ) {
290: u1 = d1[i]; u2 = d2[i];
291: d[i] = u1>u2?u1:u2;
292: }
293: break;
294: default:
295: for ( i = 0; i < nd_wpd; i++ ) {
296: u1 = d1[i]; u2 = d2[i];
297: for ( j = 0, u = 0; j < nd_epw; j++ ) {
298: t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2;
299: }
300: d[i] = u;
301: }
302: break;
303: }
304: }
305:
306: int ndl_td(unsigned int *d)
307: {
308: unsigned int t,u;
309: int i,j;
310:
311: for ( t = 0, i = 0; i < nd_wpd; i++ ) {
312: u = d[i];
313: for ( j = 0; j < nd_epw; j++, u>>=nd_bpe )
314: t += (u&nd_mask0);
315: }
316: return t;
317: }
318:
319: INLINE int ndl_compare(unsigned int *d1,unsigned int *d2)
320: {
321: int i;
322:
323: for ( i = 0; i < nd_wpd; i++, d1++, d2++ )
324: if ( *d1 > *d2 )
325: return is_rlex ? -1 : 1;
326: else if ( *d1 < *d2 )
327: return is_rlex ? 1 : -1;
328: return 0;
329: }
330:
331: INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)
332: {
333: int i;
334:
335: for ( i = 0; i < nd_wpd; i++ )
336: if ( d1[i] != d2[i] )
337: return 0;
338: return 1;
339: }
340:
341: INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d)
342: {
343: int i;
344:
345: for ( i = 0; i < nd_wpd; i++ ) {
346: d[i] = d1[i]+d2[i];
347: }
348: }
349:
350: void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
351: {
352: int i;
353:
354: for ( i = 0; i < nd_wpd; i++ )
355: d[i] = d1[i]-d2[i];
356: }
357:
358: int ndl_disjoint(unsigned int *d1,unsigned int *d2)
359: {
360: unsigned int t1,t2,u,u1,u2;
361: int i,j;
362:
363: switch ( nd_bpe ) {
364: case 4:
365: for ( i = 0; i < nd_wpd; i++ ) {
366: u1 = d1[i]; u2 = d2[i];
367: t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0;
368: t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0;
369: t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0;
370: t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0;
371: t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0;
372: t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0;
373: t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0;
374: t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0;
375: }
376: return 1;
377: break;
378: case 6:
379: for ( i = 0; i < nd_wpd; i++ ) {
380: u1 = d1[i]; u2 = d2[i];
381: t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0;
382: t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0;
383: t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0;
384: t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0;
385: t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0;
386: }
387: return 1;
388: break;
389: case 8:
390: for ( i = 0; i < nd_wpd; i++ ) {
391: u1 = d1[i]; u2 = d2[i];
392: t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0;
393: t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0;
394: t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0;
395: t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0;
396: }
397: return 1;
398: break;
399: case 16:
400: for ( i = 0; i < nd_wpd; i++ ) {
401: u1 = d1[i]; u2 = d2[i];
402: t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0;
403: t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0;
404: }
405: return 1;
406: break;
407: case 32:
408: for ( i = 0; i < nd_wpd; i++ )
409: if ( d1[i] && d2[i] ) return 0;
410: return 1;
411: break;
412: default:
413: for ( i = 0; i < nd_wpd; i++ ) {
414: u1 = d1[i]; u2 = d2[i];
415: for ( j = 0; j < nd_epw; j++ ) {
416: if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0;
417: u1 >>= nd_bpe; u2 >>= nd_bpe;
418: }
419: }
420: return 1;
421: break;
422: }
423: }
424:
425: ND nd_reduce(ND p1,ND p2)
426: {
427: int c,c1,c2,t,td,td2,mul;
428: NM m2,prev,head,cur,new;
429: unsigned int *d;
430:
431: if ( !p1 )
432: return 0;
433: else {
434: c2 = invm(HC(p2),nd_mod);
435: c1 = nd_mod-HC(p1);
436: DMAR(c1,c2,0,nd_mod,mul);
437: td = HTD(p1)-HTD(p2);
438: d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
439: ndl_sub(HDL(p1),HDL(p2),d);
440: prev = 0; head = cur = BDY(p1);
441: NEWNM(new);
442: for ( m2 = BDY(p2); m2; ) {
443: td2 = new->td = m2->td+td;
444: ndl_add(m2->dl,d,new->dl);
445: if ( !cur ) {
446: c1 = C(m2);
447: DMAR(c1,mul,0,nd_mod,c2);
448: C(new) = c2;
449: if ( !prev ) {
450: prev = new;
451: NEXT(prev) = 0;
452: head = prev;
453: } else {
454: NEXT(prev) = new;
455: NEXT(new) = 0;
456: prev = new;
457: }
458: m2 = NEXT(m2);
459: NEWNM(new);
460: continue;
461: }
462: if ( cur->td > td2 )
463: c = 1;
464: else if ( cur->td < td2 )
465: c = -1;
466: else
467: c = ndl_compare(cur->dl,new->dl);
468: switch ( c ) {
469: case 0:
470: c2 = C(m2);
471: c1 = C(cur);
472: DMAR(c2,mul,c1,nd_mod,t);
473: if ( t )
474: C(cur) = t;
475: else if ( !prev ) {
476: head = NEXT(cur);
477: FREENM(cur);
478: cur = head;
479: } else {
480: NEXT(prev) = NEXT(cur);
481: FREENM(cur);
482: cur = NEXT(prev);
483: }
484: m2 = NEXT(m2);
485: break;
486: case 1:
487: prev = cur;
488: cur = NEXT(cur);
489: break;
490: case -1:
491: if ( !prev ) {
492: /* cur = head */
493: prev = new;
494: c2 = C(m2);
495: DMAR(c2,mul,0,nd_mod,c1);
496: C(prev) = c1;
497: NEXT(prev) = head;
498: head = prev;
499: } else {
500: c2 = C(m2);
501: DMAR(c2,mul,0,nd_mod,c1);
502: C(new) = c1;
503: NEXT(prev) = new;
504: NEXT(new) = cur;
505: prev = new;
506: }
507: NEWNM(new);
508: m2 = NEXT(m2);
509: break;
510: }
511: }
512: FREENM(new);
513: if ( head ) {
514: BDY(p1) = head;
515: p1->sugar = MAX(p1->sugar,p2->sugar+td);
516: return p1;
517: } else {
518: FREEND(p1);
519: return 0;
520: }
521:
522: }
523: }
524:
525: /* HDL(p1) = HDL(p2) */
526:
527: ND nd_reduce_special(ND p1,ND p2)
528: {
529: int c,c1,c2,t,td,td2,mul;
530: NM m2,prev,head,cur,new;
531:
532: if ( !p1 )
533: return 0;
534: else {
535: c2 = invm(HC(p2),nd_mod);
536: c1 = nd_mod-HC(p1);
537: DMAR(c1,c2,0,nd_mod,mul);
538: prev = 0; head = cur = BDY(p1);
539: NEWNM(new);
540: for ( m2 = BDY(p2); m2; ) {
541: td2 = new->td = m2->td;
542: if ( !cur ) {
543: c1 = C(m2);
544: DMAR(c1,mul,0,nd_mod,c2);
545: C(new) = c2;
546: bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
547: if ( !prev ) {
548: prev = new;
549: NEXT(prev) = 0;
550: head = prev;
551: } else {
552: NEXT(prev) = new;
553: NEXT(new) = 0;
554: prev = new;
555: }
556: m2 = NEXT(m2);
557: NEWNM(new);
558: continue;
559: }
560: if ( cur->td > td2 )
561: c = 1;
562: else if ( cur->td < td2 )
563: c = -1;
564: else
565: c = ndl_compare(cur->dl,m2->dl);
566: switch ( c ) {
567: case 0:
568: c2 = C(m2);
569: c1 = C(cur);
570: DMAR(c2,mul,c1,nd_mod,t);
571: if ( t )
572: C(cur) = t;
573: else if ( !prev ) {
574: head = NEXT(cur);
575: FREENM(cur);
576: cur = head;
577: } else {
578: NEXT(prev) = NEXT(cur);
579: FREENM(cur);
580: cur = NEXT(prev);
581: }
582: m2 = NEXT(m2);
583: break;
584: case 1:
585: prev = cur;
586: cur = NEXT(cur);
587: break;
588: case -1:
589: bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
590: if ( !prev ) {
591: /* cur = head */
592: prev = new;
593: c2 = C(m2);
594: DMAR(c2,mul,0,nd_mod,c1);
595: C(prev) = c1;
596: NEXT(prev) = head;
597: head = prev;
598: } else {
599: c2 = C(m2);
600: DMAR(c2,mul,0,nd_mod,c1);
601: C(new) = c1;
602: NEXT(prev) = new;
603: NEXT(new) = cur;
604: prev = new;
605: }
606: NEWNM(new);
607: m2 = NEXT(m2);
608: break;
609: }
610: }
611: FREENM(new);
612: if ( head ) {
613: BDY(p1) = head;
614: p1->sugar = MAX(p1->sugar,p2->sugar+td);
615: return p1;
616: } else {
617: FREEND(p1);
618: return 0;
619: }
620:
621: }
622: }
623:
624: INLINE int ndl_check_bound(unsigned int *d)
625: {
626: int i;
627:
628: for ( i = 0; i < nd_wpd; i++ )
629: if ( d[i] & nd_mask1 )
630: return 1;
631: return 0;
632: }
633:
634: int nd_sp(ND_pairs p,ND *rp)
635: {
636: NM m;
637: ND p1,p2,t1,t2;
638: unsigned int *lcm,*check;
639: int td;
640:
641: check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
642: p1 = nd_ps[p->i1];
643: p2 = nd_ps[p->i2];
644: lcm = p->lcm;
645: td = p->td;
646: NEWNM(m);
647: C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl); NEXT(m) = 0;
648: ndl_add(nd_bound[p->i1],m->dl,check);
649: if ( ndl_check_bound(check) )
650: return 0;
651: t1 = nd_mul_nm(p1,m);
652: C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
653: ndl_add(nd_bound[p->i2],m->dl,check);
654: if ( ndl_check_bound(check) ) {
655: nd_free(t1);
656: return 0;
657: }
658: t2 = nd_mul_nm(p2,m);
659: FREENM(m);
660: *rp = nd_add(t1,t2);
661: return 1;
662: }
663:
664: int ndl_hash_value(int td,unsigned int *d)
665: {
666: int i;
667: int r;
668:
669: r = td;
670: for ( i = 0; i < nd_wpd; i++ )
671: r = ((r<<16)+d[i])%REDTAB_LEN;
672: return r;
673: }
674:
675: int nd_find_reducer(ND g, ND *rp)
676: {
677: NM m;
678: ND r,p;
679: int i,c1,c2,c;
680: int d,k,append,index;
681: unsigned int *check;
682: NM t;
683:
684: d = ndl_hash_value(HTD(g),HDL(g));
685: for ( m = nd_red[d], k = 0; m; m = NEXT(m), k++ ) {
686: if ( HTD(g) == m->td && ndl_equal(HDL(g),m->dl) ) {
687: if ( k > 0 ) nd_notfirst++;
688: index = m->c;
689: append = 0;
690: nd_found++;
691: goto found;
692: }
693: }
694:
695: for ( i = 0; i < nd_psn; i++ ) {
696: p = nd_ps[i];
697: if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
698: index = i;
699: append = 1;
700: nd_create++;
701: goto found;
702: }
703: }
704: return 0;
705:
706: found:
707: NEWNM(m);
708: p = nd_ps[index];
709: ndl_sub(HDL(g),HDL(p),m->dl);
710:
711: check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
712: ndl_add(nd_bound[index],m->dl,check);
713: if ( ndl_check_bound(check) ) {
714: FREENM(m);
715: return -1;
716: }
717:
718: c1 = invm(HC(p),nd_mod);
719: c2 = nd_mod-HC(g);
720: DMAR(c1,c2,0,nd_mod,c);
721: C(m) = c;
722: m->td = HTD(g)-HTD(p);
723: NEXT(m) = 0;
724: *rp = r = nd_mul_nm(p,m);
725: FREENM(m);
726:
727: if ( append ) nd_append_red(HDL(g),HTD(g),i);
728: return 1;
729: }
730:
731: ND nd_find_monic_reducer(ND g)
732: {
733: int *d;
734: ND p,r;
735: int i;
736:
737: for ( i = 0; i < nd_psn; i++ ) {
738: p = nd_ps[i];
739: if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
740: d = (int *)ALLOCA(nd_wpd*sizeof(int));
741: ndl_sub(HDL(g),HDL(p),d);
742: r = nd_mul_term(p,HTD(g)-HTD(p),d);
743: return r;
744: }
745: }
746: return 0;
747: }
748:
749: ND nd_add(ND p1,ND p2)
750: {
751: int n,c;
752: int t;
753: ND r;
754: NM m1,m2,mr0,mr,s;
755:
756: if ( !p1 )
757: return p2;
758: else if ( !p2 )
759: return p1;
760: else {
761: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
762: if ( m1->td > m2->td )
763: c = 1;
764: else if ( m1->td < m2->td )
765: c = -1;
766: else
767: c = ndl_compare(m1->dl,m2->dl);
768: switch ( c ) {
769: case 0:
770: t = ((C(m1))+(C(m2))) - nd_mod;
771: if ( t < 0 )
772: t += nd_mod;
773: s = m1; m1 = NEXT(m1);
774: if ( t ) {
775: NEXTNM2(mr0,mr,s); C(mr) = (t);
776: } else {
777: FREENM(s);
778: }
779: s = m2; m2 = NEXT(m2); FREENM(s);
780: break;
781: case 1:
782: s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
783: break;
784: case -1:
785: s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
786: break;
787: }
788: }
789: if ( !mr0 )
790: if ( m1 )
791: mr0 = m1;
792: else if ( m2 )
793: mr0 = m2;
794: else
795: return 0;
796: else if ( m1 )
797: NEXT(mr) = m1;
798: else if ( m2 )
799: NEXT(mr) = m2;
800: else
801: NEXT(mr) = 0;
802: BDY(p1) = mr0;
803: p1->sugar = MAX(p1->sugar,p2->sugar);
804: FREEND(p2);
805: return p1;
806: }
807: }
808:
809: ND nd_mul_nm(ND p,NM m0)
810: {
811: NM m,mr,mr0;
812: unsigned int *d,*dt,*dm;
813: int c,n,td,i,c1,c2;
814: int *pt,*p1,*p2;
815: ND r;
816:
817: if ( !p )
818: return 0;
819: else {
820: n = NV(p); m = BDY(p);
821: d = m0->dl; td = m0->td; c = C(m0);
822: mr0 = 0;
823: for ( ; m; m = NEXT(m) ) {
824: NEXTNM(mr0,mr);
825: c1 = C(m);
826: DMAR(c1,c,0,nd_mod,c2);
827: C(mr) = c2;
828: mr->td = m->td+td;
829: ndl_add(m->dl,d,mr->dl);
830: }
831: NEXT(mr) = 0;
832: MKND(NV(p),mr0,r);
833: r->sugar = p->sugar + td;
834: return r;
835: }
836: }
837:
838: ND nd_mul_term(ND p,int td,unsigned int *d)
839: {
840: NM m,mr,mr0;
841: int c,n;
842: ND r;
843:
844: if ( !p )
845: return 0;
846: else {
847: n = NV(p); m = BDY(p);
848: for ( mr0 = 0; m; m = NEXT(m) ) {
849: NEXTNM(mr0,mr);
850: C(mr) = C(m);
851: mr->td = m->td+td;
852: ndl_add(m->dl,d,mr->dl);
853: }
854: NEXT(mr) = 0;
855: MKND(NV(p),mr0,r);
856: r->sugar = p->sugar + td;
857: return r;
858: }
859: }
860:
861: #if 1
862: /* ret=1 : success, ret=0 : overflow */
863: int nd_nf(ND g,int full,ND *rp)
864: {
865: ND p,d,red;
866: NM m,mrd,tail;
867: int n,sugar,psugar,stat;
868:
869: if ( !g ) {
870: *rp = 0;
871: return 1;
872: }
873: sugar = g->sugar;
874: n = NV(g);
875: for ( d = 0; g; ) {
876: /* stat=1 : found, stat=0 : not found, stat=-1 : overflow */
877: stat = nd_find_reducer(g,&red);
878: if ( stat == 1 ) {
879: #if 1
880: g = nd_add(g,red);
881: sugar = MAX(sugar,red->sugar);
882: #else
883: psugar = (HTD(g)-HTD(red))+red->sugar;
884: g = nd_reduce(g,red);
885: sugar = MAX(sugar,psugar);
886: #endif
887: } else if ( stat == -1 ) {
888: nd_free(g);
889: nd_free(d);
890: return 0;
891: } else if ( !full ) {
892: *rp = g;
893: return 1;
894: } else {
895: m = BDY(g);
896: if ( NEXT(m) ) {
897: BDY(g) = NEXT(m); NEXT(m) = 0;
898: } else {
899: FREEND(g); g = 0;
900: }
901: if ( d ) {
902: NEXT(tail)=m;
903: tail=m;
904: } else {
905: MKND(n,m,d);
906: tail = BDY(d);
907: }
908: }
909: }
910: if ( d )
911: d->sugar = sugar;
912: *rp = d;
913: return 1;
914: }
915: #else
916:
917: ND nd_remove_head(ND p)
918: {
919: NM m;
920:
921: m = BDY(p);
922: if ( !NEXT(m) ) {
923: FREEND(p);
924: p = 0;
925: } else
926: BDY(p) = NEXT(m);
927: FREENM(m);
928: return p;
929: }
930:
931: PGeoBucket create_pbucket()
932: {
933: PGeoBucket g;
934:
935: g = CALLOC(1,sizeof(struct oPGeoBucket));
936: g->m = -1;
937: return g;
938: }
939:
940: void add_pbucket(PGeoBucket g,ND d)
941: {
942: int l,k,m;
943:
944: l = nd_length(d);
945: for ( k = 0, m = 1; l > m; k++, m <<= 2 );
946: /* 4^(k-1) < l <= 4^k */
947: d = nd_add(g->body[k],d);
948: for ( ; d && nd_length(d) > 1<<(2*k); k++ ) {
949: g->body[k] = 0;
950: d = nd_add(g->body[k+1],d);
951: }
952: g->body[k] = d;
953: g->m = MAX(g->m,k);
954: }
955:
956: int head_pbucket(PGeoBucket g)
957: {
958: int j,i,c,k,nv,sum;
959: unsigned int *di,*dj;
960: ND gi,gj;
961:
962: k = g->m;
963: while ( 1 ) {
964: j = -1;
965: for ( i = 0; i <= k; i++ ) {
966: if ( !(gi = g->body[i]) )
967: continue;
968: if ( j < 0 ) {
969: j = i;
970: gj = g->body[j];
971: dj = HDL(gj);
972: sum = HC(gj);
973: } else {
974: di = HDL(gi);
975: nv = NV(gi);
976: if ( HTD(gi) > HTD(gj) )
977: c = 1;
978: else if ( HTD(gi) < HTD(gj) )
979: c = -1;
980: else
981: c = ndl_compare(di,dj);
982: if ( c > 0 ) {
983: if ( sum )
984: HC(gj) = sum;
985: else
986: g->body[j] = nd_remove_head(gj);
987: j = i;
988: gj = g->body[j];
989: dj = HDL(gj);
990: sum = HC(gj);
991: } else if ( c == 0 ) {
992: sum = sum+HC(gi)-nd_mod;
993: if ( sum < 0 )
994: sum += nd_mod;
995: g->body[i] = nd_remove_head(gi);
996: }
997: }
998: }
999: if ( j < 0 )
1000: return -1;
1001: else if ( sum ) {
1002: HC(gj) = sum;
1003: return j;
1004: } else
1005: g->body[j] = nd_remove_head(gj);
1006: }
1007: }
1008:
1009: ND normalize_pbucket(PGeoBucket g)
1010: {
1011: int i;
1012: ND r,t;
1013:
1014: r = 0;
1015: for ( i = 0; i <= g->m; i++ )
1016: r = nd_add(r,g->body[i]);
1017: return r;
1018: }
1019:
1020: ND nd_nf(ND g,int full)
1021: {
1022: ND u,p,d,red;
1023: NODE l;
1024: NM m,mrd;
1025: int sugar,psugar,n,h_reducible,h;
1026: PGeoBucket bucket;
1027:
1028: if ( !g ) {
1029: return 0;
1030: }
1031: sugar = g->sugar;
1032: n = g->nv;
1033: bucket = create_pbucket();
1034: add_pbucket(bucket,g);
1035: d = 0;
1036: while ( 1 ) {
1037: h = head_pbucket(bucket);
1038: if ( h < 0 ) {
1039: if ( d )
1040: d->sugar = sugar;
1041: return d;
1042: }
1043: g = bucket->body[h];
1044: red = nd_find_reducer(g);
1045: if ( red ) {
1046: bucket->body[h] = nd_remove_head(g);
1047: red = nd_remove_head(red);
1048: add_pbucket(bucket,red);
1049: sugar = MAX(sugar,red->sugar);
1050: } else if ( !full ) {
1051: g = normalize_pbucket(bucket);
1052: if ( g )
1053: g->sugar = sugar;
1054: return g;
1055: } else {
1056: m = BDY(g);
1057: if ( NEXT(m) ) {
1058: BDY(g) = NEXT(m); NEXT(m) = 0;
1059: } else {
1060: FREEND(g); g = 0;
1061: }
1062: bucket->body[h] = g;
1063: NEXT(m) = 0;
1064: if ( d ) {
1065: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1066: NEXT(mrd) = m;
1067: } else {
1068: MKND(n,m,d);
1069: }
1070: }
1071: }
1072: }
1073: #endif
1074:
1075: NODE nd_gb(NODE f)
1076: {
1077: int i,nh,sugar,stat;
1078: NODE r,g,gall;
1079: ND_pairs d;
1080: ND_pairs l;
1081: ND h,nf;
1082:
1083: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1084: i = (int)BDY(r);
1085: d = update_pairs(d,g,i);
1086: g = update_base(g,i);
1087: gall = append_one(gall,i);
1088: }
1089: sugar = 0;
1090: while ( d ) {
1091: again:
1092: l = nd_minp(d,&d);
1093: if ( l->sugar != sugar ) {
1094: sugar = l->sugar;
1095: fprintf(asir_out,"%d",sugar);
1096: }
1097: stat = nd_sp(l,&h);
1098: if ( !stat ) {
1099: NEXT(l) = d; d = l;
1100: d = nd_reconstruct(d);
1101: goto again;
1102: }
1103: stat = nd_nf(h,!Top,&nf);
1104: if ( !stat ) {
1105: NEXT(l) = d; d = l;
1106: d = nd_reconstruct(d);
1107: goto again;
1108: } else if ( nf ) {
1109: printf("+"); fflush(stdout);
1110: nh = nd_newps(nf);
1111: d = update_pairs(d,g,nh);
1112: g = update_base(g,nh);
1113: gall = append_one(gall,nh);
1114: FREENDP(l);
1115: } else {
1116: printf("."); fflush(stdout);
1117: FREENDP(l);
1118: }
1119: }
1120: return g;
1121: }
1122:
1123: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t)
1124: {
1125: ND_pairs d1,nd,cur,head,prev,remove;
1126:
1127: if ( !g ) return d;
1128: d = crit_B(d,t);
1129: d1 = nd_newpairs(g,t);
1130: d1 = crit_M(d1);
1131: d1 = crit_F(d1);
1132: prev = 0; cur = head = d1;
1133: while ( cur ) {
1134: if ( crit_2( cur->i1,cur->i2 ) ) {
1135: remove = cur;
1136: if ( !prev ) {
1137: head = cur = NEXT(cur);
1138: } else {
1139: cur = NEXT(prev) = NEXT(cur);
1140: }
1141: FREENDP(remove);
1142: } else {
1143: prev = cur;
1144: cur = NEXT(cur);
1145: }
1146: }
1147: if ( !d )
1148: return head;
1149: else {
1150: nd = d;
1151: while ( NEXT(nd) )
1152: nd = NEXT(nd);
1153: NEXT(nd) = head;
1154: return d;
1155: }
1156: }
1157:
1158: ND_pairs nd_newpairs( NODE g, int t )
1159: {
1160: NODE h;
1161: unsigned int *dl;
1162: int td,ts,s;
1163: ND_pairs r,r0;
1164:
1165: dl = HDL(nd_ps[t]);
1166: td = HTD(nd_ps[t]);
1167: ts = nd_ps[t]->sugar - td;
1168: for ( r0 = 0, h = g; h; h = NEXT(h) ) {
1169: NEXTND_pairs(r0,r);
1170: r->i1 = (int)BDY(h);
1171: r->i2 = t;
1172: ndl_lcm(HDL(nd_ps[r->i1]),dl,r->lcm);
1173: r->td = ndl_td(r->lcm);
1174: s = nd_ps[r->i1]->sugar-HTD(nd_ps[r->i1]);
1175: r->sugar = MAX(s,ts) + r->td;
1176: }
1177: NEXT(r) = 0;
1178: return r0;
1179: }
1180:
1181: ND_pairs crit_B( ND_pairs d, int s )
1182: {
1183: ND_pairs cur,head,prev,remove;
1184: unsigned int *t,*tl,*lcm;
1185: int td,tdl;
1186:
1187: if ( !d ) return 0;
1188: t = HDL(nd_ps[s]);
1189: prev = 0;
1190: head = cur = d;
1191: lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1192: while ( cur ) {
1193: tl = cur->lcm;
1194: if ( ndl_reducible(tl,t)
1195: && (ndl_lcm(HDL(nd_ps[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
1196: && (ndl_lcm(HDL(nd_ps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
1197: remove = cur;
1198: if ( !prev ) {
1199: head = cur = NEXT(cur);
1200: } else {
1201: cur = NEXT(prev) = NEXT(cur);
1202: }
1203: FREENDP(remove);
1204: } else {
1205: prev = cur;
1206: cur = NEXT(cur);
1207: }
1208: }
1209: return head;
1210: }
1211:
1212: ND_pairs crit_M( ND_pairs d1 )
1213: {
1214: ND_pairs e,d2,d3,dd,p;
1215: unsigned int *id,*jd;
1216: int itd,jtd;
1217:
1218: for ( dd = 0, e = d1; e; e = d3 ) {
1219: if ( !(d2 = NEXT(e)) ) {
1220: NEXT(e) = dd;
1221: return e;
1222: }
1223: id = e->lcm;
1224: itd = e->td;
1225: for ( d3 = 0; d2; d2 = p ) {
1226: p = NEXT(d2),
1227: jd = d2->lcm;
1228: jtd = d2->td;
1229: if ( jtd == itd )
1230: if ( id == jd );
1231: else if ( ndl_reducible(jd,id) ) continue;
1232: else if ( ndl_reducible(id,jd) ) goto delit;
1233: else ;
1234: else if ( jtd > itd )
1235: if ( ndl_reducible(jd,id) ) continue;
1236: else ;
1237: else if ( ndl_reducible(id,jd ) ) goto delit;
1238: NEXT(d2) = d3;
1239: d3 = d2;
1240: }
1241: NEXT(e) = dd;
1242: dd = e;
1243: continue;
1244: /**/
1245: delit: NEXT(d2) = d3;
1246: d3 = d2;
1247: for ( ; p; p = d2 ) {
1248: d2 = NEXT(p);
1249: NEXT(p) = d3;
1250: d3 = p;
1251: }
1252: FREENDP(e);
1253: }
1254: return dd;
1255: }
1256:
1257: ND_pairs crit_F( ND_pairs d1 )
1258: {
1259: ND_pairs rest, head,remove;
1260: ND_pairs last, p, r, w;
1261: int s;
1262:
1263: for ( head = last = 0, p = d1; NEXT(p); ) {
1264: r = w = equivalent_pairs(p,&rest);
1265: s = r->sugar;
1266: w = NEXT(w);
1267: while ( w ) {
1268: if ( crit_2(w->i1,w->i2) ) {
1269: r = w;
1270: w = NEXT(w);
1271: while ( w ) {
1272: remove = w;
1273: w = NEXT(w);
1274: FREENDP(remove);
1275: }
1276: break;
1277: } else if ( w->sugar < s ) {
1278: FREENDP(r);
1279: r = w;
1280: s = r->sugar;
1281: w = NEXT(w);
1282: } else {
1283: remove = w;
1284: w = NEXT(w);
1285: FREENDP(remove);
1286: }
1287: }
1288: if ( last ) NEXT(last) = r;
1289: else head = r;
1290: NEXT(last = r) = 0;
1291: p = rest;
1292: if ( !p ) return head;
1293: }
1294: if ( !last ) return p;
1295: NEXT(last) = p;
1296: return head;
1297: }
1298:
1299: int crit_2( int dp1, int dp2 )
1300: {
1301: return ndl_disjoint(HDL(nd_ps[dp1]),HDL(nd_ps[dp2]));
1302: }
1303:
1304: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest )
1305: {
1306: ND_pairs w,p,r,s;
1307: unsigned int *d;
1308: int td;
1309:
1310: w = d1;
1311: d = w->lcm;
1312: td = w->td;
1313: s = NEXT(w);
1314: NEXT(w) = 0;
1315: for ( r = 0; s; s = p ) {
1316: p = NEXT(s);
1317: if ( td == s->td && ndl_equal(d,s->lcm) ) {
1318: NEXT(s) = w;
1319: w = s;
1320: } else {
1321: NEXT(s) = r;
1322: r = s;
1323: }
1324: }
1325: *prest = r;
1326: return w;
1327: }
1328:
1329: NODE update_base(NODE nd,int ndp)
1330: {
1331: unsigned int *dl, *dln;
1332: NODE last, p, head;
1333: int td,tdn;
1334:
1335: dl = HDL(nd_ps[ndp]);
1336: td = HTD(nd_ps[ndp]);
1337: for ( head = last = 0, p = nd; p; ) {
1338: dln = HDL(nd_ps[(int)BDY(p)]);
1339: tdn = HTD(nd_ps[(int)BDY(p)]);
1340: if ( tdn >= td && ndl_reducible( dln, dl ) ) {
1341: p = NEXT(p);
1342: if ( last ) NEXT(last) = p;
1343: } else {
1344: if ( !last ) head = p;
1345: p = NEXT(last = p);
1346: }
1347: }
1348: head = append_one(head,ndp);
1349: return head;
1350: }
1351:
1352: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
1353: {
1354: ND_pairs m,ml,p,l;
1355: unsigned int *lcm;
1356: int s,td,len,tlen,c;
1357:
1358: if ( !(p = NEXT(m = d)) ) {
1359: *prest = p;
1360: NEXT(m) = 0;
1361: return m;
1362: }
1363: lcm = m->lcm;
1364: s = m->sugar;
1365: td = m->td;
1366: len = nd_length(nd_ps[m->i1])+nd_length(nd_ps[m->i2]);
1367: for ( ml = 0, l = m; p; p = NEXT(l = p) ) {
1368: if (p->sugar < s)
1369: goto find;
1370: else if ( p->sugar == s ) {
1371: if ( p->td < td )
1372: goto find;
1373: else if ( p->td == td ) {
1374: c = ndl_compare(p->lcm,lcm);
1375: if ( c < 0 )
1376: goto find;
1377: else if ( c == 0 ) {
1378: tlen = nd_length(nd_ps[p->i1])+nd_length(nd_ps[p->i2]);
1379: if ( tlen < len )
1380: goto find;
1381: }
1382: }
1383: }
1384: continue;
1385: find:
1386: ml = l;
1387: m = p;
1388: lcm = m->lcm;
1389: s = m->sugar;
1390: td = m->td;
1391: len = tlen;
1392: }
1393: if ( !ml ) *prest = NEXT(m);
1394: else {
1395: NEXT(ml) = NEXT(m);
1396: *prest = d;
1397: }
1398: NEXT(m) = 0;
1399: return m;
1400: }
1401:
1402: int nd_newps(ND a)
1403: {
1404: if ( nd_psn == nd_pslen ) {
1405: nd_pslen *= 2;
1406: nd_ps = (ND *)REALLOC((char *)nd_ps,nd_pslen*sizeof(ND));
1407: nd_bound = (unsigned int **)
1408: REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
1409: }
1410: nd_monic(a);
1411: nd_ps[nd_psn] = a;
1412: nd_bound[nd_psn] = nd_compute_bound(a);
1413: return nd_psn++;
1414: }
1415:
1416: NODE NODE_sortb(NODE f,int);
1417: ND dptond(DP);
1418: DP ndtodp(ND);
1419:
1420: NODE nd_setup(NODE f)
1421: {
1422: int i,td;
1423: NODE s,s0,f0;
1424:
1425: nd_found = 0;
1426: nd_notfirst = 0;
1427: nd_create = 0;
1428: #if 0
1429: f0 = f = NODE_sortb(f,1);
1430: #endif
1431: nd_psn = length(f); nd_pslen = 2*nd_psn;
1432: nd_ps = (ND *)MALLOC(nd_pslen*sizeof(ND));
1433: nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
1434: nd_bpe = 4;
1435: nd_setup_parameters();
1436: nd_free_private_storage();
1437: for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
1438: nd_ps[i] = dptond((DP)BDY(f));
1439: nd_monic(nd_ps[i]);
1440: nd_bound[i] = nd_compute_bound(nd_ps[i]);
1441: }
1442: nd_red = (NM *)MALLOC(REDTAB_LEN*sizeof(NM));
1443: for ( s0 = 0, i = 0; i < nd_psn; i++ ) {
1444: NEXTNODE(s0,s); BDY(s) = (pointer)i;
1445: }
1446: if ( s0 ) NEXT(s) = 0;
1447: return s0;
1448: }
1449:
1450: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
1451: {
1452: struct order_spec ord1;
1453: VL fv,vv,vc;
1454: NODE fd,fd0,r,r0,t,x,s,xx;
1455: DP a,b,c;
1456:
1457: get_vars((Obj)f,&fv); pltovl(v,&vv);
1458: nd_nvar = length(vv);
1459: if ( ord->id )
1460: error("nd_gr : unsupported order");
1461: switch ( ord->ord.simple ) {
1462: case 0:
1463: is_rlex = 1;
1464: break;
1465: case 1:
1466: is_rlex = 0;
1467: break;
1468: default:
1469: error("nd_gr : unsupported order");
1470: }
1471: initd(ord);
1472: nd_mod = m;
1473: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
1474: ptod(CO,vv,(P)BDY(t),&b);
1475: _dp_mod(b,m,0,&c);
1476: if ( c ) {
1477: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
1478: }
1479: }
1480: if ( fd0 ) NEXT(fd) = 0;
1481: s = nd_setup(fd0);
1482: x = nd_gb(s);
1483: #if 0
1484: x = nd_reduceall(x,m);
1485: #endif
1486: for ( r0 = 0; x; x = NEXT(x) ) {
1487: NEXTNODE(r0,r);
1488: a = ndtodp(nd_ps[(int)BDY(x)]);
1489: _dtop_mod(CO,vv,a,(P *)&BDY(r));
1490: }
1491: if ( r0 ) NEXT(r) = 0;
1492: MKLIST(*rp,r0);
1493: fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
1494: nd_found,nd_notfirst,nd_create);
1495: }
1496:
1497: void dltondl(int n,DL dl,unsigned int *r)
1498: {
1499: unsigned int *d;
1500: int i;
1501:
1502: d = dl->d;
1503: bzero(r,nd_wpd*sizeof(unsigned int));
1504: if ( is_rlex )
1505: for ( i = 0; i < n; i++ )
1506: r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe));
1507: else
1508: for ( i = 0; i < n; i++ )
1509: r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe);
1510: }
1511:
1512: DL ndltodl(int n,int td,unsigned int *ndl)
1513: {
1514: DL dl;
1515: int *d;
1516: int i;
1517:
1518: NEWDL(dl,n);
1519: dl->td = td;
1520: d = dl->d;
1521: if ( is_rlex )
1522: for ( i = 0; i < n; i++ )
1523: d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
1524: &((1<<nd_bpe)-1);
1525: else
1526: for ( i = 0; i < n; i++ )
1527: d[i] = (ndl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
1528: &((1<<nd_bpe)-1);
1529: return dl;
1530: }
1531:
1532: ND dptond(DP p)
1533: {
1534: ND d;
1535: NM m0,m;
1536: MP t;
1537: int n;
1538:
1539: if ( !p )
1540: return 0;
1541: n = NV(p);
1542: m0 = 0;
1543: for ( t = BDY(p); t; t = NEXT(t) ) {
1544: NEXTNM(m0,m);
1545: m->c = ITOS(t->c);
1546: m->td = t->dl->td;
1547: dltondl(n,t->dl,m->dl);
1548: }
1549: NEXT(m) = 0;
1550: MKND(n,m0,d);
1551: d->nv = n;
1552: d->sugar = p->sugar;
1553: return d;
1554: }
1555:
1556: DP ndtodp(ND p)
1557: {
1558: DP d;
1559: MP m0,m;
1560: NM t;
1561: int n;
1562:
1563: if ( !p )
1564: return 0;
1565: n = NV(p);
1566: m0 = 0;
1567: for ( t = BDY(p); t; t = NEXT(t) ) {
1568: NEXTMP(m0,m);
1569: m->c = STOI(t->c);
1570: m->dl = ndltodl(n,t->td,t->dl);
1571: }
1572: NEXT(m) = 0;
1573: MKDP(n,m0,d);
1574: d->sugar = p->sugar;
1575: return d;
1576: }
1577:
1578: void ndl_print(unsigned int *dl)
1579: {
1580: int n;
1581: int i;
1582:
1583: n = nd_nvar;
1584: printf("<<");
1585: if ( is_rlex )
1586: for ( i = 0; i < n; i++ )
1587: printf(i==n-1?"%d":"%d,",
1588: (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
1589: &((1<<nd_bpe)-1));
1590: else
1591: for ( i = 0; i < n; i++ )
1592: printf(i==n-1?"%d":"%d,",
1593: (dl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
1594: &((1<<nd_bpe)-1));
1595: printf(">>");
1596: }
1597:
1598: void nd_print(ND p)
1599: {
1600: NM m;
1601:
1602: if ( !p )
1603: printf("0\n");
1604: else {
1605: for ( m = BDY(p); m; m = NEXT(m) ) {
1606: printf("+%d*",m->c);
1607: ndl_print(m->dl);
1608: }
1609: printf("\n");
1610: }
1611: }
1612:
1613: void ndp_print(ND_pairs d)
1614: {
1615: ND_pairs t;
1616:
1617: for ( t = d; t; t = NEXT(t) ) {
1618: printf("%d,%d ",t->i1,t->i2);
1619: }
1620: printf("\n");
1621: }
1622:
1623: void nd_monic(ND p)
1624: {
1625: if ( !p )
1626: return;
1627: else
1628: nd_mul_c(p,invm(HC(p),nd_mod));
1629: }
1630:
1631: void nd_mul_c(ND p,int mul)
1632: {
1633: NM m;
1634: int c,c1;
1635:
1636: if ( !p )
1637: return;
1638: for ( m = BDY(p); m; m = NEXT(m) ) {
1639: c1 = C(m);
1640: DMAR(c1,mul,0,nd_mod,c);
1641: C(m) = c;
1642: }
1643: }
1644:
1645: void nd_free(ND p)
1646: {
1647: NM t,s;
1648:
1649: if ( !p )
1650: return;
1651: t = BDY(p);
1652: while ( t ) {
1653: s = NEXT(t);
1654: FREENM(t);
1655: t = s;
1656: }
1657: FREEND(p);
1658: }
1659:
1660: void nd_append_red(unsigned int *d,int td,int i)
1661: {
1662: NM m,m0;
1663: int h;
1664:
1665: NEWNM(m);
1666: h = ndl_hash_value(td,d);
1667: m->c = i;
1668: m->td = td;
1669: bcopy(d,m->dl,nd_wpd*sizeof(unsigned int));
1670: NEXT(m) = nd_red[h];
1671: nd_red[h] = m;
1672: }
1673:
1674: unsigned int *nd_compute_bound(ND p)
1675: {
1676: unsigned int *d1,*d2,*t;
1677: NM m;
1678:
1679: if ( !p )
1680: return 0;
1681: d1 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1682: d2 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1683: bcopy(HDL(p),d1,nd_wpd*sizeof(unsigned int));
1684: for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
1685: ndl_lcm(m->dl,d1,d2);
1686: t = d1; d1 = d2; d2 = t;
1687: }
1688: t = (unsigned int *)MALLOC_ATOMIC(nd_wpd*sizeof(unsigned int));
1689: bcopy(d1,t,nd_wpd*sizeof(unsigned int));
1690: return t;
1691: }
1692:
1693: void nd_setup_parameters() {
1694: int i;
1695:
1696: nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
1697: nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
1698: if ( nd_bpe < 32 ) {
1699: nd_mask0 = (1<<nd_bpe)-1;
1700: } else {
1701: nd_mask0 = 0xffffffff;
1702: }
1703: bzero(nd_mask,sizeof(nd_mask));
1704: nd_mask1 = 0;
1705: for ( i = 0; i < nd_epw; i++ ) {
1706: nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
1707: nd_mask1 |= (1<<(nd_bpe-1))<<(i*nd_bpe);
1708: }
1709: }
1710:
1711: ND_pairs nd_reconstruct(ND_pairs d)
1712: {
1713: int i,obpe;
1714: NM prev_nm_free_list;
1715: ND_pairs s0,s,t,prev_ndp_free_list;
1716:
1717: obpe = nd_bpe;
1718: switch ( nd_bpe ) {
1719: case 4: nd_bpe = 6; break;
1720: case 6: nd_bpe = 8; break;
1721: case 8: nd_bpe = 16; break;
1722: case 16: nd_bpe = 32; break;
1723: }
1724: nd_setup_parameters();
1725: prev_nm_free_list = _nm_free_list;
1726: prev_ndp_free_list = _ndp_free_list;
1727: _nm_free_list = 0;
1728: _ndp_free_list = 0;
1729: for ( i = 0; i < nd_psn; i++ ) {
1730: nd_ps[i] = nd_dup(nd_ps[i],obpe);
1731: nd_bound[i] = nd_compute_bound(nd_ps[i]);
1732: }
1733: s0 = 0;
1734: for ( t = d; t; t = NEXT(t) ) {
1735: NEXTND_pairs(s0,s);
1736: s->i1 = t->i1;
1737: s->i2 = t->i2;
1738: s->td = t->td;
1739: s->sugar = t->sugar;
1740: ndl_dup(obpe,t->lcm,s->lcm);
1741: }
1742: if ( s0 ) NEXT(s) = 0;
1743: prev_nm_free_list = 0;
1744: prev_ndp_free_list = 0;
1745: GC_gcollect();
1746: return s0;
1747: }
1748:
1749: void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
1750: {
1751: int n,i,ei,oepw,cepw,cbpe;
1752:
1753: n = nd_nvar;
1754: oepw = (sizeof(unsigned int)*8)/obpe;
1755: cepw = nd_epw;
1756: cbpe = nd_bpe;
1757: if ( is_rlex )
1758: for ( i = 0; i < n; i++ ) {
1759: ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe))
1760: &((1<<obpe)-1);
1761: r[(n-1-i)/cepw] |= (ei<<((cepw-((n-1-i)%cepw)-1)*cbpe));
1762: }
1763: else
1764: for ( i = 0; i < n; i++ ) {
1765: ei = (d[i/oepw]>>((oepw-(i%oepw)-1)*obpe))
1766: &((1<<obpe)-1);
1767: r[i/cepw] |= (ei<<((cepw-(i%cepw)-1)*cbpe));
1768: }
1769: }
1770:
1771: ND nd_dup(ND p,int obpe)
1772: {
1773: NM m,mr,mr0;
1774: int c,n;
1775: ND r;
1776:
1777: if ( !p )
1778: return 0;
1779: else {
1780: n = NV(p); m = BDY(p);
1781: for ( mr0 = 0; m; m = NEXT(m) ) {
1782: NEXTNM(mr0,mr);
1783: C(mr) = C(m);
1784: mr->td = m->td;
1785: ndl_dup(obpe,m->dl,mr->dl);
1786: }
1787: NEXT(mr) = 0;
1788: MKND(NV(p),mr0,r);
1789: r->sugar = p->sugar;
1790: return r;
1791: }
1792: }