![[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.6 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.5 2000/05/29 08:54:46 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.6 ! noro 557: int i,j,a,b,k,l,n2,s,min;
1.5 noro 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 */
1.6 ! noro 571: NEWDL(d,n);
! 572: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1.5 noro 573: NEWMP(mr); mr->c = (P)ONE; mr->dl = d;
1.6 ! noro 574: MKDP(n,mr,r); r->sugar = d->td;
! 575: } else
! 576: r = (DP)ONE;
! 577: for ( i = 0; i < n2; i++ ) {
! 578: a = d0->d[i]; b = d1->d[n2+i];
! 579: k = d0->d[n2+i]; l = d1->d[i];
! 580: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 581: s = a+k+l+b;
! 582: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 583: min = MIN(k,l);
! 584:
! 585: if ( min+1 > tablen ) {
! 586: if ( tab ) GC_free(tab);
! 587: tab = (Q *)MALLOC((min+1)*sizeof(Q));
! 588: tablen = min+1;
! 589: }
! 590: mkwc(k,l,tab);
! 591: if ( n & 1 )
1.5 noro 592: for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.6 ! noro 593: NEXTMP(mr0,mr); NEWDL(d,n);
1.5 noro 594: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.6 ! noro 595: d->td = s;
! 596: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
! 597: mr->c = (P)tab[j]; mr->dl = d;
1.5 noro 598: }
1.6 ! noro 599: else
1.5 noro 600: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.6 ! noro 601: NEXTMP(mr0,mr); NEWDL(d,n);
1.5 noro 602: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
603: d->td = d->d[i]+d->d[n2+i]; /* XXX */
604: s = MAX(s,d->td); /* XXX */
1.6 ! noro 605: mr->c = (P)tab[j]; mr->dl = d;
1.5 noro 606: }
1.6 ! noro 607: bzero(tab,(min+1)*sizeof(Q));
! 608: if ( mr0 )
! 609: NEXT(mr) = 0;
! 610: MKDP(n,mr0,t);
! 611: if ( t )
! 612: t->sugar = s;
! 613: comm_muld(vl,r,t,&t1); r = t1;
! 614: }
1.2 noro 615: muldc(vl,r,c,pr);
1.1 noro 616: }
617: }
618:
619: void muldc(vl,p,c,pr)
620: VL vl;
621: DP p;
622: P c;
623: DP *pr;
624: {
625: MP m,mr,mr0;
626:
627: if ( !p || !c )
628: *pr = 0;
629: else if ( NUM(c) && UNIQ((Q)c) )
630: *pr = p;
631: else if ( NUM(c) && MUNIQ((Q)c) )
632: chsgnd(p,pr);
633: else {
634: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
635: NEXTMP(mr0,mr);
636: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
637: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
638: else
639: mulp(vl,C(m),c,&C(mr));
640: mr->dl = m->dl;
641: }
642: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
643: if ( *pr )
644: (*pr)->sugar = p->sugar;
645: }
646: }
647:
648: void divsdc(vl,p,c,pr)
649: VL vl;
650: DP p;
651: P c;
652: DP *pr;
653: {
654: MP m,mr,mr0;
655:
656: if ( !c )
657: error("disvsdc : division by 0");
658: else if ( !p )
659: *pr = 0;
660: else {
661: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
662: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
663: }
664: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
665: if ( *pr )
666: (*pr)->sugar = p->sugar;
667: }
668: }
669:
670: void adddl(n,d1,d2,dr)
671: int n;
672: DL d1,d2;
673: DL *dr;
674: {
675: DL dt;
676: int i;
677:
678: if ( !d1->td )
679: *dr = d2;
680: else if ( !d2->td )
681: *dr = d1;
682: else {
683: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
684: dt->td = d1->td + d2->td;
685: for ( i = 0; i < n; i++ )
686: dt->d[i] = d1->d[i]+d2->d[i];
687: }
688: }
689:
690: int compd(vl,p1,p2)
691: VL vl;
692: DP p1,p2;
693: {
694: int n,t;
695: MP m1,m2;
696:
697: if ( !p1 )
698: return p2 ? -1 : 0;
699: else if ( !p2 )
700: return 1;
701: else {
702: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
703: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
704: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
705: (t = compp(vl,C(m1),C(m2)) ) )
706: return t;
707: if ( m1 )
708: return 1;
709: else if ( m2 )
710: return -1;
711: else
712: return 0;
713: }
714: }
715:
716: int cmpdl_lex(n,d1,d2)
717: int n;
718: DL d1,d2;
719: {
720: int i;
721:
722: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
723: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
724: }
725:
726: int cmpdl_revlex(n,d1,d2)
727: int n;
728: DL d1,d2;
729: {
730: int i;
731:
732: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
733: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
734: }
735:
736: int cmpdl_gradlex(n,d1,d2)
737: int n;
738: DL d1,d2;
739: {
740: if ( d1->td > d2->td )
741: return 1;
742: else if ( d1->td < d2->td )
743: return -1;
744: else
745: return cmpdl_lex(n,d1,d2);
746: }
747:
748: int cmpdl_revgradlex(n,d1,d2)
749: int n;
750: DL d1,d2;
751: {
752: if ( d1->td > d2->td )
753: return 1;
754: else if ( d1->td < d2->td )
755: return -1;
756: else
757: return cmpdl_revlex(n,d1,d2);
758: }
759:
760: int cmpdl_blex(n,d1,d2)
761: int n;
762: DL d1,d2;
763: {
764: int c;
765:
766: if ( c = cmpdl_lex(n-1,d1,d2) )
767: return c;
768: else {
769: c = d1->d[n-1] - d2->d[n-1];
770: return c > 0 ? 1 : c < 0 ? -1 : 0;
771: }
772: }
773:
774: int cmpdl_bgradlex(n,d1,d2)
775: int n;
776: DL d1,d2;
777: {
778: int e1,e2,c;
779:
780: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
781: if ( e1 > e2 )
782: return 1;
783: else if ( e1 < e2 )
784: return -1;
785: else {
786: c = cmpdl_lex(n-1,d1,d2);
787: if ( c )
788: return c;
789: else
790: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
791: }
792: }
793:
794: int cmpdl_brevgradlex(n,d1,d2)
795: int n;
796: DL d1,d2;
797: {
798: int e1,e2,c;
799:
800: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
801: if ( e1 > e2 )
802: return 1;
803: else if ( e1 < e2 )
804: return -1;
805: else {
806: c = cmpdl_revlex(n-1,d1,d2);
807: if ( c )
808: return c;
809: else
810: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
811: }
812: }
813:
814: int cmpdl_brevrev(n,d1,d2)
815: int n;
816: DL d1,d2;
817: {
818: int e1,e2,f1,f2,c,i;
819:
820: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
821: e1 += d1->d[i]; e2 += d2->d[i];
822: }
823: f1 = d1->td - e1; f2 = d2->td - e2;
824: if ( e1 > e2 )
825: return 1;
826: else if ( e1 < e2 )
827: return -1;
828: else {
829: c = cmpdl_revlex(dp_nelim,d1,d2);
830: if ( c )
831: return c;
832: else if ( f1 > f2 )
833: return 1;
834: else if ( f1 < f2 )
835: return -1;
836: else {
837: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
838: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
839: }
840: }
841: }
842:
843: int cmpdl_bgradrev(n,d1,d2)
844: int n;
845: DL d1,d2;
846: {
847: int e1,e2,f1,f2,c,i;
848:
849: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
850: e1 += d1->d[i]; e2 += d2->d[i];
851: }
852: f1 = d1->td - e1; f2 = d2->td - e2;
853: if ( e1 > e2 )
854: return 1;
855: else if ( e1 < e2 )
856: return -1;
857: else {
858: c = cmpdl_lex(dp_nelim,d1,d2);
859: if ( c )
860: return c;
861: else if ( f1 > f2 )
862: return 1;
863: else if ( f1 < f2 )
864: return -1;
865: else {
866: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
867: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
868: }
869: }
870: }
871:
872: int cmpdl_blexrev(n,d1,d2)
873: int n;
874: DL d1,d2;
875: {
876: int e1,e2,f1,f2,c,i;
877:
878: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
879: e1 += d1->d[i]; e2 += d2->d[i];
880: }
881: f1 = d1->td - e1; f2 = d2->td - e2;
882: c = cmpdl_lex(dp_nelim,d1,d2);
883: if ( c )
884: return c;
885: else if ( f1 > f2 )
886: return 1;
887: else if ( f1 < f2 )
888: return -1;
889: else {
890: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
891: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
892: }
893: }
894:
895: int cmpdl_elim(n,d1,d2)
896: int n;
897: DL d1,d2;
898: {
899: int e1,e2,i;
900:
901: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
902: e1 += d1->d[i]; e2 += d2->d[i];
903: }
904: if ( e1 > e2 )
905: return 1;
906: else if ( e1 < e2 )
907: return -1;
908: else
909: return cmpdl_revgradlex(n,d1,d2);
910: }
911:
912: int cmpdl_order_pair(n,d1,d2)
913: int n;
914: DL d1,d2;
915: {
916: int e1,e2,i,j,l;
917: int *t1,*t2;
918: int len;
919: struct order_pair *pair;
920:
921: len = dp_current_spec.ord.block.length;
922: pair = dp_current_spec.ord.block.order_pair;
923:
924: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
925: l = pair[i].length;
926: switch ( pair[i].order ) {
927: case 0:
928: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
929: e1 += t1[j]; e2 += t2[j];
930: }
931: if ( e1 > e2 )
932: return 1;
933: else if ( e1 < e2 )
934: return -1;
935: else {
936: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
937: if ( j >= 0 )
938: return t1[j] < t2[j] ? 1 : -1;
939: }
940: break;
941: case 1:
942: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
943: e1 += t1[j]; e2 += t2[j];
944: }
945: if ( e1 > e2 )
946: return 1;
947: else if ( e1 < e2 )
948: return -1;
949: else {
950: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
951: if ( j < l )
952: return t1[j] > t2[j] ? 1 : -1;
953: }
954: break;
955: case 2:
956: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
957: if ( j < l )
958: return t1[j] > t2[j] ? 1 : -1;
959: break;
960: default:
961: error("cmpdl_order_pair : invalid order"); break;
962: }
963: t1 += l; t2 += l;
964: }
965: return 0;
966: }
967:
968: int cmpdl_matrix(n,d1,d2)
969: int n;
970: DL d1,d2;
971: {
972: int *v,*w,*t1,*t2;
973: int s,i,j,len;
974: int **matrix;
975:
976: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
977: w[i] = t1[i]-t2[i];
978: len = dp_current_spec.ord.matrix.row;
979: matrix = dp_current_spec.ord.matrix.matrix;
980: for ( j = 0; j < len; j++ ) {
981: v = matrix[j];
982: for ( i = 0, s = 0; i < n; i++ )
983: s += v[i]*w[i];
984: if ( s > 0 )
985: return 1;
986: else if ( s < 0 )
987: return -1;
988: }
989: return 0;
990: }