![[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.11 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.10 2002/01/28 00:54:43 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:
! 899: typedef struct oPGeoBucket {
! 900: int m;
! 901: struct oND *body[32];
! 902: } *PGeoBucket;
! 903:
! 904: typedef struct oND {
! 905: struct oNM *body;
! 906: int nv;
! 907: int sugar;
! 908: } *ND;
! 909:
! 910: typedef struct oNM {
! 911: struct oNM *next;
! 912: int td;
! 913: int c;
! 914: unsigned int dl[1];
! 915: } *NM;
! 916:
! 917: typedef struct oND_pairs {
! 918: struct oND_pairs *next;
! 919: int i1,i2;
! 920: int td,sugar;
! 921: unsigned int lcm[1];
! 922: } *ND_pairs;
! 923:
! 924: static ND *nps;
! 925: int nd_mod,nd_nvar;
! 926: int is_rlex;
! 927: int nd_epw,nd_bpe,nd_wpd;
! 928: NM _nm_free_list;
! 929: ND _nd_free_list;
! 930: ND_pairs _ndp_free_list;
! 931:
! 932: extern int Top,Reverse;
! 933: int nd_psn,nd_pslen;
! 934:
! 935: void GC_gcollect();
! 936: NODE append_one(NODE,int);
! 937:
! 938: #define HTD(d) ((d)->body->td)
! 939: #define HDL(d) ((d)->body->dl)
! 940: #define HC(d) ((d)->body->c)
! 941:
! 942: #define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list)
! 943: #define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list)
! 944: #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)
! 945:
! 946: #define NEXTNM(r,c) \
! 947: if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);}
! 948: #define NEXTNM2(r,c,s) \
! 949: if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
! 950: #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)
! 951: #define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m)
! 952:
! 953: #define NEXTND_pairs(r,c) \
! 954: if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
! 955:
! 956: ND_pairs crit_B( ND_pairs d, int s );
! 957: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp);
! 958: NODE nd_setup(NODE f);
! 959: int nd_newps(ND a);
! 960: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest );
! 961: NODE update_base(NODE nd,int ndp);
! 962: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest );
! 963: int crit_2( int dp1, int dp2 );
! 964: ND_pairs crit_F( ND_pairs d1 );
! 965: ND_pairs crit_M( ND_pairs d1 );
! 966: ND_pairs nd_newpairs( NODE g, int t );
! 967: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t);
! 968: NODE nd_gb(NODE f);
! 969: void nd_free_private_storage();
! 970: void _NM_alloc();
! 971: void _ND_alloc();
! 972: int ndl_td(unsigned int *d);
! 973: ND nd_add(ND p1,ND p2);
! 974: ND nd_mul_nm(ND p,NM m0);
! 975: ND nd_sp(ND_pairs p);
! 976: ND nd_reducer(ND p1,ND p2);
! 977: ND nd_nf(NODE b,ND g,ND *ps,int full);
! 978: void ndl_print(unsigned int *dl);
! 979: void nd_print(ND p);
! 980: void ndp_print(ND_pairs d);
! 981: int nd_length(ND p);
! 982: void nd_monic(ND p);
! 983:
! 984: void nd_free_private_storage()
! 985: {
! 986: _nd_free_list = 0;
! 987: _nm_free_list = 0;
! 988: GC_gcollect();
! 989: }
! 990:
! 991: void _NM_alloc()
! 992: {
! 993: NM p;
! 994: int i;
! 995:
! 996: for ( i = 0; i < 1024; i++ ) {
! 997: p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
! 998: p->next = _nm_free_list; _nm_free_list = p;
! 999: }
! 1000: }
! 1001:
! 1002: void _ND_alloc()
! 1003: {
! 1004: ND p;
! 1005: int i;
! 1006:
! 1007: for ( i = 0; i < 1024; i++ ) {
! 1008: p = (ND)GC_malloc(sizeof(struct oND));
! 1009: p->body = (NM)_nd_free_list; _nd_free_list = p;
! 1010: }
! 1011: }
! 1012:
! 1013: void _NDP_alloc()
! 1014: {
! 1015: ND_pairs p;
! 1016: int i;
! 1017:
! 1018: for ( i = 0; i < 1024; i++ ) {
! 1019: p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs)
! 1020: +(nd_wpd-1)*sizeof(unsigned int));
! 1021: p->next = _ndp_free_list; _ndp_free_list = p;
! 1022: }
! 1023: }
! 1024:
! 1025: INLINE nd_length(ND p)
! 1026: {
! 1027: NM m;
! 1028: int i;
! 1029:
! 1030: if ( !p )
! 1031: return 0;
! 1032: else {
! 1033: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
! 1034: return i;
! 1035: }
! 1036: }
! 1037:
! 1038: INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2)
! 1039: {
! 1040: unsigned int u1,u2;
! 1041: int i;
! 1042:
! 1043: switch ( nd_bpe ) {
! 1044: case 4:
! 1045: for ( i = 0; i < nd_wpd; i++ ) {
! 1046: u1 = d1[i]; u2 = d2[i];
! 1047: if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0;
! 1048: if ( (u1&0xf000000) < (u2&0xf000000) ) return 0;
! 1049: if ( (u1&0xf00000) < (u2&0xf00000) ) return 0;
! 1050: if ( (u1&0xf0000) < (u2&0xf0000) ) return 0;
! 1051: if ( (u1&0xf000) < (u2&0xf000) ) return 0;
! 1052: if ( (u1&0xf00) < (u2&0xf00) ) return 0;
! 1053: if ( (u1&0xf0) < (u2&0xf0) ) return 0;
! 1054: if ( (u1&0xf) < (u2&0xf) ) return 0;
! 1055: }
! 1056: return 1;
! 1057: break;
! 1058: case 8:
! 1059: for ( i = 0; i < nd_wpd; i++ ) {
! 1060: u1 = d1[i]; u2 = d2[i];
! 1061: if ( (u1&0xff000000) < (u2&0xff000000) ) return 0;
! 1062: if ( (u1&0xff0000) < (u2&0xff0000) ) return 0;
! 1063: if ( (u1&0xff00) < (u2&0xff00) ) return 0;
! 1064: if ( (u1&0xff) < (u2&0xff) ) return 0;
! 1065: }
! 1066: return 1;
! 1067: break;
! 1068: case 16:
! 1069: for ( i = 0; i < nd_wpd; i++ ) {
! 1070: u1 = d1[i]; u2 = d2[i];
! 1071: if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0;
! 1072: if ( (u1&0xffff) < (u2&0xffff) ) return 0;
! 1073: }
! 1074: return 1;
! 1075: break;
! 1076: case 32:
! 1077: for ( i = 0; i < nd_wpd; i++ )
! 1078: if ( d1[i] < d2[i] ) return 0;
! 1079: return 1;
! 1080: break;
! 1081: default:
! 1082: error("ndl_reducible : not implemented yet");
! 1083: }
! 1084: }
! 1085:
! 1086: INLINE void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
! 1087: {
! 1088: unsigned int t1,t2,u,u1,u2;
! 1089: int i;
! 1090:
! 1091: switch ( nd_bpe ) {
! 1092: case 4:
! 1093: for ( i = 0; i < nd_wpd; i++ ) {
! 1094: u1 = d1[i]; u2 = d2[i];
! 1095: t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2;
! 1096: t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2;
! 1097: t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2;
! 1098: t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2;
! 1099: t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2;
! 1100: t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2;
! 1101: t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2;
! 1102: t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2;
! 1103: d[i] = u;
! 1104: }
! 1105: break;
! 1106: case 8:
! 1107: for ( i = 0; i < nd_wpd; i++ ) {
! 1108: u1 = d1[i]; u2 = d2[i];
! 1109: t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2;
! 1110: t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2;
! 1111: t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2;
! 1112: t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2;
! 1113: d[i] = u;
! 1114: }
! 1115: break;
! 1116: case 16:
! 1117: for ( i = 0; i < nd_wpd; i++ ) {
! 1118: u1 = d1[i]; u2 = d2[i];
! 1119: t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2;
! 1120: t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2;
! 1121: d[i] = u;
! 1122: }
! 1123: break;
! 1124: case 32:
! 1125: for ( i = 0; i < nd_wpd; i++ ) {
! 1126: u1 = d1[i]; u2 = d2[i];
! 1127: d[i] = u1>u2?u1:u2;
! 1128: }
! 1129: break;
! 1130: default:
! 1131: error("ndl_lcm : not implemented yet");
! 1132: }
! 1133: }
! 1134:
! 1135: INLINE int ndl_td(unsigned int *d)
! 1136: {
! 1137: unsigned int t,u;
! 1138: int i;
! 1139:
! 1140: switch ( nd_bpe ) {
! 1141: case 4:
! 1142: for ( t = 0, i = 0; i < nd_wpd; i++ ) {
! 1143: u = d[i];
! 1144: t += ((u&0xf0000000)>>28)+((u&0xf000000)>>24)
! 1145: +((u&0xf00000)>>20)+((u&0xf0000)>>16)
! 1146: +((u&0xf000)>>12)+((u&0xf00)>>8)+((u&0xf0)>>4)+(u&0xf);
! 1147: }
! 1148: break;
! 1149: case 8:
! 1150: for ( t = 0, i = 0; i < nd_wpd; i++ ) {
! 1151: u = d[i];
! 1152: t += ((u&0xff000000)>>24)+((u&0xff0000)>>16)
! 1153: +((u&0xff00)>>8)+(u&0xff);
! 1154: }
! 1155: break;
! 1156: case 16:
! 1157: for ( t = 0, i = 0; i < nd_wpd; i++ ) {
! 1158: u = d[i];
! 1159: t += ((u&0xffff0000)>>16)+(u&0xffff);
! 1160: }
! 1161: break;
! 1162: case 32:
! 1163: for ( t = 0, i = 0; i < nd_wpd; i++ )
! 1164: t += d[i];
! 1165: break;
! 1166: default:
! 1167: error("ndl_td : not implemented yet");
! 1168: }
! 1169: return t;
! 1170: }
! 1171:
! 1172: INLINE int ndl_compare(unsigned int *d1,unsigned int *d2)
! 1173: {
! 1174: int i;
! 1175:
! 1176: for ( i = 0; i < nd_wpd; i++ )
! 1177: if ( d1[i] > d2[i] )
! 1178: return is_rlex ? -1 : 1;
! 1179: else if ( d1[i] < d2[i] )
! 1180: return is_rlex ? 1 : -1;
! 1181: return 0;
! 1182: }
! 1183:
! 1184: INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)
! 1185: {
! 1186: int i;
! 1187:
! 1188: for ( i = 0; i < nd_wpd; i++ )
! 1189: if ( d1[i] != d2[i] )
! 1190: return 0;
! 1191: return 1;
! 1192: }
! 1193:
! 1194: INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d)
! 1195: {
! 1196: int i;
! 1197:
! 1198: for ( i = 0; i < nd_wpd; i++ )
! 1199: d[i] = d1[i]+d2[i];
! 1200: }
! 1201:
! 1202: INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
! 1203: {
! 1204: int i;
! 1205:
! 1206: for ( i = 0; i < nd_wpd; i++ )
! 1207: d[i] = d1[i]-d2[i];
! 1208: }
! 1209:
! 1210: INLINE int ndl_disjoint(unsigned int *d1,unsigned int *d2)
! 1211: {
! 1212: unsigned int t1,t2,u,u1,u2;
! 1213: int i;
! 1214:
! 1215: switch ( nd_bpe ) {
! 1216: case 4:
! 1217: for ( i = 0; i < nd_wpd; i++ ) {
! 1218: u1 = d1[i]; u2 = d2[i];
! 1219: t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0;
! 1220: t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0;
! 1221: t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0;
! 1222: t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0;
! 1223: t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0;
! 1224: t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0;
! 1225: t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0;
! 1226: t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0;
! 1227: }
! 1228: return 1;
! 1229: break;
! 1230: case 8:
! 1231: for ( i = 0; i < nd_wpd; i++ ) {
! 1232: u1 = d1[i]; u2 = d2[i];
! 1233: t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0;
! 1234: t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0;
! 1235: t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0;
! 1236: t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0;
! 1237: }
! 1238: return 1;
! 1239: break;
! 1240: case 16:
! 1241: for ( i = 0; i < nd_wpd; i++ ) {
! 1242: u1 = d1[i]; u2 = d2[i];
! 1243: t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0;
! 1244: t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0;
! 1245: }
! 1246: return 1;
! 1247: break;
! 1248: case 32:
! 1249: for ( i = 0; i < nd_wpd; i++ )
! 1250: if ( d1[i] && d2[i] ) return 0;
! 1251: return 1;
! 1252: break;
! 1253: default:
! 1254: error("ndl_disjoint : not implemented yet");
! 1255: }
! 1256: }
! 1257:
! 1258: ND nd_add(ND p1,ND p2)
! 1259: {
! 1260: int n,c;
! 1261: int t;
! 1262: ND r;
! 1263: NM m1,m2,mr0,mr,s;
! 1264:
! 1265: if ( !p1 )
! 1266: return p2;
! 1267: else if ( !p2 )
! 1268: return p1;
! 1269: else {
! 1270: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
! 1271: if ( m1->td > m2->td )
! 1272: c = 1;
! 1273: else if ( m1->td < m2->td )
! 1274: c = -1;
! 1275: else
! 1276: c = ndl_compare(m1->dl,m2->dl);
! 1277: switch ( c ) {
! 1278: case 0:
! 1279: t = ((C(m1))+(C(m2))) - nd_mod;
! 1280: if ( t < 0 )
! 1281: t += nd_mod;
! 1282: s = m1; m1 = NEXT(m1);
! 1283: if ( t ) {
! 1284: NEXTNM2(mr0,mr,s); C(mr) = (t);
! 1285: } else {
! 1286: FREENM(s);
! 1287: }
! 1288: s = m2; m2 = NEXT(m2); FREENM(s);
! 1289: break;
! 1290: case 1:
! 1291: s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
! 1292: break;
! 1293: case -1:
! 1294: s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
! 1295: break;
! 1296: }
! 1297: }
! 1298: if ( !mr0 )
! 1299: if ( m1 )
! 1300: mr0 = m1;
! 1301: else if ( m2 )
! 1302: mr0 = m2;
! 1303: else
! 1304: return 0;
! 1305: else if ( m1 )
! 1306: NEXT(mr) = m1;
! 1307: else if ( m2 )
! 1308: NEXT(mr) = m2;
! 1309: else
! 1310: NEXT(mr) = 0;
! 1311: BDY(p1) = mr0;
! 1312: p1->sugar = MAX(p1->sugar,p2->sugar);
! 1313: FREEND(p2);
! 1314: return p1;
! 1315: }
! 1316: }
! 1317:
! 1318: INLINE ND nd_mul_nm(ND p,NM m0)
! 1319: {
! 1320: NM m,mr,mr0;
! 1321: unsigned int *d,*dt,*dm;
! 1322: int c,n,td;
! 1323: int *pt,*p1,*p2;
! 1324: ND r;
! 1325:
! 1326: if ( !p )
! 1327: return 0;
! 1328: else {
! 1329: n = NV(p); m = BDY(p);
! 1330: d = m0->dl; td = m0->td; c = C(m0);
! 1331: mr0 = 0;
! 1332: for ( ; m; m = NEXT(m) ) {
! 1333: NEXTNM(mr0,mr);
! 1334: C(mr) = (C(m)*c)%nd_mod;
! 1335: mr->td = m->td+td;
! 1336: ndl_add(m->dl,d,mr->dl);
! 1337: }
! 1338: NEXT(mr) = 0;
! 1339: MKND(NV(p),mr0,r);
! 1340: r->sugar = p->sugar + td;
! 1341: return r;
! 1342: }
! 1343: }
! 1344:
! 1345: ND nd_reduce(ND p1,ND p2)
! 1346: {
! 1347: int c,t,td,td2,mul;
! 1348: NM m2,prev,head,cur,new;
! 1349: unsigned int *d;
! 1350:
! 1351: if ( !p1 )
! 1352: return 0;
! 1353: else {
! 1354: mul = ((nd_mod-HC(p1))*invm(HC(p2),nd_mod))%nd_mod;
! 1355: td = HTD(p1)-HTD(p2);
! 1356: d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
! 1357: ndl_sub(HDL(p1),HDL(p2),d);
! 1358: prev = 0; head = cur = BDY(p1);
! 1359: NEWNM(new);
! 1360: for ( m2 = BDY(p2); m2; ) {
! 1361: td2 = new->td = m2->td+td;
! 1362: ndl_add(m2->dl,d,new->dl);
! 1363: if ( !cur ) {
! 1364: C(new) = (C(m2)*mul)%nd_mod;
! 1365: if ( !prev ) {
! 1366: prev = new;
! 1367: NEXT(prev) = 0;
! 1368: head = prev;
! 1369: } else {
! 1370: NEXT(prev) = new;
! 1371: NEXT(new) = 0;
! 1372: prev = new;
! 1373: }
! 1374: m2 = NEXT(m2);
! 1375: NEWNM(new);
! 1376: continue;
! 1377: }
! 1378: if ( cur->td > td2 )
! 1379: c = 1;
! 1380: else if ( cur->td < td2 )
! 1381: c = -1;
! 1382: else
! 1383: c = ndl_compare(cur->dl,new->dl);
! 1384: switch ( c ) {
! 1385: case 0:
! 1386: t = (C(cur)+C(m2)*mul)%nd_mod;
! 1387: if ( t )
! 1388: C(cur) = t;
! 1389: else if ( !prev ) {
! 1390: head = NEXT(cur);
! 1391: FREENM(cur);
! 1392: cur = head;
! 1393: } else {
! 1394: NEXT(prev) = NEXT(cur);
! 1395: FREENM(cur);
! 1396: cur = NEXT(prev);
! 1397: }
! 1398: m2 = NEXT(m2);
! 1399: break;
! 1400: case 1:
! 1401: prev = cur;
! 1402: cur = NEXT(cur);
! 1403: break;
! 1404: case -1:
! 1405: if ( !prev ) {
! 1406: /* cur = head */
! 1407: prev = new;
! 1408: C(prev) = (C(m2)*mul)%nd_mod;
! 1409: NEXT(prev) = head;
! 1410: head = prev;
! 1411: } else {
! 1412: C(new) = (C(m2)*mul)%nd_mod;
! 1413: NEXT(prev) = new;
! 1414: NEXT(new) = cur;
! 1415: prev = new;
! 1416: }
! 1417: NEWNM(new);
! 1418: m2 = NEXT(m2);
! 1419: break;
! 1420: }
! 1421: }
! 1422: FREENM(new);
! 1423: if ( head ) {
! 1424: BDY(p1) = head;
! 1425: p1->sugar = MAX(p1->sugar,p2->sugar+td);
! 1426: return p1;
! 1427: } else {
! 1428: FREEND(p1);
! 1429: return 0;
! 1430: }
! 1431:
! 1432: }
! 1433: }
! 1434:
! 1435: ND nd_sp(ND_pairs p)
! 1436: {
! 1437: NM m;
! 1438: ND p1,p2,t1,t2;
! 1439: unsigned int *lcm;
! 1440: int td;
! 1441:
! 1442: p1 = nps[p->i1];
! 1443: p2 = nps[p->i2];
! 1444: lcm = p->lcm;
! 1445: td = p->td;
! 1446: NEWNM(m);
! 1447: C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl); NEXT(m) = 0;
! 1448: t1 = nd_mul_nm(p1,m);
! 1449: C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
! 1450: t2 = nd_mul_nm(p2,m);
! 1451: FREENM(m);
! 1452: return nd_add(t1,t2);
! 1453: }
! 1454:
! 1455: ND nd_reducer(ND p1,ND p2)
! 1456: {
! 1457: NM m;
! 1458: ND r;
! 1459:
! 1460: NEWNM(m);
! 1461: C(m) = ((nd_mod-HC(p1))*invm(HC(p2),nd_mod))%nd_mod;
! 1462: m->td = HTD(p1)-HTD(p2);
! 1463: ndl_sub(HDL(p1),HDL(p2),m->dl);
! 1464: NEXT(m) = 0;
! 1465: r = nd_mul_nm(p2,m);
! 1466: FREENM(m);
! 1467: return r;
! 1468: }
! 1469:
! 1470: #if 1
! 1471: ND nd_nf(NODE b,ND g,ND *ps,int full)
! 1472: {
! 1473: ND u,p,d,red;
! 1474: NODE l;
! 1475: NM m,mrd;
! 1476: int sugar,psugar,n,h_reducible;
! 1477:
! 1478: if ( !g ) {
! 1479: return 0;
! 1480: }
! 1481: sugar = g->sugar;
! 1482: n = g->nv;
! 1483: for ( d = 0; g; ) {
! 1484: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
! 1485: p = ps[(int)BDY(l)];
! 1486: if ( HTD(g)>=HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
! 1487: h_reducible = 1;
! 1488: psugar = HTD(g)-HTD(p) + p->sugar;
! 1489: #if 0
! 1490: red = nd_reducer(g,p);
! 1491: g = nd_add(g,red);
! 1492: #else
! 1493: g = nd_reduce(g,p);
! 1494: #endif
! 1495: sugar = MAX(sugar,psugar);
! 1496: if ( !g ) {
! 1497: if ( d )
! 1498: d->sugar = sugar;
! 1499: return d;
! 1500: }
! 1501: break;
! 1502: }
! 1503: }
! 1504: if ( !h_reducible ) {
! 1505: /* head term is not reducible */
! 1506: if ( !full ) {
! 1507: if ( g )
! 1508: g->sugar = sugar;
! 1509: return g;
! 1510: } else {
! 1511: m = BDY(g);
! 1512: if ( NEXT(m) ) {
! 1513: BDY(g) = NEXT(m); NEXT(m) = 0;
! 1514: } else {
! 1515: FREEND(g); g = 0;
! 1516: }
! 1517: if ( d ) {
! 1518: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
! 1519: NEXT(mrd) = m;
! 1520: } else {
! 1521: MKND(n,m,d);
! 1522: }
! 1523: }
! 1524: }
! 1525: }
! 1526: if ( d )
! 1527: d->sugar = sugar;
! 1528: return d;
! 1529: }
! 1530: #else
! 1531:
! 1532: ND nd_remove_head(ND p)
! 1533: {
! 1534: NM m;
! 1535:
! 1536: m = BDY(p);
! 1537: if ( !NEXT(m) ) {
! 1538: FREEND(p);
! 1539: p = 0;
! 1540: } else
! 1541: BDY(p) = NEXT(m);
! 1542: FREENM(m);
! 1543: return p;
! 1544: }
! 1545:
! 1546: PGeoBucket create_pbucket()
! 1547: {
! 1548: PGeoBucket g;
! 1549:
! 1550: g = CALLOC(1,sizeof(struct oPGeoBucket));
! 1551: g->m = -1;
! 1552: return g;
! 1553: }
! 1554:
! 1555: void add_pbucket(PGeoBucket g,ND d)
! 1556: {
! 1557: int l,k,m;
! 1558:
! 1559: l = nd_length(d);
! 1560: for ( k = 0, m = 1; l > m; k++, m <<= 2 );
! 1561: /* 4^(k-1) < l <= 4^k */
! 1562: d = nd_add(g->body[k],d);
! 1563: for ( ; d && nd_length(d) > 1<<(2*k); k++ ) {
! 1564: g->body[k] = 0;
! 1565: d = nd_add(g->body[k+1],d);
! 1566: }
! 1567: g->body[k] = d;
! 1568: g->m = MAX(g->m,k);
! 1569: }
! 1570:
! 1571: int head_pbucket(PGeoBucket g)
! 1572: {
! 1573: int j,i,c,k,nv,sum;
! 1574: unsigned int *di,*dj;
! 1575: ND gi,gj;
! 1576:
! 1577: k = g->m;
! 1578: while ( 1 ) {
! 1579: j = -1;
! 1580: for ( i = 0; i <= k; i++ ) {
! 1581: if ( !(gi = g->body[i]) )
! 1582: continue;
! 1583: if ( j < 0 ) {
! 1584: j = i;
! 1585: gj = g->body[j];
! 1586: dj = HDL(gj);
! 1587: sum = HC(gj);
! 1588: } else {
! 1589: di = HDL(gi);
! 1590: nv = NV(gi);
! 1591: if ( HTD(gi) > HTD(gj) )
! 1592: c = 1;
! 1593: else if ( HTD(gi) < HTD(gj) )
! 1594: c = -1;
! 1595: else
! 1596: c = ndl_compare(di,dj);
! 1597: if ( c > 0 ) {
! 1598: if ( sum )
! 1599: HC(gj) = sum;
! 1600: else
! 1601: g->body[j] = nd_remove_head(gj);
! 1602: j = i;
! 1603: gj = g->body[j];
! 1604: dj = HDL(gj);
! 1605: sum = HC(gj);
! 1606: } else if ( c == 0 ) {
! 1607: sum = sum+HC(gi)-nd_mod;
! 1608: if ( sum < 0 )
! 1609: sum += nd_mod;
! 1610: g->body[i] = nd_remove_head(gi);
! 1611: }
! 1612: }
! 1613: }
! 1614: if ( j < 0 )
! 1615: return -1;
! 1616: else if ( sum ) {
! 1617: HC(gj) = sum;
! 1618: return j;
! 1619: } else
! 1620: g->body[j] = nd_remove_head(gj);
! 1621: }
! 1622: }
! 1623:
! 1624: ND normalize_pbucket(PGeoBucket g)
! 1625: {
! 1626: int i;
! 1627: ND r,t;
! 1628:
! 1629: r = 0;
! 1630: for ( i = 0; i <= g->m; i++ )
! 1631: r = nd_add(r,g->body[i]);
! 1632: return r;
! 1633: }
! 1634:
! 1635: ND nd_nf(NODE b,ND g,ND *ps,int full)
! 1636: {
! 1637: ND u,p,d,red;
! 1638: NODE l;
! 1639: NM m,mrd;
! 1640: int sugar,psugar,n,h_reducible,h;
! 1641: PGeoBucket bucket;
! 1642:
! 1643: if ( !g ) {
! 1644: return 0;
! 1645: }
! 1646: sugar = g->sugar;
! 1647: n = g->nv;
! 1648: bucket = create_pbucket();
! 1649: add_pbucket(bucket,g);
! 1650: d = 0;
! 1651: while ( 1 ) {
! 1652: h = head_pbucket(bucket);
! 1653: if ( h < 0 ) {
! 1654: if ( d )
! 1655: d->sugar = sugar;
! 1656: return d;
! 1657: }
! 1658: g = bucket->body[h];
! 1659: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
! 1660: p = ps[(int)BDY(l)];
! 1661: if ( ndl_reducible(HDL(g),HDL(p)) ) {
! 1662: h_reducible = 1;
! 1663: psugar = HTD(g)-HTD(p) + p->sugar;
! 1664: red = nd_reducer(g,p);
! 1665: bucket->body[h] = nd_remove_head(g);
! 1666: red = nd_remove_head(red);
! 1667: add_pbucket(bucket,red);
! 1668: sugar = MAX(sugar,psugar);
! 1669: break;
! 1670: }
! 1671: }
! 1672: if ( !h_reducible ) {
! 1673: /* head term is not reducible */
! 1674: if ( !full ) {
! 1675: g = normalize_pbucket(bucket);
! 1676: if ( g )
! 1677: g->sugar = sugar;
! 1678: return g;
! 1679: } else {
! 1680: m = BDY(g);
! 1681: if ( NEXT(m) ) {
! 1682: BDY(g) = NEXT(m); NEXT(m) = 0;
! 1683: } else {
! 1684: FREEND(g); g = 0;
! 1685: }
! 1686: bucket->body[h] = g;
! 1687: NEXT(m) = 0;
! 1688: if ( d ) {
! 1689: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
! 1690: NEXT(mrd) = m;
! 1691: } else {
! 1692: MKND(n,m,d);
! 1693: }
! 1694: }
! 1695: }
! 1696: }
! 1697: }
! 1698: #endif
! 1699:
! 1700: NODE nd_gb(NODE f)
! 1701: {
! 1702: int i,nh;
! 1703: NODE r,g,gall;
! 1704: ND_pairs d;
! 1705: ND_pairs l;
! 1706: ND h,nf;
! 1707:
! 1708: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
! 1709: i = (int)BDY(r);
! 1710: d = update_pairs(d,g,i);
! 1711: g = update_base(g,i);
! 1712: gall = append_one(gall,i);
! 1713: }
! 1714: while ( d ) {
! 1715: #if 0
! 1716: ndp_print(d);
! 1717: #endif
! 1718: l = nd_minp(d,&d);
! 1719: h = nd_sp(l);
! 1720: nf = nd_nf(gall,h,nps,!Top);
! 1721: if ( nf ) {
! 1722: printf("+"); fflush(stdout);
! 1723: #if 0
! 1724: ndl_print(HDL(nf)); fflush(stdout);
! 1725: #endif
! 1726: nh = nd_newps(nf);
! 1727: d = update_pairs(d,g,nh);
! 1728: g = update_base(g,nh);
! 1729: gall = append_one(gall,nh);
! 1730: } else {
! 1731: printf("."); fflush(stdout);
! 1732: }
! 1733: }
! 1734: return g;
! 1735: }
! 1736:
! 1737: ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t)
! 1738: {
! 1739: ND_pairs d1,nd,cur,head,prev;
! 1740:
! 1741: if ( !g ) return d;
! 1742: d = crit_B(d,t);
! 1743: d1 = nd_newpairs(g,t);
! 1744: d1 = crit_M(d1);
! 1745: d1 = crit_F(d1);
! 1746: prev = 0; cur = head = d1;
! 1747: while ( cur ) {
! 1748: if ( crit_2( cur->i1,cur->i2 ) ) {
! 1749: if ( !prev ) {
! 1750: head = cur = NEXT(cur);
! 1751: } else {
! 1752: cur = NEXT(prev) = NEXT(cur);
! 1753: }
! 1754: } else {
! 1755: prev = cur;
! 1756: cur = NEXT(cur);
! 1757: }
! 1758: }
! 1759: if ( !d )
! 1760: return head;
! 1761: else {
! 1762: nd = d;
! 1763: while ( NEXT(nd) )
! 1764: nd = NEXT(nd);
! 1765: NEXT(nd) = head;
! 1766: return d;
! 1767: }
! 1768: }
! 1769:
! 1770: ND_pairs nd_newpairs( NODE g, int t )
! 1771: {
! 1772: NODE h;
! 1773: unsigned int *dl;
! 1774: int td,ts,s;
! 1775: ND_pairs r,r0;
! 1776:
! 1777: dl = HDL(nps[t]);
! 1778: td = HTD(nps[t]);
! 1779: ts = nps[t]->sugar - td;
! 1780: for ( r0 = 0, h = g; h; h = NEXT(h) ) {
! 1781: NEXTND_pairs(r0,r);
! 1782: r->i1 = (int)BDY(h);
! 1783: r->i2 = t;
! 1784: ndl_lcm(HDL(nps[r->i1]),dl,r->lcm);
! 1785: r->td = ndl_td(r->lcm);
! 1786: s = nps[r->i1]->sugar-HTD(nps[r->i1]);
! 1787: r->sugar = MAX(s,ts) + r->td;
! 1788: }
! 1789: NEXT(r) = 0;
! 1790: return r0;
! 1791: }
! 1792:
! 1793: ND_pairs crit_B( ND_pairs d, int s )
! 1794: {
! 1795: ND_pairs cur,head,prev;
! 1796: unsigned int *t,*tl,*lcm;
! 1797: int td,tdl;
! 1798:
! 1799: if ( !d ) return 0;
! 1800: t = HDL(nps[s]);
! 1801: prev = 0;
! 1802: head = cur = d;
! 1803: lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
! 1804: while ( cur ) {
! 1805: tl = cur->lcm;
! 1806: if ( ndl_reducible(tl,t)
! 1807: && (ndl_lcm(HDL(nps[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
! 1808: && (ndl_lcm(HDL(nps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
! 1809: if ( !prev ) {
! 1810: head = cur = NEXT(cur);
! 1811: } else {
! 1812: cur = NEXT(prev) = NEXT(cur);
! 1813: }
! 1814: } else {
! 1815: prev = cur;
! 1816: cur = NEXT(cur);
! 1817: }
! 1818: }
! 1819: return head;
! 1820: }
! 1821:
! 1822: ND_pairs crit_M( ND_pairs d1 )
! 1823: {
! 1824: ND_pairs e,d2,d3,dd,p;
! 1825: unsigned int *id,*jd;
! 1826: int itd,jtd;
! 1827:
! 1828: for ( dd = 0, e = d1; e; e = d3 ) {
! 1829: if ( !(d2 = NEXT(e)) ) {
! 1830: NEXT(e) = dd;
! 1831: return e;
! 1832: }
! 1833: id = e->lcm;
! 1834: itd = e->td;
! 1835: for ( d3 = 0; d2; d2 = p ) {
! 1836: p = NEXT(d2),
! 1837: jd = d2->lcm;
! 1838: jtd = d2->td;
! 1839: if ( jtd == itd )
! 1840: if ( id == jd );
! 1841: else if ( ndl_reducible(jd,id) ) continue;
! 1842: else if ( ndl_reducible(id,jd) ) goto delit;
! 1843: else ;
! 1844: else if ( jtd > itd )
! 1845: if ( ndl_reducible(jd,id) ) continue;
! 1846: else ;
! 1847: else if ( ndl_reducible(id,jd ) ) goto delit;
! 1848: NEXT(d2) = d3;
! 1849: d3 = d2;
! 1850: }
! 1851: NEXT(e) = dd;
! 1852: dd = e;
! 1853: continue;
! 1854: /**/
! 1855: delit: NEXT(d2) = d3;
! 1856: d3 = d2;
! 1857: for ( ; p; p = d2 ) {
! 1858: d2 = NEXT(p);
! 1859: NEXT(p) = d3;
! 1860: d3 = p;
! 1861: }
! 1862: }
! 1863: return dd;
! 1864: }
! 1865:
! 1866: ND_pairs crit_F( ND_pairs d1 )
! 1867: {
! 1868: ND_pairs rest, head;
! 1869: ND_pairs last, p, r, w;
! 1870: int s;
! 1871:
! 1872: for ( head = last = 0, p = d1; NEXT(p); ) {
! 1873: r = w = equivalent_pairs(p,&rest);
! 1874: s = r->sugar;
! 1875: while ( w = NEXT(w) )
! 1876: if ( crit_2(w->i1,w->i2) ) {
! 1877: r = w;
! 1878: break;
! 1879: } else if ( w->sugar < s ) {
! 1880: r = w;
! 1881: s = r->sugar;
! 1882: }
! 1883: if ( last ) NEXT(last) = r;
! 1884: else head = r;
! 1885: NEXT(last = r) = 0;
! 1886: p = rest;
! 1887: if ( !p ) return head;
! 1888: }
! 1889: if ( !last ) return p;
! 1890: NEXT(last) = p;
! 1891: return head;
! 1892: }
! 1893:
! 1894: int crit_2( int dp1, int dp2 )
! 1895: {
! 1896: return ndl_disjoint(HDL(nps[dp1]),HDL(nps[dp2]));
! 1897: }
! 1898:
! 1899: static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest )
! 1900: {
! 1901: ND_pairs w,p,r,s;
! 1902: unsigned int *d;
! 1903: int td;
! 1904:
! 1905: w = d1;
! 1906: d = w->lcm;
! 1907: td = w->td;
! 1908: s = NEXT(w);
! 1909: NEXT(w) = 0;
! 1910: for ( r = 0; s; s = p ) {
! 1911: p = NEXT(s);
! 1912: if ( td == s->td && ndl_equal(d,s->lcm) ) {
! 1913: NEXT(s) = w;
! 1914: w = s;
! 1915: } else {
! 1916: NEXT(s) = r;
! 1917: r = s;
! 1918: }
! 1919: }
! 1920: *prest = r;
! 1921: return w;
! 1922: }
! 1923:
! 1924: NODE update_base(NODE nd,int ndp)
! 1925: {
! 1926: unsigned int *dl, *dln;
! 1927: NODE last, p, head;
! 1928: int td,tdn;
! 1929:
! 1930: dl = HDL(nps[ndp]);
! 1931: td = HTD(nps[ndp]);
! 1932: for ( head = last = 0, p = nd; p; ) {
! 1933: dln = HDL(nps[(int)BDY(p)]);
! 1934: tdn = HTD(nps[(int)BDY(p)]);
! 1935: if ( tdn >= td && ndl_reducible( dln, dl ) ) {
! 1936: p = NEXT(p);
! 1937: if ( last ) NEXT(last) = p;
! 1938: } else {
! 1939: if ( !last ) head = p;
! 1940: p = NEXT(last = p);
! 1941: }
! 1942: }
! 1943: head = append_one(head,ndp);
! 1944: return head;
! 1945: }
! 1946:
! 1947: ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
! 1948: {
! 1949: ND_pairs m,ml,p,l;
! 1950: unsigned int *lcm;
! 1951: int s,td,len,tlen,c;
! 1952:
! 1953: if ( !(p = NEXT(m = d)) ) {
! 1954: *prest = p;
! 1955: NEXT(m) = 0;
! 1956: return m;
! 1957: }
! 1958: lcm = m->lcm;
! 1959: s = m->sugar;
! 1960: td = m->td;
! 1961: len = nd_length(nps[m->i1])+nd_length(nps[m->i2]);
! 1962: for ( ml = 0, l = m; p; p = NEXT(l = p) ) {
! 1963: if (p->sugar < s)
! 1964: goto find;
! 1965: else if ( p->sugar == s ) {
! 1966: if ( p->td < td )
! 1967: goto find;
! 1968: else if ( p->td == td ) {
! 1969: c = ndl_compare(p->lcm,lcm);
! 1970: if ( c < 0 )
! 1971: goto find;
! 1972: else if ( c == 0 ) {
! 1973: tlen = nd_length(nps[p->i1])+nd_length(nps[p->i2]);
! 1974: if ( tlen < len )
! 1975: goto find;
! 1976: }
! 1977: }
! 1978: }
! 1979: continue;
! 1980: find:
! 1981: ml = l;
! 1982: m = p;
! 1983: lcm = m->lcm;
! 1984: s = m->sugar;
! 1985: td = m->td;
! 1986: len = tlen;
! 1987: }
! 1988: if ( !ml ) *prest = NEXT(m);
! 1989: else {
! 1990: NEXT(ml) = NEXT(m);
! 1991: *prest = d;
! 1992: }
! 1993: NEXT(m) = 0;
! 1994: return m;
! 1995: }
! 1996:
! 1997: int nd_newps(ND a)
! 1998: {
! 1999: if ( nd_psn == nd_pslen ) {
! 2000: nd_pslen *= 2;
! 2001: nps = (ND *)REALLOC((char *)nps,nd_pslen*sizeof(ND));
! 2002: }
! 2003: nd_monic(a);
! 2004: nps[nd_psn] = a;
! 2005: return nd_psn++;
! 2006: }
! 2007:
! 2008: NODE NODE_sortb(NODE f,int);
! 2009: ND dptond(DP);
! 2010: DP ndtodp(ND);
! 2011:
! 2012: NODE nd_setup(NODE f)
! 2013: {
! 2014: int i;
! 2015: NODE s,s0,f0;
! 2016:
! 2017: #if 0
! 2018: f0 = f = NODE_sortb(f,1);
! 2019: #endif
! 2020: nd_psn = length(f); nd_pslen = 2*nd_psn;
! 2021: nps = (ND *)MALLOC(nd_pslen*sizeof(ND));
! 2022: nd_bpe = 4;
! 2023: nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
! 2024: nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
! 2025: nd_free_private_storage();
! 2026: for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
! 2027: nps[i] = dptond((DP)BDY(f));
! 2028: nd_monic(nps[i]);
! 2029: }
! 2030: for ( s0 = 0, i = 0; i < nd_psn; i++ ) {
! 2031: NEXTNODE(s0,s); BDY(s) = (pointer)i;
! 2032: }
! 2033: if ( s0 ) NEXT(s) = 0;
! 2034: return s0;
! 2035: }
! 2036:
! 2037: void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
! 2038: {
! 2039: struct order_spec ord1;
! 2040: VL fv,vv,vc;
! 2041: NODE fd,fd0,r,r0,t,x,s,xx;
! 2042: DP a,b,c;
! 2043:
! 2044: get_vars((Obj)f,&fv); pltovl(v,&vv);
! 2045: nd_nvar = length(vv);
! 2046: if ( ord->id )
! 2047: error("nd_gr : unsupported order");
! 2048: switch ( ord->ord.simple ) {
! 2049: case 0:
! 2050: is_rlex = 1;
! 2051: break;
! 2052: case 1:
! 2053: is_rlex = 0;
! 2054: break;
! 2055: default:
! 2056: error("nd_gr : unsupported order");
! 2057: }
! 2058: initd(ord);
! 2059: nd_mod = m;
! 2060: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
! 2061: ptod(CO,vv,(P)BDY(t),&b);
! 2062: _dp_mod(b,m,0,&c);
! 2063: if ( c ) {
! 2064: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
! 2065: }
! 2066: }
! 2067: if ( fd0 ) NEXT(fd) = 0;
! 2068: s = nd_setup(fd0);
! 2069: x = nd_gb(s);
! 2070: #if 0
! 2071: x = nd_reduceall(x,m);
! 2072: #endif
! 2073: for ( r0 = 0; x; x = NEXT(x) ) {
! 2074: NEXTNODE(r0,r);
! 2075: a = ndtodp(nps[(int)BDY(x)]);
! 2076: _dtop_mod(CO,vv,a,(P *)&BDY(r));
! 2077: }
! 2078: if ( r0 ) NEXT(r) = 0;
! 2079: MKLIST(*rp,r0);
! 2080: }
! 2081:
! 2082: void dltondl(int n,DL dl,unsigned int *r)
! 2083: {
! 2084: unsigned int *d;
! 2085: int i;
! 2086:
! 2087: d = dl->d;
! 2088: if ( is_rlex )
! 2089: for ( i = 0; i < n; i++ )
! 2090: r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe));
! 2091: else
! 2092: for ( i = 0; i < n; i++ )
! 2093: r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe);
! 2094: }
! 2095:
! 2096: DL ndltodl(int n,int td,unsigned int *ndl)
! 2097: {
! 2098: DL dl;
! 2099: int *d;
! 2100: int i;
! 2101:
! 2102: NEWDL(dl,n);
! 2103: dl->td = td;
! 2104: d = dl->d;
! 2105: if ( is_rlex )
! 2106: for ( i = 0; i < n; i++ )
! 2107: d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
! 2108: &((1<<nd_bpe)-1);
! 2109: else
! 2110: for ( i = 0; i < n; i++ )
! 2111: d[i] = (ndl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
! 2112: &((1<<nd_bpe)-1);
! 2113: return dl;
! 2114: }
! 2115:
! 2116: ND dptond(DP p)
! 2117: {
! 2118: ND d;
! 2119: NM m0,m;
! 2120: MP t;
! 2121: int n;
! 2122:
! 2123: if ( !p )
! 2124: return 0;
! 2125: n = NV(p);
! 2126: m0 = 0;
! 2127: for ( t = BDY(p); t; t = NEXT(t) ) {
! 2128: NEXTNM(m0,m);
! 2129: m->c = ITOS(t->c);
! 2130: m->td = t->dl->td;
! 2131: dltondl(n,t->dl,m->dl);
! 2132: }
! 2133: NEXT(m) = 0;
! 2134: MKND(n,m0,d);
! 2135: d->nv = n;
! 2136: d->sugar = p->sugar;
! 2137: return d;
! 2138: }
! 2139:
! 2140: DP ndtodp(ND p)
! 2141: {
! 2142: DP d;
! 2143: MP m0,m;
! 2144: NM t;
! 2145: int n;
! 2146:
! 2147: if ( !p )
! 2148: return 0;
! 2149: n = NV(p);
! 2150: m0 = 0;
! 2151: for ( t = BDY(p); t; t = NEXT(t) ) {
! 2152: NEXTMP(m0,m);
! 2153: m->c = STOI(t->c);
! 2154: m->dl = ndltodl(n,t->td,t->dl);
! 2155: }
! 2156: NEXT(m) = 0;
! 2157: MKDP(n,m0,d);
! 2158: d->sugar = p->sugar;
! 2159: return d;
! 2160: }
! 2161:
! 2162: void ndl_print(unsigned int *dl)
! 2163: {
! 2164: int n;
! 2165: int i;
! 2166:
! 2167: n = nd_nvar;
! 2168: printf("<<");
! 2169: if ( is_rlex )
! 2170: for ( i = 0; i < n; i++ )
! 2171: printf(i==n-1?"%d":"%d,",
! 2172: (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe))
! 2173: &((1<<nd_bpe)-1));
! 2174: else
! 2175: for ( i = 0; i < n; i++ )
! 2176: printf(i==n-1?"%d":"%d,",
! 2177: (dl[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))
! 2178: &((1<<nd_bpe)-1));
! 2179: printf(">>");
! 2180: }
! 2181:
! 2182: void nd_print(ND p)
! 2183: {
! 2184: NM m;
! 2185:
! 2186: if ( !p )
! 2187: printf("0\n");
! 2188: else {
! 2189: for ( m = BDY(p); m; m = NEXT(m) ) {
! 2190: printf("+%d*",m->c);
! 2191: ndl_print(m->dl);
! 2192: }
! 2193: printf("\n");
! 2194: }
! 2195: }
! 2196:
! 2197: void ndp_print(ND_pairs d)
! 2198: {
! 2199: ND_pairs t;
! 2200:
! 2201: for ( t = d; t; t = NEXT(t) ) {
! 2202: printf("%d,%d ",t->i1,t->i2);
! 2203: }
! 2204: printf("\n");
! 2205: }
! 2206:
! 2207: void nd_monic(ND p)
! 2208: {
! 2209: int mul;
! 2210: NM m;
! 2211:
! 2212: if ( !p )
! 2213: return;
! 2214: mul = invm(HC(p),nd_mod);
! 2215: for ( m = BDY(p); m; m = NEXT(m) )
! 2216: C(m) = (C(m)*mul)%nd_mod;
1.1 noro 2217: }