![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dist.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 dist.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.5 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.4 2000/04/25 04:07:59 noro Exp $ */
1.1 noro 2: #include "ca.h"
3:
4: #define NV(p) ((p)->nv)
5: #define C(p) ((p)->c)
6:
7: #define ORD_REVGRADLEX 0
8: #define ORD_GRADLEX 1
9: #define ORD_LEX 2
10: #define ORD_BREVGRADLEX 3
11: #define ORD_BGRADLEX 4
12: #define ORD_BLEX 5
13: #define ORD_BREVREV 6
14: #define ORD_BGRADREV 7
15: #define ORD_BLEXREV 8
16: #define ORD_ELIM 9
17:
18: int (*cmpdl)()=cmpdl_revgradlex;
19: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
20:
1.2 noro 21: void comm_muld(VL,DP,DP,DP *);
22: void weyl_muld(VL,DP,DP,DP *);
23: void weyl_muldm(VL,DP,MP,DP *);
24: void weyl_mulmm(VL,MP,MP,int,DP *);
25: void mkwc(int,int,Q *);
26:
27: int do_weyl;
28:
1.1 noro 29: int dp_nelim,dp_fcoeffs;
30: struct order_spec dp_current_spec;
31: int *dp_dl_work;
32:
33: int has_fcoef(DP);
34: int has_fcoef_p(P);
35:
36: int has_fcoef(f)
37: DP f;
38: {
39: MP t;
40:
41: if ( !f )
42: return 0;
43: for ( t = BDY(f); t; t = NEXT(t) )
44: if ( has_fcoef_p(t->c) )
45: break;
46: return t ? 1 : 0;
47: }
48:
49: int has_fcoef_p(f)
50: P f;
51: {
52: DCP dc;
53:
54: if ( !f )
55: return 0;
56: else if ( NUM(f) )
57: return (NID((Num)f) == N_LM || NID((Num)f) == N_GF2N) ? 1 : 0;
58: else {
59: for ( dc = DC(f); dc; dc = NEXT(dc) )
60: if ( has_fcoef_p(COEF(dc)) )
61: return 1;
62: return 0;
63: }
64: }
65:
66: void initd(spec)
67: struct order_spec *spec;
68: {
69: switch ( spec->id ) {
70: case 2:
71: cmpdl = cmpdl_matrix;
72: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
73: break;
74: case 1:
75: cmpdl = cmpdl_order_pair;
76: break;
77: default:
78: switch ( spec->ord.simple ) {
79: case ORD_REVGRADLEX:
80: cmpdl = cmpdl_revgradlex; break;
81: case ORD_GRADLEX:
82: cmpdl = cmpdl_gradlex; break;
83: case ORD_BREVGRADLEX:
84: cmpdl = cmpdl_brevgradlex; break;
85: case ORD_BGRADLEX:
86: cmpdl = cmpdl_bgradlex; break;
87: case ORD_BLEX:
88: cmpdl = cmpdl_blex; break;
89: case ORD_BREVREV:
90: cmpdl = cmpdl_brevrev; break;
91: case ORD_BGRADREV:
92: cmpdl = cmpdl_bgradrev; break;
93: case ORD_BLEXREV:
94: cmpdl = cmpdl_blexrev; break;
95: case ORD_ELIM:
96: cmpdl = cmpdl_elim; break;
97: case ORD_LEX: default:
98: cmpdl = cmpdl_lex; break;
99: }
100: break;
101: }
102: dp_current_spec = *spec;
103: }
104:
105: void ptod(vl,dvl,p,pr)
106: VL vl,dvl;
107: P p;
108: DP *pr;
109: {
110: int isconst = 0;
111: int n,i;
112: VL tvl;
113: V v;
114: DL d;
115: MP m;
116: DCP dc;
117: DP r,s,t,u;
118: P x,c;
119:
120: if ( !p )
121: *pr = 0;
122: else {
123: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
124: if ( NUM(p) ) {
125: NEWDL(d,n);
126: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
127: } else {
128: for ( i = 0, tvl = dvl, v = VR(p);
129: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
130: if ( !tvl ) {
131: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
132: ptod(vl,dvl,COEF(dc),&t); pwrp(vl,x,DEG(dc),&c);
133: muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
134: }
135: *pr = s;
136: } else {
137: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
138: ptod(vl,dvl,COEF(dc),&t);
139: NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
140: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
1.2 noro 141: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
1.1 noro 142: }
143: *pr = s;
144: }
145: }
146: }
147: if ( !dp_fcoeffs && has_fcoef(*pr) )
148: dp_fcoeffs = 1;
149: }
150:
151: void dtop(vl,dvl,p,pr)
152: VL vl,dvl;
153: DP p;
154: P *pr;
155: {
156: int n,i;
157: DL d;
158: MP m;
159: P r,s,t,u,w;
160: Q q;
161: VL tvl;
162:
163: if ( !p )
164: *pr = 0;
165: else {
166: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
167: t = C(m);
168: if ( NUM(t) && NID((Num)t) == N_M ) {
169: mptop(t,&u); t = u;
170: }
171: for ( i = 0, d = m->dl, tvl = dvl;
172: i < n; tvl = NEXT(tvl), i++ ) {
173: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
174: mulp(vl,t,u,&w); t = w;
175: }
176: addp(vl,s,t,&u); s = u;
177: }
178: *pr = s;
179: }
180: }
181:
182: void nodetod(node,dp)
183: NODE node;
184: DP *dp;
185: {
186: NODE t;
187: int len,i,td;
188: Q e;
189: DL d;
190: MP m;
191: DP u;
192:
193: for ( t = node, len = 0; t; t = NEXT(t), len++ );
194: NEWDL(d,len);
195: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
196: e = (Q)BDY(t);
197: if ( !e )
198: d->d[i] = 0;
199: else if ( !NUM(e) || !RATN(e) || !INT(e) )
200: error("nodetod : invalid input");
201: else {
202: d->d[i] = QTOS((Q)e); td += d->d[i];
203: }
204: }
205: d->td = td;
206: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0;
207: MKDP(len,m,u); u->sugar = td; *dp = u;
208: }
209:
210: int sugard(m)
211: MP m;
212: {
213: int s;
214:
215: for ( s = 0; m; m = NEXT(m) )
216: s = MAX(s,m->dl->td);
217: return s;
218: }
219:
220: void addd(vl,p1,p2,pr)
221: VL vl;
222: DP p1,p2,*pr;
223: {
224: int n;
225: MP m1,m2,mr,mr0;
226: P t;
227:
228: if ( !p1 )
229: *pr = p2;
230: else if ( !p2 )
231: *pr = p1;
232: else {
233: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
234: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
235: case 0:
236: addp(vl,C(m1),C(m2),&t);
237: if ( t ) {
238: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
239: }
240: m1 = NEXT(m1); m2 = NEXT(m2); break;
241: case 1:
242: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
243: m1 = NEXT(m1); break;
244: case -1:
245: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
246: m2 = NEXT(m2); break;
247: }
248: if ( !mr0 )
249: if ( m1 )
250: mr0 = m1;
251: else if ( m2 )
252: mr0 = m2;
253: else {
254: *pr = 0;
255: return;
256: }
257: else if ( m1 )
258: NEXT(mr) = m1;
259: else if ( m2 )
260: NEXT(mr) = m2;
261: else
262: NEXT(mr) = 0;
263: MKDP(NV(p1),mr0,*pr);
264: if ( *pr )
265: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
266: }
267: }
268:
269: /* for F4 symbolic reduction */
270:
271: void symb_addd(p1,p2,pr)
272: DP p1,p2,*pr;
273: {
274: int n;
275: MP m1,m2,mr,mr0;
276: P t;
277:
278: if ( !p1 )
279: *pr = p2;
280: else if ( !p2 )
281: *pr = p1;
282: else {
283: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
284: NEXTMP(mr0,mr); C(mr) = (P)ONE;
285: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
286: case 0:
287: mr->dl = m1->dl;
288: m1 = NEXT(m1); m2 = NEXT(m2); break;
289: case 1:
290: mr->dl = m1->dl;
291: m1 = NEXT(m1); break;
292: case -1:
293: mr->dl = m2->dl;
294: m2 = NEXT(m2); break;
295: }
296: }
297: if ( !mr0 )
298: if ( m1 )
299: mr0 = m1;
300: else if ( m2 )
301: mr0 = m2;
302: else {
303: *pr = 0;
304: return;
305: }
306: else if ( m1 )
307: NEXT(mr) = m1;
308: else if ( m2 )
309: NEXT(mr) = m2;
310: else
311: NEXT(mr) = 0;
312: MKDP(NV(p1),mr0,*pr);
313: if ( *pr )
314: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.3 noro 315: }
316: }
317:
318: /*
319: * destructive merge of two list
320: *
321: * p1, p2 : list of DL
322: * return : a merged list
323: */
324:
325: NODE symb_merge(m1,m2,n)
326: NODE m1,m2;
327: int n;
328: {
329: NODE top,prev,cur,m,t;
330:
331: if ( !m1 )
332: return m2;
333: else if ( !m2 )
334: return m1;
335: else {
336: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
337: case 0:
338: top = m1; m = NEXT(m2);
339: break;
340: case 1:
341: top = m1; m = m2;
342: break;
343: case -1:
344: top = m2; m = m1;
345: break;
346: }
347: prev = top; cur = NEXT(top);
348: /* BDY(prev) > BDY(m) always holds */
349: while ( cur && m ) {
350: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
351: case 0:
352: m = NEXT(m);
353: prev = cur; cur = NEXT(cur);
354: break;
355: case 1:
356: t = NEXT(cur); NEXT(cur) = m; m = t;
357: prev = cur; cur = NEXT(cur);
358: break;
359: case -1:
360: NEXT(prev) = m; m = cur;
361: prev = NEXT(prev); cur = NEXT(prev);
362: break;
363: }
364: }
365: if ( !cur )
366: NEXT(prev) = m;
367: return top;
1.1 noro 368: }
369: }
370:
371: void subd(vl,p1,p2,pr)
372: VL vl;
373: DP p1,p2,*pr;
374: {
375: DP t;
376:
377: if ( !p2 )
378: *pr = p1;
379: else {
380: chsgnd(p2,&t); addd(vl,p1,t,pr);
381: }
382: }
383:
384: void chsgnd(p,pr)
385: DP p,*pr;
386: {
387: MP m,mr,mr0;
388:
389: if ( !p )
390: *pr = 0;
391: else {
392: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
393: NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
394: }
395: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
396: if ( *pr )
397: (*pr)->sugar = p->sugar;
398: }
399: }
400:
401: void muld(vl,p1,p2,pr)
402: VL vl;
403: DP p1,p2,*pr;
404: {
1.2 noro 405: if ( ! do_weyl )
406: comm_muld(vl,p1,p2,pr);
407: else
408: weyl_muld(vl,p1,p2,pr);
409: }
410:
411: void comm_muld(vl,p1,p2,pr)
412: VL vl;
413: DP p1,p2,*pr;
414: {
1.1 noro 415: MP m;
416: DP s,t,u;
1.5 ! noro 417: int i,l,l1;
! 418: static MP *w;
! 419: static int wlen;
1.1 noro 420:
421: if ( !p1 || !p2 )
422: *pr = 0;
423: else if ( OID(p1) <= O_P )
424: muldc(vl,p2,(P)p1,pr);
425: else if ( OID(p2) <= O_P )
426: muldc(vl,p1,(P)p2,pr);
427: else {
1.5 ! noro 428: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
1.4 noro 429: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5 ! noro 430: if ( l1 < l ) {
! 431: t = p1; p1 = p2; p2 = t;
! 432: l = l1;
! 433: }
! 434: if ( l > wlen ) {
! 435: if ( w ) GC_free(w);
! 436: w = (MP *)MALLOC(l*sizeof(MP));
! 437: wlen = l;
! 438: }
1.4 noro 439: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
440: w[i] = m;
441: for ( s = 0, i = l-1; i >= 0; i-- ) {
442: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1 noro 443: }
1.5 ! noro 444: bzero(w,l*sizeof(MP));
1.1 noro 445: *pr = s;
446: }
447: }
448:
449: void muldm(vl,p,m0,pr)
450: VL vl;
451: DP p;
452: MP m0;
453: DP *pr;
454: {
455: MP m,mr,mr0;
456: P c;
457: DL d;
458: int n;
459:
460: if ( !p )
461: *pr = 0;
462: else {
463: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
464: m; m = NEXT(m) ) {
465: NEXTMP(mr0,mr);
466: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
467: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
468: else
469: mulp(vl,C(m),c,&C(mr));
470: adddl(n,m->dl,d,&mr->dl);
471: }
472: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
473: if ( *pr )
474: (*pr)->sugar = p->sugar + m0->dl->td;
1.2 noro 475: }
476: }
477:
478: void weyl_muld(vl,p1,p2,pr)
479: VL vl;
480: DP p1,p2,*pr;
481: {
482: MP m;
483: DP s,t,u;
1.4 noro 484: int i,l;
1.5 ! noro 485: static MP *w;
! 486: static int wlen;
1.2 noro 487:
488: if ( !p1 || !p2 )
489: *pr = 0;
490: else if ( OID(p1) <= O_P )
491: muldc(vl,p2,(P)p1,pr);
492: else if ( OID(p2) <= O_P )
493: muldc(vl,p1,(P)p2,pr);
494: else {
1.4 noro 495: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5 ! noro 496: if ( l > wlen ) {
! 497: if ( w ) GC_free(w);
! 498: w = (MP *)MALLOC(l*sizeof(MP));
! 499: wlen = l;
! 500: }
1.4 noro 501: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
502: w[i] = m;
503: for ( s = 0, i = l-1; i >= 0; i-- ) {
504: weyl_muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.2 noro 505: }
1.5 ! noro 506: bzero(w,l*sizeof(MP));
1.2 noro 507: *pr = s;
508: }
509: }
510:
511: void weyl_muldm(vl,p,m0,pr)
512: VL vl;
513: DP p;
514: MP m0;
515: DP *pr;
516: {
517: DP r,t,t1;
518: MP m;
1.4 noro 519: int n,l,i;
1.5 ! noro 520: static MP *w;
! 521: static int wlen;
1.2 noro 522:
523: if ( !p )
524: *pr = 0;
525: else {
1.4 noro 526: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5 ! noro 527: if ( l > wlen ) {
! 528: if ( w ) GC_free(w);
! 529: w = (MP *)MALLOC(l*sizeof(MP));
! 530: wlen = l;
! 531: }
1.4 noro 532: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
533: w[i] = m;
534: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
535: weyl_mulmm(vl,w[i],m0,n,&t);
1.2 noro 536: addd(vl,r,t,&t1); r = t1;
537: }
1.5 ! noro 538: bzero(w,l*sizeof(MP));
1.2 noro 539: if ( r )
540: r->sugar = p->sugar + m0->dl->td;
541: *pr = r;
542: }
543: }
544:
545: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
546:
547: void weyl_mulmm(vl,m0,m1,n,pr)
548: VL vl;
549: MP m0,m1;
550: int n;
551: DP *pr;
552: {
553: MP m,mr,mr0;
554: DP r,t,t1;
555: P c,c0,c1,cc;
556: DL d,d0,d1;
1.5 ! noro 557: int i,j,a,b,k,l,n2,s,min,h;
! 558: static Q *tab;
! 559: static int tablen;
1.2 noro 560:
561: if ( !m0 || !m1 )
562: *pr = 0;
563: else {
564: c0 = C(m0); c1 = C(m1);
565: mulp(vl,c0,c1,&c);
566: d0 = m0->dl; d1 = m1->dl;
567: n2 = n>>1;
1.5 ! noro 568: if ( n & 1 ) {
! 569: /* homogenized computation; dx-xd=h^2 */
! 570: /* offset of h-degree */
! 571: NEWDL(d,n); d->d[n-1] = d0->d[n-1]+d1->d[n-1]; d->td = d->d[n-1];
! 572: NEWMP(mr); mr->c = (P)ONE; mr->dl = d;
! 573: MKDP(n,mr,r); r->sugar = d->d[n-1];
! 574:
! 575: for ( i = 0; i < n2; i++ ) {
! 576: a = d0->d[i]; b = d1->d[n2+i];
! 577: k = d0->d[n2+i]; l = d1->d[i];
! 578: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 579: h = a+k+l+b;
! 580: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 581: min = MIN(k,l);
! 582:
! 583: if ( min+1 > tablen ) {
! 584: if ( tab ) GC_free(tab);
! 585: tab = (Q *)MALLOC((min+1)*sizeof(Q));
! 586: tablen = min+1;
! 587: }
! 588: mkwc(k,l,tab);
! 589: for ( mr0 = 0, j = 0; j <= min; j++ ) {
! 590: NEXTMP(mr0,mr);
! 591: NEWDL(d,n);
! 592: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 593: d->td = h;
! 594: d->d[n-1] = h-(d->d[i]+d->d[n2+i]);
! 595: mr->c = (P)tab[j];
! 596: mr->dl = d;
! 597: }
! 598: bzero(tab,(min+1)*sizeof(Q));
! 599: if ( mr0 )
! 600: NEXT(mr) = 0;
! 601: MKDP(n,mr0,t);
! 602: if ( t )
! 603: t->sugar = h;
! 604: comm_muld(vl,r,t,&t1); r = t1;
! 605: }
! 606: } else
! 607: for ( i = 0, r = (DP)ONE; i < n2; i++ ) {
! 608: a = d0->d[i]; b = d1->d[n2+i];
! 609: k = d0->d[n2+i]; l = d1->d[i];
! 610: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 611: min = MIN(k,l);
! 612: tab = (Q *)ALLOCA((min+1)*sizeof(Q));
! 613: mkwc(k,l,tab);
! 614: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
! 615: NEXTMP(mr0,mr);
! 616: NEWDL(d,n);
! 617: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 618: d->td = d->d[i]+d->d[n2+i]; /* XXX */
! 619: s = MAX(s,d->td); /* XXX */
! 620: mr->c = (P)tab[j];
! 621: mr->dl = d;
! 622: }
! 623: if ( mr0 )
! 624: NEXT(mr) = 0;
! 625: MKDP(n,mr0,t);
! 626: if ( t )
! 627: t->sugar = s;
! 628: comm_muld(vl,r,t,&t1); r = t1;
! 629: }
1.2 noro 630: muldc(vl,r,c,pr);
1.1 noro 631: }
632: }
633:
634: void muldc(vl,p,c,pr)
635: VL vl;
636: DP p;
637: P c;
638: DP *pr;
639: {
640: MP m,mr,mr0;
641:
642: if ( !p || !c )
643: *pr = 0;
644: else if ( NUM(c) && UNIQ((Q)c) )
645: *pr = p;
646: else if ( NUM(c) && MUNIQ((Q)c) )
647: chsgnd(p,pr);
648: else {
649: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
650: NEXTMP(mr0,mr);
651: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
652: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
653: else
654: mulp(vl,C(m),c,&C(mr));
655: mr->dl = m->dl;
656: }
657: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
658: if ( *pr )
659: (*pr)->sugar = p->sugar;
660: }
661: }
662:
663: void divsdc(vl,p,c,pr)
664: VL vl;
665: DP p;
666: P c;
667: DP *pr;
668: {
669: MP m,mr,mr0;
670:
671: if ( !c )
672: error("disvsdc : division by 0");
673: else if ( !p )
674: *pr = 0;
675: else {
676: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
677: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
678: }
679: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
680: if ( *pr )
681: (*pr)->sugar = p->sugar;
682: }
683: }
684:
685: void adddl(n,d1,d2,dr)
686: int n;
687: DL d1,d2;
688: DL *dr;
689: {
690: DL dt;
691: int i;
692:
693: if ( !d1->td )
694: *dr = d2;
695: else if ( !d2->td )
696: *dr = d1;
697: else {
698: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
699: dt->td = d1->td + d2->td;
700: for ( i = 0; i < n; i++ )
701: dt->d[i] = d1->d[i]+d2->d[i];
702: }
703: }
704:
705: int compd(vl,p1,p2)
706: VL vl;
707: DP p1,p2;
708: {
709: int n,t;
710: MP m1,m2;
711:
712: if ( !p1 )
713: return p2 ? -1 : 0;
714: else if ( !p2 )
715: return 1;
716: else {
717: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
718: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
719: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
720: (t = compp(vl,C(m1),C(m2)) ) )
721: return t;
722: if ( m1 )
723: return 1;
724: else if ( m2 )
725: return -1;
726: else
727: return 0;
728: }
729: }
730:
731: int cmpdl_lex(n,d1,d2)
732: int n;
733: DL d1,d2;
734: {
735: int i;
736:
737: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
738: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
739: }
740:
741: int cmpdl_revlex(n,d1,d2)
742: int n;
743: DL d1,d2;
744: {
745: int i;
746:
747: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
748: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
749: }
750:
751: int cmpdl_gradlex(n,d1,d2)
752: int n;
753: DL d1,d2;
754: {
755: if ( d1->td > d2->td )
756: return 1;
757: else if ( d1->td < d2->td )
758: return -1;
759: else
760: return cmpdl_lex(n,d1,d2);
761: }
762:
763: int cmpdl_revgradlex(n,d1,d2)
764: int n;
765: DL d1,d2;
766: {
767: if ( d1->td > d2->td )
768: return 1;
769: else if ( d1->td < d2->td )
770: return -1;
771: else
772: return cmpdl_revlex(n,d1,d2);
773: }
774:
775: int cmpdl_blex(n,d1,d2)
776: int n;
777: DL d1,d2;
778: {
779: int c;
780:
781: if ( c = cmpdl_lex(n-1,d1,d2) )
782: return c;
783: else {
784: c = d1->d[n-1] - d2->d[n-1];
785: return c > 0 ? 1 : c < 0 ? -1 : 0;
786: }
787: }
788:
789: int cmpdl_bgradlex(n,d1,d2)
790: int n;
791: DL d1,d2;
792: {
793: int e1,e2,c;
794:
795: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
796: if ( e1 > e2 )
797: return 1;
798: else if ( e1 < e2 )
799: return -1;
800: else {
801: c = cmpdl_lex(n-1,d1,d2);
802: if ( c )
803: return c;
804: else
805: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
806: }
807: }
808:
809: int cmpdl_brevgradlex(n,d1,d2)
810: int n;
811: DL d1,d2;
812: {
813: int e1,e2,c;
814:
815: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
816: if ( e1 > e2 )
817: return 1;
818: else if ( e1 < e2 )
819: return -1;
820: else {
821: c = cmpdl_revlex(n-1,d1,d2);
822: if ( c )
823: return c;
824: else
825: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
826: }
827: }
828:
829: int cmpdl_brevrev(n,d1,d2)
830: int n;
831: DL d1,d2;
832: {
833: int e1,e2,f1,f2,c,i;
834:
835: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
836: e1 += d1->d[i]; e2 += d2->d[i];
837: }
838: f1 = d1->td - e1; f2 = d2->td - e2;
839: if ( e1 > e2 )
840: return 1;
841: else if ( e1 < e2 )
842: return -1;
843: else {
844: c = cmpdl_revlex(dp_nelim,d1,d2);
845: if ( c )
846: return c;
847: else if ( f1 > f2 )
848: return 1;
849: else if ( f1 < f2 )
850: return -1;
851: else {
852: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
853: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
854: }
855: }
856: }
857:
858: int cmpdl_bgradrev(n,d1,d2)
859: int n;
860: DL d1,d2;
861: {
862: int e1,e2,f1,f2,c,i;
863:
864: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
865: e1 += d1->d[i]; e2 += d2->d[i];
866: }
867: f1 = d1->td - e1; f2 = d2->td - e2;
868: if ( e1 > e2 )
869: return 1;
870: else if ( e1 < e2 )
871: return -1;
872: else {
873: c = cmpdl_lex(dp_nelim,d1,d2);
874: if ( c )
875: return c;
876: else if ( f1 > f2 )
877: return 1;
878: else if ( f1 < f2 )
879: return -1;
880: else {
881: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
882: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
883: }
884: }
885: }
886:
887: int cmpdl_blexrev(n,d1,d2)
888: int n;
889: DL d1,d2;
890: {
891: int e1,e2,f1,f2,c,i;
892:
893: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
894: e1 += d1->d[i]; e2 += d2->d[i];
895: }
896: f1 = d1->td - e1; f2 = d2->td - e2;
897: c = cmpdl_lex(dp_nelim,d1,d2);
898: if ( c )
899: return c;
900: else if ( f1 > f2 )
901: return 1;
902: else if ( f1 < f2 )
903: return -1;
904: else {
905: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
906: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
907: }
908: }
909:
910: int cmpdl_elim(n,d1,d2)
911: int n;
912: DL d1,d2;
913: {
914: int e1,e2,i;
915:
916: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
917: e1 += d1->d[i]; e2 += d2->d[i];
918: }
919: if ( e1 > e2 )
920: return 1;
921: else if ( e1 < e2 )
922: return -1;
923: else
924: return cmpdl_revgradlex(n,d1,d2);
925: }
926:
927: int cmpdl_order_pair(n,d1,d2)
928: int n;
929: DL d1,d2;
930: {
931: int e1,e2,i,j,l;
932: int *t1,*t2;
933: int len;
934: struct order_pair *pair;
935:
936: len = dp_current_spec.ord.block.length;
937: pair = dp_current_spec.ord.block.order_pair;
938:
939: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
940: l = pair[i].length;
941: switch ( pair[i].order ) {
942: case 0:
943: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
944: e1 += t1[j]; e2 += t2[j];
945: }
946: if ( e1 > e2 )
947: return 1;
948: else if ( e1 < e2 )
949: return -1;
950: else {
951: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
952: if ( j >= 0 )
953: return t1[j] < t2[j] ? 1 : -1;
954: }
955: break;
956: case 1:
957: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
958: e1 += t1[j]; e2 += t2[j];
959: }
960: if ( e1 > e2 )
961: return 1;
962: else if ( e1 < e2 )
963: return -1;
964: else {
965: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
966: if ( j < l )
967: return t1[j] > t2[j] ? 1 : -1;
968: }
969: break;
970: case 2:
971: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
972: if ( j < l )
973: return t1[j] > t2[j] ? 1 : -1;
974: break;
975: default:
976: error("cmpdl_order_pair : invalid order"); break;
977: }
978: t1 += l; t2 += l;
979: }
980: return 0;
981: }
982:
983: int cmpdl_matrix(n,d1,d2)
984: int n;
985: DL d1,d2;
986: {
987: int *v,*w,*t1,*t2;
988: int s,i,j,len;
989: int **matrix;
990:
991: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
992: w[i] = t1[i]-t2[i];
993: len = dp_current_spec.ord.matrix.row;
994: matrix = dp_current_spec.ord.matrix.matrix;
995: for ( j = 0; j < len; j++ ) {
996: v = matrix[j];
997: for ( i = 0, s = 0; i < n; i++ )
998: s += v[i]*w[i];
999: if ( s > 0 )
1000: return 1;
1001: else if ( s < 0 )
1002: return -1;
1003: }
1004: return 0;
1005: }