![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to distm.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 distm.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.5 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.6 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.5 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.15 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.14 2003/07/22 07:12:41 noro Exp $
1.5 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "inline.h"
52:
53: extern int (*cmpdl)();
1.2 noro 54: extern int do_weyl;
55:
1.9 noro 56: void ptomd(VL vl,int mod,VL dvl,P p,DP *pr)
1.1 noro 57: {
58: P t;
59:
60: ptomp(mod,p,&t);
61: mptomd(vl,mod,dvl,t,pr);
62: }
63:
1.9 noro 64: void mptomd(VL vl,int mod,VL dvl,P p,DP *pr)
1.1 noro 65: {
66: int n,i;
67: VL tvl;
68: V v;
69: DL d;
70: MP m;
71: DCP dc;
72: DP r,s,t,u;
73: P x,c;
74:
75: if ( !p )
76: *pr = 0;
77: else {
78: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
79: if ( NUM(p) ) {
80: NEWDL(d,n);
81: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr);
82: } else {
83: for ( i = 0, tvl = dvl, v = VR(p);
84: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
85: if ( !tvl ) {
86: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
87: mptomd(vl,mod,dvl,COEF(dc),&t); pwrmp(vl,mod,x,DEG(dc),&c);
88: mulmdc(vl,mod,t,c,&r); addmd(vl,mod,r,s,&t); s = t;
89: }
90: *pr = s;
91: } else {
92: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
93: mptomd(vl,mod,dvl,COEF(dc),&t);
1.10 noro 94: NEWDL(d,n); d->d[i] = QTOS(DEG(dc));
95: d->td = MUL_WEIGHT(d->d[i],i);
1.1 noro 96: NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u);
1.2 noro 97: comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t;
1.1 noro 98: }
99: *pr = s;
100: }
101: }
102: }
103: }
104:
1.9 noro 105: void mdtop(VL vl,int mod,VL dvl,DP p,P *pr)
1.1 noro 106: {
107: int n,i;
108: DL d;
109: MP m;
110: P r,s,t,u,w;
111: Q q;
112: VL tvl;
113:
114: if ( !p )
115: *pr = 0;
116: else {
117: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
118: for ( i = 0, t = C(m), d = m->dl, tvl = dvl;
119: i < n; tvl = NEXT(tvl), i++ ) {
120: MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
121: mulmp(vl,mod,t,u,&w); t = w;
122: }
123: addmp(vl,mod,s,t,&u); s = u;
124: }
125: mptop(s,pr);
126: }
127: }
128:
1.9 noro 129: void addmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.1 noro 130: {
131: int n;
132: MP m1,m2,mr,mr0;
133: P t;
134: int tmp;
135: MQ q;
136:
137: if ( !p1 )
138: *pr = p2;
139: else if ( !p2 )
140: *pr = p1;
141: else {
142: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
143: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
144: case 0:
145: if ( NUM(C(m1)) && NUM(C(m2)) ) {
146: tmp = (CONT((MQ)C(m1))+CONT((MQ)C(m2))) - mod;
147: if ( tmp < 0 )
148: tmp += mod;
149: if ( tmp ) {
150: STOMQ(tmp,q); t = (P)q;
151: } else
152: t = 0;
153: } else
154: addmp(vl,mod,C(m1),C(m2),&t);
155: if ( t ) {
156: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
157: }
158: m1 = NEXT(m1); m2 = NEXT(m2); break;
159: case 1:
160: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
161: m1 = NEXT(m1); break;
162: case -1:
163: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
164: m2 = NEXT(m2); break;
165: }
166: if ( !mr0 )
167: if ( m1 )
168: mr0 = m1;
169: else if ( m2 )
170: mr0 = m2;
171: else {
172: *pr = 0;
173: return;
174: }
175: else if ( m1 )
176: NEXT(mr) = m1;
177: else if ( m2 )
178: NEXT(mr) = m2;
179: else
180: NEXT(mr) = 0;
181: MKDP(NV(p1),mr0,*pr);
182: if ( *pr )
183: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
184: }
185: }
186:
1.9 noro 187: void submd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.1 noro 188: {
189: DP t;
190:
191: if ( !p2 )
192: *pr = p1;
193: else {
194: chsgnmd(mod,p2,&t); addmd(vl,mod,p1,t,pr);
195: }
196: }
197:
1.9 noro 198: void chsgnmd(int mod,DP p,DP *pr)
1.1 noro 199: {
200: MP m,mr,mr0;
201:
202: if ( !p )
203: *pr = 0;
204: else {
205: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
206: NEXTMP(mr0,mr); chsgnmp(mod,C(m),&C(mr)); mr->dl = m->dl;
207: }
208: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
209: if ( *pr )
210: (*pr)->sugar = p->sugar;
211: }
212: }
213:
1.9 noro 214: void mulmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.1 noro 215: {
1.2 noro 216: if ( !do_weyl )
217: comm_mulmd(vl,mod,p1,p2,pr);
218: else
219: weyl_mulmd(vl,mod,p1,p2,pr);
220: }
221:
1.9 noro 222: void comm_mulmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.2 noro 223: {
224: MP m;
225: DP s,t,u;
226: int i,l,l1;
227: static MP *w;
228: static int wlen;
229:
230: if ( !p1 || !p2 )
231: *pr = 0;
232: else if ( OID(p1) <= O_P )
233: mulmdc(vl,mod,p2,(P)p1,pr);
234: else if ( OID(p2) <= O_P )
235: mulmdc(vl,mod,p1,(P)p2,pr);
236: else {
237: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
238: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
239: if ( l1 < l ) {
240: t = p1; p1 = p2; p2 = t;
241: l = l1;
242: }
243: if ( l > wlen ) {
244: if ( w ) GC_free(w);
245: w = (MP *)MALLOC(l*sizeof(MP));
246: wlen = l;
247: }
248: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
249: w[i] = m;
250: for ( s = 0, i = l-1; i >= 0; i-- ) {
251: mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
252: }
253: bzero(w,l*sizeof(MP));
254: *pr = s;
255: }
256: }
257:
1.9 noro 258: void weyl_mulmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.2 noro 259: {
1.1 noro 260: MP m;
261: DP s,t,u;
1.9 noro 262: int i,l;
1.2 noro 263: static MP *w;
264: static int wlen;
1.1 noro 265:
266: if ( !p1 || !p2 )
267: *pr = 0;
268: else if ( OID(p1) <= O_P )
269: mulmdc(vl,mod,p2,(P)p1,pr);
270: else if ( OID(p2) <= O_P )
271: mulmdc(vl,mod,p1,(P)p2,pr);
272: else {
1.2 noro 273: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
274: if ( l > wlen ) {
275: if ( w ) GC_free(w);
276: w = (MP *)MALLOC(l*sizeof(MP));
277: wlen = l;
278: }
279: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
280: w[i] = m;
281: for ( s = 0, i = l-1; i >= 0; i-- ) {
282: weyl_mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
1.1 noro 283: }
1.2 noro 284: bzero(w,l*sizeof(MP));
1.1 noro 285: *pr = s;
286: }
287: }
288:
1.9 noro 289: void mulmdm(VL vl,int mod,DP p,MP m0,DP *pr)
1.1 noro 290: {
291: MP m,mr,mr0;
292: P c;
293: DL d;
294: int n,t;
295: MQ q;
296:
297: if ( !p )
298: *pr = 0;
299: else {
300: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
301: m; m = NEXT(m) ) {
302: NEXTMP(mr0,mr);
303: if ( NUM(C(m)) && NUM(c) ) {
304: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
305: STOMQ(t,q); C(mr) = (P)q;
306: } else
307: mulmp(vl,mod,C(m),c,&C(mr));
308: adddl(n,m->dl,d,&mr->dl);
309: }
310: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
311: if ( *pr )
312: (*pr)->sugar = p->sugar + m0->dl->td;
313: }
314: }
315:
1.9 noro 316: void weyl_mulmdm(VL vl,int mod,DP p,MP m0,DP *pr)
1.2 noro 317: {
318: DP r,t,t1;
319: MP m;
320: int n,l,i;
321: static MP *w;
322: static int wlen;
323:
324: if ( !p )
325: *pr = 0;
326: else {
327: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
328: if ( l > wlen ) {
329: if ( w ) GC_free(w);
330: w = (MP *)MALLOC(l*sizeof(MP));
331: wlen = l;
332: }
333: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
334: w[i] = m;
335: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
336: weyl_mulmmm(vl,mod,w[i],m0,n,&t);
337: addmd(vl,mod,r,t,&t1); r = t1;
338: }
339: bzero(w,l*sizeof(MP));
340: if ( r )
341: r->sugar = p->sugar + m0->dl->td;
342: *pr = r;
343: }
344: }
345:
346: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
347:
1.9 noro 348: void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *pr)
1.2 noro 349: {
1.9 noro 350: MP mr,mr0;
1.2 noro 351: MQ mq;
352: DP r,t,t1;
1.9 noro 353: P c,c0,c1;
1.2 noro 354: DL d,d0,d1;
1.3 noro 355: int i,j,a,b,k,l,n2,s,min;
1.2 noro 356: static int *tab;
357: static int tablen;
358:
359: if ( !m0 || !m1 )
360: *pr = 0;
361: else {
362: c0 = C(m0); c1 = C(m1);
363: mulmp(vl,mod,c0,c1,&c);
364: d0 = m0->dl; d1 = m1->dl;
365: n2 = n>>1;
366: if ( n & 1 ) {
367: /* homogenized computation; dx-xd=h^2 */
368: /* offset of h-degree */
1.3 noro 369: NEWDL(d,n);
370: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1.4 noro 371: NEWMP(mr); mr->c = (P)ONEM; mr->dl = d; NEXT(mr) = 0;
1.3 noro 372: MKDP(n,mr,r); r->sugar = d->td;
373: } else
374: r = (DP)ONEM;
375: for ( i = 0; i < n2; i++ ) {
376: a = d0->d[i]; b = d1->d[n2+i];
377: k = d0->d[n2+i]; l = d1->d[i];
1.10 noro 378:
1.3 noro 379: /* degree of xi^a*(Di^k*xi^l)*Di^b */
1.10 noro 380: a += l;
381: b += k;
382: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
383:
1.3 noro 384: /* compute xi^a*(Di^k*xi^l)*Di^b */
385: min = MIN(k,l);
386:
387: if ( min+1 > tablen ) {
388: if ( tab ) GC_free(tab);
389: tab = (int *)MALLOC((min+1)*sizeof(int));
390: tablen = min+1;
391: }
392: mkwcm(k,l,mod,tab);
393: if ( n & 1 )
1.2 noro 394: for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3 noro 395: NEXTMP(mr0,mr); NEWDL(d,n);
1.10 noro 396: d->d[i] = a-j; d->d[n2+i] = b-j;
1.3 noro 397: d->td = s;
1.10 noro 398: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.3 noro 399: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 400: }
1.3 noro 401: else
1.2 noro 402: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3 noro 403: NEXTMP(mr0,mr); NEWDL(d,n);
1.10 noro 404: d->d[i] = a-j; d->d[n2+i] = b-j;
405: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.2 noro 406: s = MAX(s,d->td); /* XXX */
1.3 noro 407: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 408: }
1.3 noro 409: bzero(tab,(min+1)*sizeof(int));
410: if ( mr0 )
411: NEXT(mr) = 0;
412: MKDP(n,mr0,t);
413: if ( t )
414: t->sugar = s;
415: comm_mulmd(vl,mod,r,t,&t1); r = t1;
416: }
1.2 noro 417: mulmdc(vl,mod,r,c,pr);
418: }
419: }
420:
1.9 noro 421: void mulmdc(VL vl,int mod,DP p,P c,DP *pr)
1.1 noro 422: {
423: MP m,mr,mr0;
424: int t;
425: MQ q;
426:
427: if ( !p || !c )
428: *pr = 0;
429: else {
430: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
431: NEXTMP(mr0,mr);
432: if ( NUM(C(m)) && NUM(c) ) {
433: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
434: STOMQ(t,q); C(mr) = (P)q;
435: } else
436: mulmp(vl,mod,C(m),c,&C(mr));
437: mr->dl = m->dl;
438: }
439: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
440: if ( *pr )
441: (*pr)->sugar = p->sugar;
442: }
443: }
444:
1.9 noro 445: void divsmdc(VL vl,int mod,DP p,P c,DP *pr)
1.1 noro 446: {
447: MP m,mr,mr0;
448:
449: if ( !c )
450: error("disvsdc : division by 0");
451: else if ( !p )
452: *pr = 0;
453: else {
454: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
455: NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl;
456: }
457: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
458: if ( *pr )
459: (*pr)->sugar = p->sugar;
460: }
461: }
462:
1.9 noro 463: void _dtop_mod(VL vl,VL dvl,DP p,P *pr)
1.1 noro 464: {
465: int n,i;
466: DL d;
467: MP m;
468: P r,s,t,u,w;
469: Q q;
470: VL tvl;
471:
472: if ( !p )
473: *pr = 0;
474: else {
475: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
476: i = ITOS(m->c); STOQ(i,q); t = (P)q;
477: for ( i = 0, d = m->dl, tvl = dvl;
478: i < n; tvl = NEXT(tvl), i++ ) {
479: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
480: mulp(vl,t,u,&w); t = w;
481: }
482: addp(vl,s,t,&u); s = u;
483: }
484: *pr = s;
485: }
486: }
487:
1.9 noro 488: void _dp_mod(DP p,int mod,NODE subst,DP *rp)
1.1 noro 489: {
490: MP m,mr,mr0;
491: P t,s;
492: NODE tn;
493:
494: if ( !p )
495: *rp = 0;
496: else {
497: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
498: for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) {
499: substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t);
500: s = t;
501: }
502: ptomp(mod,s,&t);
503: if ( t ) {
504: NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl;
505: }
506: }
507: if ( mr0 ) {
508: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
509: } else
510: *rp = 0;
511: }
512: }
513:
1.9 noro 514: void _dp_monic(DP p,int mod,DP *rp)
1.2 noro 515: {
516: MP m,mr,mr0;
517: int c,c1;
518:
519: if ( !p )
520: *rp = 0;
521: else {
522: c = invm(ITOS(BDY(p)->c),mod);
523: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
524: c1 = dmar(ITOS(m->c),c,0,mod);
525: NEXTMP(mr0,mr); mr->c = STOI(c1); mr->dl = m->dl;
526: }
527: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
528: }
1.1 noro 529: }
530:
1.9 noro 531: void _printdp(DP d)
1.1 noro 532: {
533: int n,i;
534: MP m;
535: DL dl;
536:
537: if ( !d ) {
538: printf("0"); return;
539: }
540: for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) {
541: printf("%d*<<",ITOS(m->c));
542: for ( i = 0, dl = m->dl; i < n-1; i++ )
543: printf("%d,",dl->d[i]);
544: printf("%d",dl->d[i]);
545: printf(">>");
546: if ( NEXT(m) )
547: printf("+");
548: }
1.8 noro 549: }
550:
551: /* merge p1 and p2 into pr */
552:
1.9 noro 553: void addmd_destructive(int mod,DP p1,DP p2,DP *pr)
1.8 noro 554: {
555: int n;
556: MP m1,m2,mr,mr0,s;
557: int t;
558:
559: if ( !p1 )
560: *pr = p2;
561: else if ( !p2 )
562: *pr = p1;
563: else {
564: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
565: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
566: case 0:
567: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
568: if ( t < 0 )
569: t += mod;
570: s = m1; m1 = NEXT(m1);
571: if ( t ) {
572: NEXTMP2(mr0,mr,s); C(mr) = STOI(t);
573: }
574: s = m2; m2 = NEXT(m2);
575: break;
576: case 1:
577: s = m1; m1 = NEXT(m1); NEXTMP2(mr0,mr,s);
578: break;
579: case -1:
580: s = m2; m2 = NEXT(m2); NEXTMP2(mr0,mr,s);
581: break;
582: }
583: if ( !mr0 )
584: if ( m1 )
585: mr0 = m1;
586: else if ( m2 )
587: mr0 = m2;
588: else {
589: *pr = 0;
590: return;
591: }
592: else if ( m1 )
593: NEXT(mr) = m1;
594: else if ( m2 )
595: NEXT(mr) = m2;
596: else
597: NEXT(mr) = 0;
598: MKDP(NV(p1),mr0,*pr);
599: if ( *pr )
600: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
601: }
602: }
603:
1.9 noro 604: void mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.8 noro 605: {
606: if ( !do_weyl )
607: comm_mulmd_dup(mod,p1,p2,pr);
608: else
609: weyl_mulmd_dup(mod,p1,p2,pr);
610: }
611:
1.9 noro 612: void comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.8 noro 613: {
614: MP m;
615: DP s,t,u;
616: int i,l,l1;
617: static MP *w;
618: static int wlen;
619:
620: if ( !p1 || !p2 )
621: *pr = 0;
622: else {
623: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
624: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
625: if ( l1 < l ) {
626: t = p1; p1 = p2; p2 = t;
627: l = l1;
628: }
629: if ( l > wlen ) {
630: if ( w ) GC_free(w);
631: w = (MP *)MALLOC(l*sizeof(MP));
632: wlen = l;
633: }
634: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
635: w[i] = m;
636: for ( s = 0, i = l-1; i >= 0; i-- ) {
637: mulmdm_dup(mod,p1,w[i],&t); addmd_destructive(mod,s,t,&u); s = u;
638: }
639: bzero(w,l*sizeof(MP));
640: *pr = s;
641: }
642: }
643:
1.9 noro 644: void weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.8 noro 645: {
646: MP m;
647: DP s,t,u;
1.9 noro 648: int i,l;
1.8 noro 649: static MP *w;
650: static int wlen;
651:
652: if ( !p1 || !p2 )
653: *pr = 0;
654: else {
655: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
656: if ( l > wlen ) {
657: if ( w ) GC_free(w);
658: w = (MP *)MALLOC(l*sizeof(MP));
659: wlen = l;
660: }
661: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
662: w[i] = m;
663: for ( s = 0, i = l-1; i >= 0; i-- ) {
664: weyl_mulmdm_dup(mod,w[i],p2,&t); addmd_destructive(mod,s,t,&u); s = u;
665: }
666: bzero(w,l*sizeof(MP));
667: *pr = s;
668: }
669: }
670:
1.9 noro 671: void mulmdm_dup(int mod,DP p,MP m0,DP *pr)
1.8 noro 672: {
673: MP m,mr,mr0;
674: DL d,dt,dm;
1.9 noro 675: int c,n,i;
1.8 noro 676: int *pt,*p1,*p2;
677:
678: if ( !p )
679: *pr = 0;
680: else {
681: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
682: m; m = NEXT(m) ) {
683: NEXTMP(mr0,mr);
684: C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
685: NEWDL_NOINIT(dt,n); mr->dl = dt;
686: dm = m->dl;
687: dt->td = d->td + dm->td;
688: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
689: *pt++ = *p1++ + *p2++;
690: }
691: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
692: if ( *pr )
693: (*pr)->sugar = p->sugar + m0->dl->td;
694: }
695: }
696:
1.9 noro 697: void weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr)
1.8 noro 698: {
699: DP r,t,t1;
700: MP m;
701: DL d0;
702: int n,n2,l,i,j,tlen;
703: static MP *w,*psum;
704: static struct cdlm *tab;
705: static int wlen;
706: static int rtlen;
707:
708: if ( !p )
709: *pr = 0;
710: else {
711: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
712: if ( l > wlen ) {
713: if ( w ) GC_free(w);
714: w = (MP *)MALLOC(l*sizeof(MP));
715: wlen = l;
716: }
717: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
718: w[i] = m;
719: n = NV(p); n2 = n>>1;
720: d0 = m0->dl;
721:
722: for ( i = 0, tlen = 1; i < n2; i++ )
723: tlen *= d0->d[n2+i]+1;
724: if ( tlen > rtlen ) {
725: if ( tab ) GC_free(tab);
726: if ( psum ) GC_free(psum);
727: rtlen = tlen;
728: tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
729: psum = (MP *)MALLOC(rtlen*sizeof(MP));
730: }
731: bzero(psum,tlen*sizeof(MP));
732: for ( i = l-1; i >= 0; i-- ) {
733: bzero(tab,tlen*sizeof(struct cdlm));
734: weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
735: for ( j = 0; j < tlen; j++ ) {
736: if ( tab[j].c ) {
737: NEWMP(m); m->dl = tab[j].d;
738: C(m) = STOI(tab[j].c); NEXT(m) = psum[j];
739: psum[j] = m;
740: }
741: }
742: }
743: for ( j = tlen-1, r = 0; j >= 0; j-- )
744: if ( psum[j] ) {
745: MKDP(n,psum[j],t); addmd_destructive(mod,r,t,&t1); r = t1;
746: }
747: if ( r )
748: r->sugar = p->sugar + m0->dl->td;
749: *pr = r;
750: }
751: }
752:
753: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
754:
1.9 noro 755: void weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct cdlm *rtab,int rtablen)
1.8 noro 756: {
1.9 noro 757: int c,c0,c1;
1.8 noro 758: DL d,d0,d1,dt;
1.9 noro 759: int i,j,a,b,k,l,n2,s,min,curlen;
1.8 noro 760: struct cdlm *p;
761: static int *ctab;
762: static struct cdlm *tab;
763: static int tablen;
764: static struct cdlm *tmptab;
765: static int tmptablen;
766:
767: if ( !m0 || !m1 ) {
768: rtab[0].c = 0;
769: rtab[0].d = 0;
770: return;
771: }
772: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
773: c = dmar(c0,c1,0,mod);
774: d0 = m0->dl; d1 = m1->dl;
775: n2 = n>>1;
776: curlen = 1;
777:
778: NEWDL(d,n);
779: if ( n & 1 )
780: /* offset of h-degree */
781: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
782: else
783: d->td = 0;
784: rtab[0].c = c;
785: rtab[0].d = d;
786:
787: if ( rtablen > tmptablen ) {
788: if ( tmptab ) GC_free(tmptab);
789: tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
790: tmptablen = rtablen;
791: }
792:
793: for ( i = 0; i < n2; i++ ) {
794: a = d0->d[i]; b = d1->d[n2+i];
795: k = d0->d[n2+i]; l = d1->d[i];
1.10 noro 796:
797: /* degree of xi^a*(Di^k*xi^l)*Di^b */
798: a += l;
799: b += k;
800: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
801:
1.8 noro 802: if ( !k || !l ) {
803: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
804: if ( p->c ) {
805: dt = p->d;
806: dt->d[i] = a;
807: dt->d[n2+i] = b;
808: dt->td += s;
809: }
810: }
811: curlen *= k+1;
812: continue;
813: }
814: if ( k+1 > tablen ) {
815: if ( tab ) GC_free(tab);
816: if ( ctab ) GC_free(ctab);
817: tablen = k+1;
818: tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
819: ctab = (int *)MALLOC(tablen*sizeof(int));
820: }
821: /* compute xi^a*(Di^k*xi^l)*Di^b */
822: min = MIN(k,l);
823: mkwcm(k,l,mod,ctab);
824: bzero(tab,(k+1)*sizeof(struct cdlm));
825: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
826: if ( n & 1 )
827: for ( j = 0; j <= min; j++ ) {
828: NEWDL(d,n);
1.10 noro 829: d->d[i] = a-j; d->d[n2+i] = b-j;
1.8 noro 830: d->td = s;
1.10 noro 831: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.8 noro 832: tab[j].d = d;
833: tab[j].c = ctab[j];
834: }
835: else
836: for ( j = 0; j <= min; j++ ) {
837: NEWDL(d,n);
1.10 noro 838: d->d[i] = a-j; d->d[n2+i] = b-j;
839: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.8 noro 840: tab[j].d = d;
841: tab[j].c = ctab[j];
842: }
843: comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
844: curlen *= k+1;
845: }
846: }
847:
1.9 noro 848: void comm_mulmd_tab_destructive(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1)
1.8 noro 849: {
850: int i,j;
851: struct cdlm *p;
852: int c;
853: DL d;
854:
855: for ( j = 1, p = t+n; j < n1; j++ ) {
856: c = t1[j].c;
857: d = t1[j].d;
858: if ( !c )
859: break;
860: for ( i = 0; i < n; i++, p++ ) {
861: if ( t[i].c ) {
862: p->c = dmar(t[i].c,c,0,mod);
863: adddl_dup(nv,t[i].d,d,&p->d);
864: }
865: }
866: }
867: c = t1[0].c;
868: d = t1[0].d;
869: for ( i = 0, p = t; i < n; i++, p++ )
870: if ( t[i].c ) {
871: p->c = dmar(t[i].c,c,0,mod);
872: /* t[i].d += d */
873: adddl_destructive(nv,t[i].d,d);
874: }
875: }
876:
1.9 noro 877: void adddl_dup(int n,DL d1,DL d2,DL *dr)
1.8 noro 878: {
879: DL dt;
880: int i;
881:
882: NEWDL(dt,n);
883: *dr = dt;
884: dt->td = d1->td + d2->td;
885: for ( i = 0; i < n; i++ )
886: dt->d[i] = d1->d[i]+d2->d[i];
1.11 noro 887: }
888:
889: /* ------------------------------------------------------ */
890:
891: #if defined(__GNUC__)
892: #define INLINE inline
893: #elif defined(VISUAL)
894: #define INLINE __inline
895: #else
896: #define INLINE
897: #endif
898:
1.14 noro 899: #define REDTAB_LEN 32003
900:
1.11 noro 901: typedef struct oPGeoBucket {
902: int m;
903: struct oND *body[32];
904: } *PGeoBucket;
905:
906: typedef struct oND {
907: struct oNM *body;
908: int nv;
909: int sugar;
910: } *ND;
911:
912: typedef struct oNM {
913: struct oNM *next;
914: int td;
915: int c;
916: unsigned int dl[1];
917: } *NM;
918:
919: typedef struct oND_pairs {
920: struct oND_pairs *next;
921: int i1,i2;
922: int td,sugar;
923: unsigned int lcm[1];
924: } *ND_pairs;
925:
1.15 ! noro 926: static ND *nd_ps;
! 927: static unsigned int **nd_bound;
1.11 noro 928: int nd_mod,nd_nvar;
929: int is_rlex;
930: int nd_epw,nd_bpe,nd_wpd;
1.12 noro 931: unsigned int nd_mask[32];
1.15 ! noro 932: unsigned int nd_mask0,nd_mask1;
1.11 noro 933: NM _nm_free_list;
934: ND _nd_free_list;
935: ND_pairs _ndp_free_list;
1.14 noro 936: NM *nd_red;
937: int nd_red_len;
1.11 noro 938:
939: extern int Top,Reverse;
940: int nd_psn,nd_pslen;
1.14 noro 941: int nd_found,nd_create,nd_notfirst;
1.11 noro 942:
943: void GC_gcollect();
944: NODE append_one(NODE,int);
945:
946: #define HTD(d) ((d)->body->td)
947: #define HDL(d) ((d)->body->dl)
948: #define HC(d) ((d)->body->c)
949:
950: #define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list)
951: #define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list)
952: #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)
953:
954: #define NEXTNM(r,c) \
955: if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);}
956: #define NEXTNM2(r,c,s) \
957: if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
958: #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)
1.15 ! noro 959: #define FREENDP(m) NEXT(m)=_ndp_free_list; _ndp_free_list=(m)
1.11 noro 960: #define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m)
961:
962: #define NEXTND_pairs(r,c) \
963: if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
964:
965: ND_pairs crit_B( ND_pairs d, int s );
966: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp);
967: NODE nd_setup(NODE f);
968: int nd_newps(ND a);
969: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest );
970: NODE update_base(NODE nd,int ndp);
971: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest );
972: int crit_2( int dp1, int dp2 );
973: ND_pairs crit_F( ND_pairs d1 );
974: ND_pairs crit_M( ND_pairs d1 );
975: ND_pairs nd_newpairs( NODE g, int t );
976: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t);
977: NODE nd_gb(NODE f);
978: void nd_free_private_storage();
979: void _NM_alloc();
980: void _ND_alloc();
981: int ndl_td(unsigned int *d);
982: ND nd_add(ND p1,ND p2);
983: ND nd_mul_nm(ND p,NM m0);
1.13 noro 984: ND nd_mul_term(ND p,int td,unsigned int *d);
1.15 ! noro 985: int nd_sp(ND_pairs p,ND *nf);
! 986: int nd_find_reducer(ND g,ND *red);
! 987: int nd_nf(ND g,int full,ND *nf);
1.14 noro 988: ND nd_reduce(ND p1,ND p2);
989: ND nd_reduce_special(ND p1,ND p2);
990: void nd_free(ND p);
1.11 noro 991: void ndl_print(unsigned int *dl);
992: void nd_print(ND p);
993: void ndp_print(ND_pairs d);
994: int nd_length(ND p);
995: void nd_monic(ND p);
1.13 noro 996: void nd_mul_c(ND p,int mul);
1.14 noro 997: void nd_free_redlist();
998: void nd_append_red(unsigned int *d,int td,int i);
1.15 ! noro 999: unsigned int *nd_compute_bound(ND p);
! 1000: ND_pairs nd_reconstruct(ND_pairs);
! 1001: void nd_setup_parameters();
! 1002: ND nd_dup(ND p,int obpe);
! 1003: void ndl_dup(int obpe,unsigned int *d,unsigned int *r);
1.11 noro 1004:
1005: void nd_free_private_storage()
1006: {
1007: _nd_free_list = 0;
1008: _nm_free_list = 0;
1.14 noro 1009: nd_red = 0;
1.11 noro 1010: GC_gcollect();
1011: }
1012:
1013: void _NM_alloc()
1014: {
1015: NM p;
1016: int i;
1017:
1.14 noro 1018: for ( i = 0; i < 16; i++ ) {
1.11 noro 1019: p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
1020: p->next = _nm_free_list; _nm_free_list = p;
1021: }
1022: }
1023:
1024: void _ND_alloc()
1025: {
1026: ND p;
1027: int i;
1028:
1.13 noro 1029: for ( i = 0; i < 1024; i++ ) {
1.11 noro 1030: p = (ND)GC_malloc(sizeof(struct oND));
1031: p->body = (NM)_nd_free_list; _nd_free_list = p;
1032: }
1033: }
1034:
1035: void _NDP_alloc()
1036: {
1037: ND_pairs p;
1038: int i;
1039:
1.12 noro 1040: for ( i = 0; i < 10240; i++ ) {
1.11 noro 1041: p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs)
1042: +(nd_wpd-1)*sizeof(unsigned int));
1043: p->next = _ndp_free_list; _ndp_free_list = p;
1044: }
1045: }
1046:
1047: INLINE nd_length(ND p)
1048: {
1049: NM m;
1050: int i;
1051:
1052: if ( !p )
1053: return 0;
1054: else {
1055: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
1056: return i;
1057: }
1058: }
1059:
1060: INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2)
1061: {
1062: unsigned int u1,u2;
1.12 noro 1063: int i,j;
1.11 noro 1064:
1065: switch ( nd_bpe ) {
1066: case 4:
1067: for ( i = 0; i < nd_wpd; i++ ) {
1068: u1 = d1[i]; u2 = d2[i];
1069: if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0;
1070: if ( (u1&0xf000000) < (u2&0xf000000) ) return 0;
1071: if ( (u1&0xf00000) < (u2&0xf00000) ) return 0;
1072: if ( (u1&0xf0000) < (u2&0xf0000) ) return 0;
1073: if ( (u1&0xf000) < (u2&0xf000) ) return 0;
1074: if ( (u1&0xf00) < (u2&0xf00) ) return 0;
1075: if ( (u1&0xf0) < (u2&0xf0) ) return 0;
1076: if ( (u1&0xf) < (u2&0xf) ) return 0;
1077: }
1078: return 1;
1079: break;
1.13 noro 1080: case 6:
1081: for ( i = 0; i < nd_wpd; i++ ) {
1082: u1 = d1[i]; u2 = d2[i];
1083: if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0;
1084: if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0;
1085: if ( (u1&0x3f000) < (u2&0x3f000) ) return 0;
1086: if ( (u1&0xfc0) < (u2&0xfc0) ) return 0;
1087: if ( (u1&0x3f) < (u2&0x3f) ) return 0;
1088: }
1089: return 1;
1090: break;
1.11 noro 1091: case 8:
1092: for ( i = 0; i < nd_wpd; i++ ) {
1093: u1 = d1[i]; u2 = d2[i];
1094: if ( (u1&0xff000000) < (u2&0xff000000) ) return 0;
1095: if ( (u1&0xff0000) < (u2&0xff0000) ) return 0;
1096: if ( (u1&0xff00) < (u2&0xff00) ) return 0;
1097: if ( (u1&0xff) < (u2&0xff) ) return 0;
1098: }
1099: return 1;
1100: break;
1101: case 16:
1102: for ( i = 0; i < nd_wpd; i++ ) {
1103: u1 = d1[i]; u2 = d2[i];
1104: if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0;
1105: if ( (u1&0xffff) < (u2&0xffff) ) return 0;
1106: }
1107: return 1;
1108: break;
1109: case 32:
1110: for ( i = 0; i < nd_wpd; i++ )
1111: if ( d1[i] < d2[i] ) return 0;
1112: return 1;
1113: break;
1114: default:
1.12 noro 1115: for ( i = 0; i < nd_wpd; i++ ) {
1116: u1 = d1[i]; u2 = d2[i];
1117: for ( j = 0; j < nd_epw; j++ )
1118: if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0;
1119: }
1120: return 1;
1.11 noro 1121: }
1122: }
1123:
1.12 noro 1124: void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
1.11 noro 1125: {
1126: unsigned int t1,t2,u,u1,u2;
1.12 noro 1127: int i,j;
1.11 noro 1128:
1129: switch ( nd_bpe ) {
1130: case 4:
1131: for ( i = 0; i < nd_wpd; i++ ) {
1132: u1 = d1[i]; u2 = d2[i];
1133: t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2;
1134: t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2;
1135: t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2;
1136: t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2;
1137: t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2;
1138: t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2;
1139: t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2;
1140: t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2;
1141: d[i] = u;
1142: }
1143: break;
1.13 noro 1144: case 6:
1145: for ( i = 0; i < nd_wpd; i++ ) {
1146: u1 = d1[i]; u2 = d2[i];
1147: t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2;
1148: t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2;
1149: t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2;
1150: t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2;
1151: t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2;
1152: d[i] = u;
1153: }
1154: break;
1.11 noro 1155: case 8:
1156: for ( i = 0; i < nd_wpd; i++ ) {
1157: u1 = d1[i]; u2 = d2[i];
1158: t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2;
1159: t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2;
1160: t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2;
1161: t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2;
1162: d[i] = u;
1163: }
1164: break;
1165: case 16:
1166: for ( i = 0; i < nd_wpd; i++ ) {
1167: u1 = d1[i]; u2 = d2[i];
1168: t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2;
1169: t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2;
1170: d[i] = u;
1171: }
1172: break;
1173: case 32:
1174: for ( i = 0; i < nd_wpd; i++ ) {
1175: u1 = d1[i]; u2 = d2[i];
1176: d[i] = u1>u2?u1:u2;
1177: }
1178: break;
1179: default:
1.12 noro 1180: for ( i = 0; i < nd_wpd; i++ ) {
1181: u1 = d1[i]; u2 = d2[i];
1182: for ( j = 0, u = 0; j < nd_epw; j++ ) {
1183: t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2;
1184: }
1185: d[i] = u;
1186: }
1187: break;
1.11 noro 1188: }
1189: }
1190:
1.12 noro 1191: int ndl_td(unsigned int *d)
1.11 noro 1192: {
1193: unsigned int t,u;
1.12 noro 1194: int i,j;
1.11 noro 1195:
1.12 noro 1196: for ( t = 0, i = 0; i < nd_wpd; i++ ) {
1197: u = d[i];
1198: for ( j = 0; j < nd_epw; j++, u>>=nd_bpe )
1199: t += (u&nd_mask0);
1.11 noro 1200: }
1201: return t;
1202: }
1203:
1204: INLINE int ndl_compare(unsigned int *d1,unsigned int *d2)
1205: {
1206: int i;
1207:
1.12 noro 1208: for ( i = 0; i < nd_wpd; i++, d1++, d2++ )
1209: if ( *d1 > *d2 )
1.11 noro 1210: return is_rlex ? -1 : 1;
1.12 noro 1211: else if ( *d1 < *d2 )
1.11 noro 1212: return is_rlex ? 1 : -1;
1213: return 0;
1214: }
1215:
1216: INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)
1217: {
1218: int i;
1219:
1220: for ( i = 0; i < nd_wpd; i++ )
1221: if ( d1[i] != d2[i] )
1222: return 0;
1223: return 1;
1224: }
1225:
1226: INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d)
1227: {
1228: int i;
1229:
1.15 ! noro 1230: for ( i = 0; i < nd_wpd; i++ ) {
1.11 noro 1231: d[i] = d1[i]+d2[i];
1.15 ! noro 1232: }
1.11 noro 1233: }
1234:
1.12 noro 1235: void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
1.11 noro 1236: {
1237: int i;
1238:
1239: for ( i = 0; i < nd_wpd; i++ )
1240: d[i] = d1[i]-d2[i];
1241: }
1242:
1.12 noro 1243: int ndl_disjoint(unsigned int *d1,unsigned int *d2)
1.11 noro 1244: {
1245: unsigned int t1,t2,u,u1,u2;
1.12 noro 1246: int i,j;
1.11 noro 1247:
1248: switch ( nd_bpe ) {
1249: case 4:
1250: for ( i = 0; i < nd_wpd; i++ ) {
1251: u1 = d1[i]; u2 = d2[i];
1252: t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0;
1253: t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0;
1254: t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0;
1255: t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0;
1256: t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0;
1257: t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0;
1258: t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0;
1259: t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0;
1260: }
1261: return 1;
1262: break;
1.13 noro 1263: case 6:
1264: for ( i = 0; i < nd_wpd; i++ ) {
1265: u1 = d1[i]; u2 = d2[i];
1266: t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0;
1267: t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0;
1268: t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0;
1269: t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0;
1270: t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0;
1271: }
1272: return 1;
1273: break;
1.11 noro 1274: case 8:
1275: for ( i = 0; i < nd_wpd; i++ ) {
1276: u1 = d1[i]; u2 = d2[i];
1277: t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0;
1278: t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0;
1279: t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0;
1280: t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0;
1281: }
1282: return 1;
1283: break;
1284: case 16:
1285: for ( i = 0; i < nd_wpd; i++ ) {
1286: u1 = d1[i]; u2 = d2[i];
1287: t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0;
1288: t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0;
1289: }
1290: return 1;
1291: break;
1292: case 32:
1293: for ( i = 0; i < nd_wpd; i++ )
1294: if ( d1[i] && d2[i] ) return 0;
1295: return 1;
1296: break;
1297: default:
1.12 noro 1298: for ( i = 0; i < nd_wpd; i++ ) {
1299: u1 = d1[i]; u2 = d2[i];
1300: for ( j = 0; j < nd_epw; j++ ) {
1301: if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0;
1302: u1 >>= nd_bpe; u2 >>= nd_bpe;
1303: }
1.11 noro 1304: }
1.12 noro 1305: return 1;
1306: break;
1.11 noro 1307: }
1308: }
1309:
1310: ND nd_reduce(ND p1,ND p2)
1311: {
1.13 noro 1312: int c,c1,c2,t,td,td2,mul;
1.11 noro 1313: NM m2,prev,head,cur,new;
1314: unsigned int *d;
1315:
1316: if ( !p1 )
1317: return 0;
1318: else {
1.13 noro 1319: c2 = invm(HC(p2),nd_mod);
1320: c1 = nd_mod-HC(p1);
1321: DMAR(c1,c2,0,nd_mod,mul);
1.11 noro 1322: td = HTD(p1)-HTD(p2);
1323: d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
1324: ndl_sub(HDL(p1),HDL(p2),d);
1325: prev = 0; head = cur = BDY(p1);
1326: NEWNM(new);
1327: for ( m2 = BDY(p2); m2; ) {
1328: td2 = new->td = m2->td+td;
1329: ndl_add(m2->dl,d,new->dl);
1330: if ( !cur ) {
1.13 noro 1331: c1 = C(m2);
1332: DMAR(c1,mul,0,nd_mod,c2);
1333: C(new) = c2;
1.11 noro 1334: if ( !prev ) {
1335: prev = new;
1336: NEXT(prev) = 0;
1337: head = prev;
1338: } else {
1339: NEXT(prev) = new;
1340: NEXT(new) = 0;
1341: prev = new;
1342: }
1343: m2 = NEXT(m2);
1344: NEWNM(new);
1345: continue;
1346: }
1347: if ( cur->td > td2 )
1348: c = 1;
1349: else if ( cur->td < td2 )
1350: c = -1;
1351: else
1352: c = ndl_compare(cur->dl,new->dl);
1353: switch ( c ) {
1354: case 0:
1.13 noro 1355: c2 = C(m2);
1356: c1 = C(cur);
1357: DMAR(c2,mul,c1,nd_mod,t);
1.11 noro 1358: if ( t )
1359: C(cur) = t;
1360: else if ( !prev ) {
1361: head = NEXT(cur);
1362: FREENM(cur);
1363: cur = head;
1364: } else {
1365: NEXT(prev) = NEXT(cur);
1366: FREENM(cur);
1367: cur = NEXT(prev);
1368: }
1369: m2 = NEXT(m2);
1370: break;
1371: case 1:
1372: prev = cur;
1373: cur = NEXT(cur);
1374: break;
1375: case -1:
1376: if ( !prev ) {
1377: /* cur = head */
1378: prev = new;
1.13 noro 1379: c2 = C(m2);
1380: DMAR(c2,mul,0,nd_mod,c1);
1381: C(prev) = c1;
1.11 noro 1382: NEXT(prev) = head;
1383: head = prev;
1384: } else {
1.13 noro 1385: c2 = C(m2);
1386: DMAR(c2,mul,0,nd_mod,c1);
1387: C(new) = c1;
1.11 noro 1388: NEXT(prev) = new;
1389: NEXT(new) = cur;
1390: prev = new;
1391: }
1392: NEWNM(new);
1393: m2 = NEXT(m2);
1394: break;
1395: }
1396: }
1397: FREENM(new);
1398: if ( head ) {
1399: BDY(p1) = head;
1400: p1->sugar = MAX(p1->sugar,p2->sugar+td);
1401: return p1;
1402: } else {
1403: FREEND(p1);
1404: return 0;
1405: }
1406:
1407: }
1408: }
1409:
1.14 noro 1410: /* HDL(p1) = HDL(p2) */
1411:
1412: ND nd_reduce_special(ND p1,ND p2)
1413: {
1414: int c,c1,c2,t,td,td2,mul;
1415: NM m2,prev,head,cur,new;
1416:
1417: if ( !p1 )
1418: return 0;
1419: else {
1420: c2 = invm(HC(p2),nd_mod);
1421: c1 = nd_mod-HC(p1);
1422: DMAR(c1,c2,0,nd_mod,mul);
1423: prev = 0; head = cur = BDY(p1);
1424: NEWNM(new);
1425: for ( m2 = BDY(p2); m2; ) {
1426: td2 = new->td = m2->td;
1427: if ( !cur ) {
1428: c1 = C(m2);
1429: DMAR(c1,mul,0,nd_mod,c2);
1430: C(new) = c2;
1431: bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
1432: if ( !prev ) {
1433: prev = new;
1434: NEXT(prev) = 0;
1435: head = prev;
1436: } else {
1437: NEXT(prev) = new;
1438: NEXT(new) = 0;
1439: prev = new;
1440: }
1441: m2 = NEXT(m2);
1442: NEWNM(new);
1443: continue;
1444: }
1445: if ( cur->td > td2 )
1446: c = 1;
1447: else if ( cur->td < td2 )
1448: c = -1;
1449: else
1450: c = ndl_compare(cur->dl,m2->dl);
1451: switch ( c ) {
1452: case 0:
1453: c2 = C(m2);
1454: c1 = C(cur);
1455: DMAR(c2,mul,c1,nd_mod,t);
1456: if ( t )
1457: C(cur) = t;
1458: else if ( !prev ) {
1459: head = NEXT(cur);
1460: FREENM(cur);
1461: cur = head;
1462: } else {
1463: NEXT(prev) = NEXT(cur);
1464: FREENM(cur);
1465: cur = NEXT(prev);
1466: }
1467: m2 = NEXT(m2);
1468: break;
1469: case 1:
1470: prev = cur;
1471: cur = NEXT(cur);
1472: break;
1473: case -1:
1474: bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
1475: if ( !prev ) {
1476: /* cur = head */
1477: prev = new;
1478: c2 = C(m2);
1479: DMAR(c2,mul,0,nd_mod,c1);
1480: C(prev) = c1;
1481: NEXT(prev) = head;
1482: head = prev;
1483: } else {
1484: c2 = C(m2);
1485: DMAR(c2,mul,0,nd_mod,c1);
1486: C(new) = c1;
1487: NEXT(prev) = new;
1488: NEXT(new) = cur;
1489: prev = new;
1490: }
1491: NEWNM(new);
1492: m2 = NEXT(m2);
1493: break;
1494: }
1495: }
1496: FREENM(new);
1497: if ( head ) {
1498: BDY(p1) = head;
1499: p1->sugar = MAX(p1->sugar,p2->sugar+td);
1500: return p1;
1501: } else {
1502: FREEND(p1);
1503: return 0;
1504: }
1505:
1506: }
1507: }
1508:
1.15 ! noro 1509: INLINE int ndl_check_bound(unsigned int *d)
! 1510: {
! 1511: int i;
! 1512:
! 1513: for ( i = 0; i < nd_wpd; i++ )
! 1514: if ( d[i] & nd_mask1 )
! 1515: return 1;
! 1516: return 0;
! 1517: }
! 1518:
! 1519: int nd_sp(ND_pairs p,ND *rp)
1.11 noro 1520: {
1521: NM m;
1522: ND p1,p2,t1,t2;
1.15 ! noro 1523: unsigned int *lcm,*check;
1.11 noro 1524: int td;
1525:
1.15 ! noro 1526: check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
! 1527: p1 = nd_ps[p->i1];
! 1528: p2 = nd_ps[p->i2];
1.11 noro 1529: lcm = p->lcm;
1530: td = p->td;
1531: NEWNM(m);
1532: C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl); NEXT(m) = 0;
1.15 ! noro 1533: ndl_add(nd_bound[p->i1],m->dl,check);
! 1534: if ( ndl_check_bound(check) )
! 1535: return 0;
1.11 noro 1536: t1 = nd_mul_nm(p1,m);
1537: C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
1.15 ! noro 1538: ndl_add(nd_bound[p->i2],m->dl,check);
! 1539: if ( ndl_check_bound(check) ) {
! 1540: nd_free(t1);
! 1541: return 0;
! 1542: }
1.11 noro 1543: t2 = nd_mul_nm(p2,m);
1544: FREENM(m);
1.15 ! noro 1545: *rp = nd_add(t1,t2);
! 1546: return 1;
1.11 noro 1547: }
1548:
1.14 noro 1549: int ndl_hash_value(int td,unsigned int *d)
1550: {
1551: int i;
1552: int r;
1553:
1554: r = td;
1555: for ( i = 0; i < nd_wpd; i++ )
1556: r = ((r<<16)+d[i])%REDTAB_LEN;
1557: return r;
1558: }
1559:
1.15 ! noro 1560: int nd_find_reducer(ND g, ND *rp)
1.12 noro 1561: {
1562: NM m;
1563: ND r,p;
1.13 noro 1564: int i,c1,c2,c;
1.15 ! noro 1565: int d,k,append,index;
! 1566: unsigned int *check;
1.14 noro 1567: NM t;
1568:
1569: d = ndl_hash_value(HTD(g),HDL(g));
1570: for ( m = nd_red[d], k = 0; m; m = NEXT(m), k++ ) {
1571: if ( HTD(g) == m->td && ndl_equal(HDL(g),m->dl) ) {
1572: if ( k > 0 ) nd_notfirst++;
1.15 ! noro 1573: index = m->c;
! 1574: append = 0;
1.14 noro 1575: nd_found++;
1.15 ! noro 1576: goto found;
1.14 noro 1577: }
1578: }
1.12 noro 1579:
1580: for ( i = 0; i < nd_psn; i++ ) {
1.15 ! noro 1581: p = nd_ps[i];
1.12 noro 1582: if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
1.15 ! noro 1583: index = i;
! 1584: append = 1;
1.14 noro 1585: nd_create++;
1.15 ! noro 1586: goto found;
1.12 noro 1587: }
1588: }
1589: return 0;
1.15 ! noro 1590:
! 1591: found:
! 1592: NEWNM(m);
! 1593: p = nd_ps[index];
! 1594: ndl_sub(HDL(g),HDL(p),m->dl);
! 1595:
! 1596: check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
! 1597: ndl_add(nd_bound[index],m->dl,check);
! 1598: if ( ndl_check_bound(check) ) {
! 1599: FREENM(m);
! 1600: return -1;
! 1601: }
! 1602:
! 1603: c1 = invm(HC(p),nd_mod);
! 1604: c2 = nd_mod-HC(g);
! 1605: DMAR(c1,c2,0,nd_mod,c);
! 1606: C(m) = c;
! 1607: m->td = HTD(g)-HTD(p);
! 1608: NEXT(m) = 0;
! 1609: *rp = r = nd_mul_nm(p,m);
! 1610: FREENM(m);
! 1611:
! 1612: if ( append ) nd_append_red(HDL(g),HTD(g),i);
! 1613: return 1;
1.12 noro 1614: }
1615:
1.13 noro 1616: ND nd_find_monic_reducer(ND g)
1617: {
1618: int *d;
1619: ND p,r;
1620: int i;
1621:
1622: for ( i = 0; i < nd_psn; i++ ) {
1.15 ! noro 1623: p = nd_ps[i];
1.13 noro 1624: if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
1625: d = (int *)ALLOCA(nd_wpd*sizeof(int));
1626: ndl_sub(HDL(g),HDL(p),d);
1627: r = nd_mul_term(p,HTD(g)-HTD(p),d);
1628: return r;
1629: }
1630: }
1631: return 0;
1632: }
1633:
1.12 noro 1634: ND nd_add(ND p1,ND p2)
1.11 noro 1635: {
1.12 noro 1636: int n,c;
1637: int t;
1638: ND r;
1639: NM m1,m2,mr0,mr,s;
1640:
1641: if ( !p1 )
1642: return p2;
1643: else if ( !p2 )
1644: return p1;
1645: else {
1646: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
1647: if ( m1->td > m2->td )
1648: c = 1;
1649: else if ( m1->td < m2->td )
1650: c = -1;
1651: else
1652: c = ndl_compare(m1->dl,m2->dl);
1653: switch ( c ) {
1654: case 0:
1655: t = ((C(m1))+(C(m2))) - nd_mod;
1656: if ( t < 0 )
1657: t += nd_mod;
1658: s = m1; m1 = NEXT(m1);
1659: if ( t ) {
1660: NEXTNM2(mr0,mr,s); C(mr) = (t);
1661: } else {
1662: FREENM(s);
1663: }
1664: s = m2; m2 = NEXT(m2); FREENM(s);
1665: break;
1666: case 1:
1667: s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
1668: break;
1669: case -1:
1670: s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
1671: break;
1672: }
1673: }
1674: if ( !mr0 )
1675: if ( m1 )
1676: mr0 = m1;
1677: else if ( m2 )
1678: mr0 = m2;
1679: else
1680: return 0;
1681: else if ( m1 )
1682: NEXT(mr) = m1;
1683: else if ( m2 )
1684: NEXT(mr) = m2;
1685: else
1686: NEXT(mr) = 0;
1687: BDY(p1) = mr0;
1688: p1->sugar = MAX(p1->sugar,p2->sugar);
1689: FREEND(p2);
1690: return p1;
1691: }
1692: }
1693:
1694: ND nd_mul_nm(ND p,NM m0)
1695: {
1696: NM m,mr,mr0;
1697: unsigned int *d,*dt,*dm;
1.13 noro 1698: int c,n,td,i,c1,c2;
1.12 noro 1699: int *pt,*p1,*p2;
1700: ND r;
1.11 noro 1701:
1.12 noro 1702: if ( !p )
1.11 noro 1703: return 0;
1.12 noro 1704: else {
1705: n = NV(p); m = BDY(p);
1706: d = m0->dl; td = m0->td; c = C(m0);
1707: mr0 = 0;
1708: for ( ; m; m = NEXT(m) ) {
1709: NEXTNM(mr0,mr);
1.13 noro 1710: c1 = C(m);
1711: DMAR(c1,c,0,nd_mod,c2);
1712: C(mr) = c2;
1713: mr->td = m->td+td;
1714: ndl_add(m->dl,d,mr->dl);
1715: }
1716: NEXT(mr) = 0;
1717: MKND(NV(p),mr0,r);
1718: r->sugar = p->sugar + td;
1719: return r;
1720: }
1721: }
1722:
1723: ND nd_mul_term(ND p,int td,unsigned int *d)
1724: {
1725: NM m,mr,mr0;
1726: int c,n;
1727: ND r;
1728:
1729: if ( !p )
1730: return 0;
1731: else {
1732: n = NV(p); m = BDY(p);
1.15 ! noro 1733: for ( mr0 = 0; m; m = NEXT(m) ) {
! 1734: NEXTNM(mr0,mr);
1.13 noro 1735: C(mr) = C(m);
1.12 noro 1736: mr->td = m->td+td;
1737: ndl_add(m->dl,d,mr->dl);
1738: }
1739: NEXT(mr) = 0;
1740: MKND(NV(p),mr0,r);
1741: r->sugar = p->sugar + td;
1742: return r;
1.11 noro 1743: }
1.12 noro 1744: }
1745:
1746: #if 1
1.15 ! noro 1747: /* ret=1 : success, ret=0 : overflow */
! 1748: int nd_nf(ND g,int full,ND *rp)
1.12 noro 1749: {
1750: ND p,d,red;
1751: NM m,mrd,tail;
1.15 ! noro 1752: int n,sugar,psugar,stat;
1.12 noro 1753:
1.15 ! noro 1754: if ( !g ) {
! 1755: *rp = 0;
! 1756: return 1;
! 1757: }
1.11 noro 1758: sugar = g->sugar;
1.12 noro 1759: n = NV(g);
1.11 noro 1760: for ( d = 0; g; ) {
1.15 ! noro 1761: /* stat=1 : found, stat=0 : not found, stat=-1 : overflow */
! 1762: stat = nd_find_reducer(g,&red);
! 1763: if ( stat == 1 ) {
1.12 noro 1764: #if 1
1765: g = nd_add(g,red);
1766: sugar = MAX(sugar,red->sugar);
1.11 noro 1767: #else
1.12 noro 1768: psugar = (HTD(g)-HTD(red))+red->sugar;
1769: g = nd_reduce(g,red);
1770: sugar = MAX(sugar,psugar);
1.11 noro 1771: #endif
1.15 ! noro 1772: } else if ( stat == -1 ) {
! 1773: nd_free(g);
! 1774: nd_free(d);
! 1775: return 0;
! 1776: } else if ( !full ) {
! 1777: *rp = g;
! 1778: return 1;
! 1779: } else {
1.12 noro 1780: m = BDY(g);
1781: if ( NEXT(m) ) {
1782: BDY(g) = NEXT(m); NEXT(m) = 0;
1783: } else {
1784: FREEND(g); g = 0;
1.11 noro 1785: }
1.12 noro 1786: if ( d ) {
1787: NEXT(tail)=m;
1788: tail=m;
1.11 noro 1789: } else {
1.12 noro 1790: MKND(n,m,d);
1791: tail = BDY(d);
1.11 noro 1792: }
1793: }
1794: }
1795: if ( d )
1796: d->sugar = sugar;
1.15 ! noro 1797: *rp = d;
! 1798: return 1;
1.11 noro 1799: }
1800: #else
1801:
1802: ND nd_remove_head(ND p)
1803: {
1804: NM m;
1805:
1806: m = BDY(p);
1807: if ( !NEXT(m) ) {
1808: FREEND(p);
1809: p = 0;
1810: } else
1811: BDY(p) = NEXT(m);
1812: FREENM(m);
1813: return p;
1814: }
1815:
1816: PGeoBucket create_pbucket()
1817: {
1818: PGeoBucket g;
1819:
1820: g = CALLOC(1,sizeof(struct oPGeoBucket));
1821: g->m = -1;
1822: return g;
1823: }
1824:
1825: void add_pbucket(PGeoBucket g,ND d)
1826: {
1827: int l,k,m;
1828:
1829: l = nd_length(d);
1830: for ( k = 0, m = 1; l > m; k++, m <<= 2 );
1831: /* 4^(k-1) < l <= 4^k */
1832: d = nd_add(g->body[k],d);
1833: for ( ; d && nd_length(d) > 1<<(2*k); k++ ) {
1834: g->body[k] = 0;
1835: d = nd_add(g->body[k+1],d);
1836: }
1837: g->body[k] = d;
1838: g->m = MAX(g->m,k);
1839: }
1840:
1841: int head_pbucket(PGeoBucket g)
1842: {
1843: int j,i,c,k,nv,sum;
1844: unsigned int *di,*dj;
1845: ND gi,gj;
1846:
1847: k = g->m;
1848: while ( 1 ) {
1849: j = -1;
1850: for ( i = 0; i <= k; i++ ) {
1851: if ( !(gi = g->body[i]) )
1852: continue;
1853: if ( j < 0 ) {
1854: j = i;
1855: gj = g->body[j];
1856: dj = HDL(gj);
1857: sum = HC(gj);
1858: } else {
1859: di = HDL(gi);
1860: nv = NV(gi);
1861: if ( HTD(gi) > HTD(gj) )
1862: c = 1;
1863: else if ( HTD(gi) < HTD(gj) )
1864: c = -1;
1865: else
1866: c = ndl_compare(di,dj);
1867: if ( c > 0 ) {
1868: if ( sum )
1869: HC(gj) = sum;
1870: else
1871: g->body[j] = nd_remove_head(gj);
1872: j = i;
1873: gj = g->body[j];
1874: dj = HDL(gj);
1875: sum = HC(gj);
1876: } else if ( c == 0 ) {
1877: sum = sum+HC(gi)-nd_mod;
1878: if ( sum < 0 )
1879: sum += nd_mod;
1880: g->body[i] = nd_remove_head(gi);
1881: }
1882: }
1883: }
1884: if ( j < 0 )
1885: return -1;
1886: else if ( sum ) {
1887: HC(gj) = sum;
1888: return j;
1889: } else
1890: g->body[j] = nd_remove_head(gj);
1891: }
1892: }
1893:
1894: ND normalize_pbucket(PGeoBucket g)
1895: {
1896: int i;
1897: ND r,t;
1898:
1899: r = 0;
1900: for ( i = 0; i <= g->m; i++ )
1901: r = nd_add(r,g->body[i]);
1902: return r;
1903: }
1904:
1.12 noro 1905: ND nd_nf(ND g,int full)
1.11 noro 1906: {
1907: ND u,p,d,red;
1908: NODE l;
1909: NM m,mrd;
1910: int sugar,psugar,n,h_reducible,h;
1911: PGeoBucket bucket;
1912:
1913: if ( !g ) {
1914: return 0;
1915: }
1916: sugar = g->sugar;
1917: n = g->nv;
1918: bucket = create_pbucket();
1919: add_pbucket(bucket,g);
1920: d = 0;
1921: while ( 1 ) {
1922: h = head_pbucket(bucket);
1923: if ( h < 0 ) {
1924: if ( d )
1925: d->sugar = sugar;
1926: return d;
1927: }
1928: g = bucket->body[h];
1.12 noro 1929: red = nd_find_reducer(g);
1930: if ( red ) {
1931: bucket->body[h] = nd_remove_head(g);
1932: red = nd_remove_head(red);
1933: add_pbucket(bucket,red);
1934: sugar = MAX(sugar,red->sugar);
1935: } else if ( !full ) {
1936: g = normalize_pbucket(bucket);
1937: if ( g )
1938: g->sugar = sugar;
1939: return g;
1940: } else {
1941: m = BDY(g);
1942: if ( NEXT(m) ) {
1943: BDY(g) = NEXT(m); NEXT(m) = 0;
1944: } else {
1945: FREEND(g); g = 0;
1.11 noro 1946: }
1.12 noro 1947: bucket->body[h] = g;
1948: NEXT(m) = 0;
1949: if ( d ) {
1950: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1951: NEXT(mrd) = m;
1.11 noro 1952: } else {
1.12 noro 1953: MKND(n,m,d);
1.11 noro 1954: }
1955: }
1956: }
1957: }
1958: #endif
1959:
1960: NODE nd_gb(NODE f)
1961: {
1.15 ! noro 1962: int i,nh,sugar,stat;
1.11 noro 1963: NODE r,g,gall;
1964: ND_pairs d;
1965: ND_pairs l;
1966: ND h,nf;
1967:
1968: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1969: i = (int)BDY(r);
1970: d = update_pairs(d,g,i);
1971: g = update_base(g,i);
1972: gall = append_one(gall,i);
1973: }
1.14 noro 1974: sugar = 0;
1.11 noro 1975: while ( d ) {
1.15 ! noro 1976: again:
1.11 noro 1977: l = nd_minp(d,&d);
1.14 noro 1978: if ( l->sugar != sugar ) {
1979: sugar = l->sugar;
1980: fprintf(asir_out,"%d",sugar);
1981: }
1.15 ! noro 1982: stat = nd_sp(l,&h);
! 1983: if ( !stat ) {
! 1984: d = nd_reconstruct(d);
! 1985: goto again;
! 1986: }
! 1987: stat = nd_nf(h,!Top,&nf);
! 1988: if ( !stat ) {
! 1989: d = nd_reconstruct(d);
! 1990: goto again;
! 1991: } else if ( nf ) {
1.11 noro 1992: printf("+"); fflush(stdout);
1993: nh = nd_newps(nf);
1994: d = update_pairs(d,g,nh);
1995: g = update_base(g,nh);
1996: gall = append_one(gall,nh);
1.15 ! noro 1997: FREENDP(l);
1.11 noro 1998: } else {
1999: printf("."); fflush(stdout);
1.15 ! noro 2000: FREENDP(l);
1.11 noro 2001: }
2002: }
2003: return g;
2004: }
2005:
2006: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t)
2007: {
1.15 ! noro 2008: ND_pairs d1,nd,cur,head,prev,remove;
1.11 noro 2009:
2010: if ( !g ) return d;
2011: d = crit_B(d,t);
2012: d1 = nd_newpairs(g,t);
2013: d1 = crit_M(d1);
2014: d1 = crit_F(d1);
2015: prev = 0; cur = head = d1;
2016: while ( cur ) {
2017: if ( crit_2( cur->i1,cur->i2 ) ) {
1.15 ! noro 2018: remove = cur;
1.11 noro 2019: if ( !prev ) {
2020: head = cur = NEXT(cur);
2021: } else {
2022: cur = NEXT(prev) = NEXT(cur);
2023: }
1.15 ! noro 2024: FREENDP(remove);
1.11 noro 2025: } else {
2026: prev = cur;
2027: cur = NEXT(cur);
2028: }
2029: }
2030: if ( !d )
2031: return head;
2032: else {
2033: nd = d;
2034: while ( NEXT(nd) )
2035: nd = NEXT(nd);
2036: NEXT(nd) = head;
2037: return d;
2038: }
2039: }
2040:
2041: ND_pairs nd_newpairs( NODE g, int t )
2042: {
2043: NODE h;
2044: unsigned int *dl;
2045: int td,ts,s;
2046: ND_pairs r,r0;
2047:
1.15 ! noro 2048: dl = HDL(nd_ps[t]);
! 2049: td = HTD(nd_ps[t]);
! 2050: ts = nd_ps[t]->sugar - td;
1.11 noro 2051: for ( r0 = 0, h = g; h; h = NEXT(h) ) {
2052: NEXTND_pairs(r0,r);
2053: r->i1 = (int)BDY(h);
2054: r->i2 = t;
1.15 ! noro 2055: ndl_lcm(HDL(nd_ps[r->i1]),dl,r->lcm);
1.11 noro 2056: r->td = ndl_td(r->lcm);
1.15 ! noro 2057: s = nd_ps[r->i1]->sugar-HTD(nd_ps[r->i1]);
1.11 noro 2058: r->sugar = MAX(s,ts) + r->td;
2059: }
2060: NEXT(r) = 0;
2061: return r0;
2062: }
2063:
2064: ND_pairs crit_B( ND_pairs d, int s )
2065: {
1.15 ! noro 2066: ND_pairs cur,head,prev,remove;
1.11 noro 2067: unsigned int *t,*tl,*lcm;
2068: int td,tdl;
2069:
2070: if ( !d ) return 0;
1.15 ! noro 2071: t = HDL(nd_ps[s]);
1.11 noro 2072: prev = 0;
2073: head = cur = d;
2074: lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
2075: while ( cur ) {
2076: tl = cur->lcm;
2077: if ( ndl_reducible(tl,t)
1.15 ! noro 2078: && (ndl_lcm(HDL(nd_ps[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
! 2079: && (ndl_lcm(HDL(nd_ps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
! 2080: remove = cur;
1.11 noro 2081: if ( !prev ) {
2082: head = cur = NEXT(cur);
2083: } else {
2084: cur = NEXT(prev) = NEXT(cur);
2085: }
1.15 ! noro 2086: FREENDP(remove);
1.11 noro 2087: } else {
2088: prev = cur;
2089: cur = NEXT(cur);
2090: }
2091: }
2092: return head;
2093: }
2094:
2095: ND_pairs crit_M( ND_pairs d1 )
2096: {
2097: ND_pairs e,d2,d3,dd,p;
2098: unsigned int *id,*jd;
2099: int itd,jtd;
2100:
2101: for ( dd = 0, e = d1; e; e = d3 ) {
2102: if ( !(d2 = NEXT(e)) ) {
2103: NEXT(e) = dd;
2104: return e;
2105: }
2106: id = e->lcm;
2107: itd = e->td;
2108: for ( d3 = 0; d2; d2 = p ) {
2109: p = NEXT(d2),
2110: jd = d2->lcm;
2111: jtd = d2->td;
2112: if ( jtd == itd )
2113: if ( id == jd );
2114: else if ( ndl_reducible(jd,id) ) continue;
2115: else if ( ndl_reducible(id,jd) ) goto delit;
2116: else ;
2117: else if ( jtd > itd )
2118: if ( ndl_reducible(jd,id) ) continue;
2119: else ;
2120: else if ( ndl_reducible(id,jd ) ) goto delit;
2121: NEXT(d2) = d3;
2122: d3 = d2;
2123: }
2124: NEXT(e) = dd;
2125: dd = e;
2126: continue;
2127: /**/
2128: delit: NEXT(d2) = d3;
2129: d3 = d2;
2130: for ( ; p; p = d2 ) {
2131: d2 = NEXT(p);
2132: NEXT(p) = d3;
2133: d3 = p;
2134: }
1.15 ! noro 2135: FREENDP(e);
1.11 noro 2136: }
2137: return dd;
2138: }
2139:
2140: ND_pairs crit_F( ND_pairs d1 )
2141: {
1.15 ! noro 2142: ND_pairs rest, head,remove;
1.11 noro 2143: ND_pairs last, p, r, w;
2144: int s;
2145:
2146: for ( head = last = 0, p = d1; NEXT(p); ) {
2147: r = w = equivalent_pairs(p,&rest);
2148: s = r->sugar;
1.15 ! noro 2149: w = NEXT(w);
! 2150: while ( w ) {
1.11 noro 2151: if ( crit_2(w->i1,w->i2) ) {
2152: r = w;
1.15 ! noro 2153: w = NEXT(w);
! 2154: while ( w ) {
! 2155: remove = w;
! 2156: w = NEXT(w);
! 2157: FREENDP(remove);
! 2158: }
1.11 noro 2159: break;
2160: } else if ( w->sugar < s ) {
1.15 ! noro 2161: FREENDP(r);
1.11 noro 2162: r = w;
2163: s = r->sugar;
1.15 ! noro 2164: w = NEXT(w);
! 2165: } else {
! 2166: remove = w;
! 2167: w = NEXT(w);
! 2168: FREENDP(remove);
1.11 noro 2169: }
1.15 ! noro 2170: }
1.11 noro 2171: if ( last ) NEXT(last) = r;
2172: else head = r;
2173: NEXT(last = r) = 0;
2174: p = rest;
2175: if ( !p ) return head;
2176: }
2177: if ( !last ) return p;
2178: NEXT(last) = p;
2179: return head;
2180: }
2181:
2182: int crit_2( int dp1, int dp2 )
2183: {
1.15 ! noro 2184: return ndl_disjoint(HDL(nd_ps[dp1]),HDL(nd_ps[dp2]));
1.11 noro 2185: }
2186:
2187: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest )
2188: {
2189: ND_pairs w,p,r,s;
2190: unsigned int *d;
2191: int td;
2192:
2193: w = d1;
2194: d = w->lcm;
2195: td = w->td;
2196: s = NEXT(w);
2197: NEXT(w) = 0;
2198: for ( r = 0; s; s = p ) {
2199: p = NEXT(s);
2200: if ( td == s->td && ndl_equal(d,s->lcm) ) {
2201: NEXT(s) = w;
2202: w = s;
2203: } else {
2204: NEXT(s) = r;
2205: r = s;
2206: }
2207: }
2208: *prest = r;
2209: return w;
2210: }
2211:
2212: NODE update_base(NODE nd,int ndp)
2213: {
2214: unsigned int *dl, *dln;
2215: NODE last, p, head;
2216: int td,tdn;
2217:
1.15 ! noro 2218: dl = HDL(nd_ps[ndp]);
! 2219: td = HTD(nd_ps[ndp]);
1.11 noro 2220: for ( head = last = 0, p = nd; p; ) {
1.15 ! noro 2221: dln = HDL(nd_ps[(int)BDY(p)]);
! 2222: tdn = HTD(nd_ps[(int)BDY(p)]);
1.11 noro 2223: if ( tdn >= td && ndl_reducible( dln, dl ) ) {
2224: p = NEXT(p);
2225: if ( last ) NEXT(last) = p;
2226: } else {
2227: if ( !last ) head = p;
2228: p = NEXT(last = p);
2229: }
2230: }
2231: head = append_one(head,ndp);
2232: return head;
2233: }
2234:
2235: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
2236: {
2237: ND_pairs m,ml,p,l;
2238: unsigned int *lcm;
2239: int s,td,len,tlen,c;
2240:
2241: if ( !(p = NEXT(m = d)) ) {
2242: *prest = p;
2243: NEXT(m) = 0;
2244: return m;
2245: }
2246: lcm = m->lcm;
2247: s = m->sugar;
2248: td = m->td;
1.15 ! noro 2249: len = nd_length(nd_ps[m->i1])+nd_length(nd_ps[m->i2]);
1.11 noro 2250: for ( ml = 0, l = m; p; p = NEXT(l = p) ) {
2251: if (p->sugar < s)
2252: goto find;
2253: else if ( p->sugar == s ) {
2254: if ( p->td < td )
2255: goto find;
2256: else if ( p->td == td ) {
2257: c = ndl_compare(p->lcm,lcm);
2258: if ( c < 0 )
2259: goto find;
2260: else if ( c == 0 ) {
1.15 ! noro 2261: tlen = nd_length(nd_ps[p->i1])+nd_length(nd_ps[p->i2]);
1.11 noro 2262: if ( tlen < len )
2263: goto find;
2264: }
2265: }
2266: }
2267: continue;
2268: find:
2269: ml = l;
2270: m = p;
2271: lcm = m->lcm;
2272: s = m->sugar;
2273: td = m->td;
2274: len = tlen;
2275: }
2276: if ( !ml ) *prest = NEXT(m);
2277: else {
2278: NEXT(ml) = NEXT(m);
2279: *prest = d;
2280: }
2281: NEXT(m) = 0;
2282: return m;
2283: }
2284:
2285: int nd_newps(ND a)
2286: {
2287: if ( nd_psn == nd_pslen ) {
2288: nd_pslen *= 2;
1.15 ! noro 2289: nd_ps = (ND *)REALLOC((char *)nd_ps,nd_pslen*sizeof(ND));
! 2290: nd_bound = (unsigned int **)
! 2291: REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
1.11 noro 2292: }
2293: nd_monic(a);
1.15 ! noro 2294: nd_ps[nd_psn] = a;
! 2295: nd_bound[nd_psn] = nd_compute_bound(a);
1.11 noro 2296: return nd_psn++;
2297: }
2298:
2299: NODE NODE_sortb(NODE f,int);
2300: ND dptond(DP);
2301: DP ndtodp(ND);
2302:
2303: NODE nd_setup(NODE f)
2304: {
1.14 noro 2305: int i,td;
1.11 noro 2306: NODE s,s0,f0;
2307:
1.14 noro 2308: nd_found = 0;
2309: nd_notfirst = 0;
2310: nd_create = 0;
1.11 noro 2311: #if 0
2312: f0 = f = NODE_sortb(f,1);
2313: #endif
2314: nd_psn = length(f); nd_pslen = 2*nd_psn;
1.15 ! noro 2315: nd_ps = (ND *)MALLOC(nd_pslen*sizeof(ND));
! 2316: nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
! 2317: nd_bpe = 4;
! 2318: nd_setup_parameters();
1.11 noro 2319: nd_free_private_storage();
2320: for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
1.15 ! noro 2321: nd_ps[i] = dptond((DP)BDY(f));
! 2322: nd_monic(nd_ps[i]);
! 2323: nd_bound[i] = nd_compute_bound(nd_ps[i]);
1.11 noro 2324: }
1.14 noro 2325: nd_red = (NM *)MALLOC(REDTAB_LEN*sizeof(NM));
1.11 noro 2326: for ( s0 = 0, i = 0; i < nd_psn; i++ ) {
2327: NEXTNODE(s0,s); BDY(s) = (pointer)i;
2328: }
2329: if ( s0 ) NEXT(s) = 0;
2330: return s0;
2331: }
2332:
2333: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
2334: {
2335: struct order_spec ord1;
2336: VL fv,vv,vc;
2337: NODE fd,fd0,r,r0,t,x,s,xx;
2338: DP a,b,c;
2339:
2340: get_vars((Obj)f,&fv); pltovl(v,&vv);
2341: nd_nvar = length(vv);
2342: if ( ord->id )
2343: error("nd_gr : unsupported order");
2344: switch ( ord->ord.simple ) {
2345: case 0:
2346: is_rlex = 1;
2347: break;
2348: case 1:
2349: is_rlex = 0;
2350: break;
2351: default:
2352: error("nd_gr : unsupported order");
2353: }
2354: initd(ord);
2355: nd_mod = m;
2356: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
2357: ptod(CO,vv,(P)BDY(t),&b);
2358: _dp_mod(b,m,0,&c);
2359: if ( c ) {
2360: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
2361: }
2362: }
2363: if ( fd0 ) NEXT(fd) = 0;
2364: s = nd_setup(fd0);
2365: x = nd_gb(s);
2366: #if 0
2367: x = nd_reduceall(x,m);
2368: #endif
2369: for ( r0 = 0; x; x = NEXT(x) ) {
2370: NEXTNODE(r0,r);
1.15 ! noro 2371: a = ndtodp(nd_ps[(int)BDY(x)]);
1.11 noro 2372: _dtop_mod(CO,vv,a,(P *)&BDY(r));
2373: }
2374: if ( r0 ) NEXT(r) = 0;
2375: MKLIST(*rp,r0);
1.14 noro 2376: fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
2377: nd_found,nd_notfirst,nd_create);
1.11 noro 2378: }
2379:
2380: void dltondl(int n,DL dl,unsigned int *r)
2381: {
2382: unsigned int *d;
2383: int i;
2384:
2385: d = dl->d;
1.12 noro 2386: bzero(r,nd_wpd*sizeof(unsigned int));
1.11 noro 2387: if ( is_rlex )
2388: for ( i = 0; i < n; i++ )
2389: r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe));
2390: else
2391: for ( i = 0; i < n; i++ )
2392: r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe);
2393: }
2394:
2395: DL ndltodl(int n,int td,unsigned int *ndl)
2396: {
2397: DL dl;
2398: int *d;
2399: int i;
2400:
2401: NEWDL(dl,n);
2402: dl->td = td;
2403: d = dl->d;
2404: if ( is_rlex )
2405: for ( i = 0; i < n; i++ )
2406: d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
2407: &((1<<nd_bpe)-1);
2408: else
2409: for ( i = 0; i < n; i++ )
2410: d[i] = (ndl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
2411: &((1<<nd_bpe)-1);
2412: return dl;
2413: }
2414:
2415: ND dptond(DP p)
2416: {
2417: ND d;
2418: NM m0,m;
2419: MP t;
2420: int n;
2421:
2422: if ( !p )
2423: return 0;
2424: n = NV(p);
2425: m0 = 0;
2426: for ( t = BDY(p); t; t = NEXT(t) ) {
2427: NEXTNM(m0,m);
2428: m->c = ITOS(t->c);
2429: m->td = t->dl->td;
2430: dltondl(n,t->dl,m->dl);
2431: }
2432: NEXT(m) = 0;
2433: MKND(n,m0,d);
2434: d->nv = n;
2435: d->sugar = p->sugar;
2436: return d;
2437: }
2438:
2439: DP ndtodp(ND p)
2440: {
2441: DP d;
2442: MP m0,m;
2443: NM t;
2444: int n;
2445:
2446: if ( !p )
2447: return 0;
2448: n = NV(p);
2449: m0 = 0;
2450: for ( t = BDY(p); t; t = NEXT(t) ) {
2451: NEXTMP(m0,m);
2452: m->c = STOI(t->c);
2453: m->dl = ndltodl(n,t->td,t->dl);
2454: }
2455: NEXT(m) = 0;
2456: MKDP(n,m0,d);
2457: d->sugar = p->sugar;
2458: return d;
2459: }
2460:
2461: void ndl_print(unsigned int *dl)
2462: {
2463: int n;
2464: int i;
2465:
2466: n = nd_nvar;
2467: printf("<<");
2468: if ( is_rlex )
2469: for ( i = 0; i < n; i++ )
2470: printf(i==n-1?"%d":"%d,",
2471: (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
2472: &((1<<nd_bpe)-1));
2473: else
2474: for ( i = 0; i < n; i++ )
2475: printf(i==n-1?"%d":"%d,",
2476: (dl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
2477: &((1<<nd_bpe)-1));
2478: printf(">>");
2479: }
2480:
2481: void nd_print(ND p)
2482: {
2483: NM m;
2484:
2485: if ( !p )
2486: printf("0\n");
2487: else {
2488: for ( m = BDY(p); m; m = NEXT(m) ) {
2489: printf("+%d*",m->c);
2490: ndl_print(m->dl);
2491: }
2492: printf("\n");
2493: }
2494: }
2495:
2496: void ndp_print(ND_pairs d)
2497: {
2498: ND_pairs t;
2499:
2500: for ( t = d; t; t = NEXT(t) ) {
2501: printf("%d,%d ",t->i1,t->i2);
2502: }
2503: printf("\n");
2504: }
2505:
2506: void nd_monic(ND p)
2507: {
1.13 noro 2508: if ( !p )
2509: return;
2510: else
2511: nd_mul_c(p,invm(HC(p),nd_mod));
2512: }
2513:
2514: void nd_mul_c(ND p,int mul)
2515: {
1.11 noro 2516: NM m;
1.13 noro 2517: int c,c1;
1.11 noro 2518:
2519: if ( !p )
2520: return;
1.13 noro 2521: for ( m = BDY(p); m; m = NEXT(m) ) {
2522: c1 = C(m);
2523: DMAR(c1,mul,0,nd_mod,c);
2524: C(m) = c;
2525: }
1.14 noro 2526: }
2527:
2528: void nd_free(ND p)
2529: {
2530: NM t,s;
2531:
2532: if ( !p )
2533: return;
2534: t = BDY(p);
2535: while ( t ) {
2536: s = NEXT(t);
2537: FREENM(t);
2538: t = s;
2539: }
2540: FREEND(p);
2541: }
2542:
2543: void nd_append_red(unsigned int *d,int td,int i)
2544: {
2545: NM m,m0;
2546: int h;
2547:
2548: NEWNM(m);
2549: h = ndl_hash_value(td,d);
2550: m->c = i;
2551: m->td = td;
2552: bcopy(d,m->dl,nd_wpd*sizeof(unsigned int));
2553: NEXT(m) = nd_red[h];
2554: nd_red[h] = m;
1.15 ! noro 2555: }
! 2556:
! 2557: unsigned int *nd_compute_bound(ND p)
! 2558: {
! 2559: unsigned int *d1,*d2,*t;
! 2560: NM m;
! 2561:
! 2562: if ( !p )
! 2563: return 0;
! 2564: d1 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
! 2565: d2 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
! 2566: bcopy(HDL(p),d1,nd_wpd*sizeof(unsigned int));
! 2567: for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
! 2568: ndl_lcm(m->dl,d1,d2);
! 2569: t = d1; d1 = d2; d2 = t;
! 2570: }
! 2571: t = (unsigned int *)MALLOC_ATOMIC(nd_wpd*sizeof(unsigned int));
! 2572: bcopy(d1,t,nd_wpd*sizeof(unsigned int));
! 2573: return t;
! 2574: }
! 2575:
! 2576: void nd_setup_parameters() {
! 2577: int i;
! 2578:
! 2579: nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
! 2580: nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
! 2581: if ( nd_bpe < 32 ) {
! 2582: nd_mask0 = (1<<nd_bpe)-1;
! 2583: } else {
! 2584: nd_mask0 = 0xffffffff;
! 2585: }
! 2586: bzero(nd_mask,sizeof(nd_mask));
! 2587: nd_mask1 = 0;
! 2588: for ( i = 0; i < nd_epw; i++ ) {
! 2589: nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
! 2590: nd_mask1 |= (1<<(nd_bpe-1))<<(i*nd_bpe);
! 2591: }
! 2592: }
! 2593:
! 2594: ND_pairs nd_reconstruct(ND_pairs d)
! 2595: {
! 2596: int i,obpe;
! 2597: NM prev_nm_free_list;
! 2598: ND_pairs s0,s,t,prev_ndp_free_list;
! 2599:
! 2600: obpe = nd_bpe;
! 2601: switch ( nd_bpe ) {
! 2602: case 4: nd_bpe = 6; break;
! 2603: case 6: nd_bpe = 8; break;
! 2604: case 8: nd_bpe = 16; break;
! 2605: case 16: nd_bpe = 32; break;
! 2606: }
! 2607: nd_setup_parameters();
! 2608: prev_nm_free_list = _nm_free_list;
! 2609: prev_ndp_free_list = _ndp_free_list;
! 2610: _nm_free_list = 0;
! 2611: _ndp_free_list = 0;
! 2612: for ( i = 0; i < nd_psn; i++ ) {
! 2613: nd_ps[i] = nd_dup(nd_ps[i],obpe);
! 2614: nd_bound[i] = nd_compute_bound(nd_ps[i]);
! 2615: }
! 2616: s0 = 0;
! 2617: for ( t = d; t; t = NEXT(t) ) {
! 2618: NEXTND_pairs(s0,s);
! 2619: s->i1 = t->i1;
! 2620: s->i2 = t->i2;
! 2621: s->td = t->td;
! 2622: s->sugar = t->sugar;
! 2623: ndl_dup(obpe,t->lcm,s->lcm);
! 2624: }
! 2625: if ( s0 ) NEXT(s) = 0;
! 2626: prev_nm_free_list = 0;
! 2627: prev_ndp_free_list = 0;
! 2628: GC_gcollect();
! 2629: return s0;
! 2630: }
! 2631:
! 2632: void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
! 2633: {
! 2634: int n,i,ei,oepw,cepw,cbpe;
! 2635:
! 2636: n = nd_nvar;
! 2637: oepw = (sizeof(unsigned int)*8)/obpe;
! 2638: cepw = nd_epw;
! 2639: cbpe = nd_bpe;
! 2640: if ( is_rlex )
! 2641: for ( i = 0; i < n; i++ ) {
! 2642: ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe))
! 2643: &((1<<obpe)-1);
! 2644: r[(n-1-i)/cepw] |= (ei<<((cepw-((n-1-i)%cepw)-1)*cbpe));
! 2645: }
! 2646: else
! 2647: for ( i = 0; i < n; i++ ) {
! 2648: ei = (d[i/oepw]>>((oepw-(i%oepw)-1)*obpe))
! 2649: &((1<<obpe)-1);
! 2650: r[i/cepw] |= (ei<<((cepw-(i%cepw)-1)*cbpe));
! 2651: }
! 2652: }
! 2653:
! 2654: ND nd_dup(ND p,int obpe)
! 2655: {
! 2656: NM m,mr,mr0;
! 2657: int c,n;
! 2658: ND r;
! 2659:
! 2660: if ( !p )
! 2661: return 0;
! 2662: else {
! 2663: n = NV(p); m = BDY(p);
! 2664: for ( mr0 = 0; m; m = NEXT(m) ) {
! 2665: NEXTNM(mr0,mr);
! 2666: C(mr) = C(m);
! 2667: mr->td = m->td;
! 2668: ndl_dup(obpe,m->dl,mr->dl);
! 2669: }
! 2670: NEXT(mr) = 0;
! 2671: MKND(NV(p),mr0,r);
! 2672: r->sugar = p->sugar;
! 2673: return r;
! 2674: }
1.1 noro 2675: }