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