![[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.7 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.6 2000/05/30 01:35:12 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: {
1.7 ! noro 752: register int i;
! 753: register int *p1,*p2;
! 754:
1.1 noro 755: if ( d1->td > d2->td )
756: return 1;
757: else if ( d1->td < d2->td )
758: return -1;
1.7 ! noro 759: else {
! 760: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
! 761: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
! 762: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
! 763: }
1.1 noro 764: }
765:
766: int cmpdl_blex(n,d1,d2)
767: int n;
768: DL d1,d2;
769: {
770: int c;
771:
772: if ( c = cmpdl_lex(n-1,d1,d2) )
773: return c;
774: else {
775: c = d1->d[n-1] - d2->d[n-1];
776: return c > 0 ? 1 : c < 0 ? -1 : 0;
777: }
778: }
779:
780: int cmpdl_bgradlex(n,d1,d2)
781: int n;
782: DL d1,d2;
783: {
784: int e1,e2,c;
785:
786: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
787: if ( e1 > e2 )
788: return 1;
789: else if ( e1 < e2 )
790: return -1;
791: else {
792: c = cmpdl_lex(n-1,d1,d2);
793: if ( c )
794: return c;
795: else
796: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
797: }
798: }
799:
800: int cmpdl_brevgradlex(n,d1,d2)
801: int n;
802: DL d1,d2;
803: {
804: int e1,e2,c;
805:
806: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
807: if ( e1 > e2 )
808: return 1;
809: else if ( e1 < e2 )
810: return -1;
811: else {
812: c = cmpdl_revlex(n-1,d1,d2);
813: if ( c )
814: return c;
815: else
816: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
817: }
818: }
819:
820: int cmpdl_brevrev(n,d1,d2)
821: int n;
822: DL d1,d2;
823: {
824: int e1,e2,f1,f2,c,i;
825:
826: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
827: e1 += d1->d[i]; e2 += d2->d[i];
828: }
829: f1 = d1->td - e1; f2 = d2->td - e2;
830: if ( e1 > e2 )
831: return 1;
832: else if ( e1 < e2 )
833: return -1;
834: else {
835: c = cmpdl_revlex(dp_nelim,d1,d2);
836: if ( c )
837: return c;
838: else if ( f1 > f2 )
839: return 1;
840: else if ( f1 < f2 )
841: return -1;
842: else {
843: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
844: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
845: }
846: }
847: }
848:
849: int cmpdl_bgradrev(n,d1,d2)
850: int n;
851: DL d1,d2;
852: {
853: int e1,e2,f1,f2,c,i;
854:
855: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
856: e1 += d1->d[i]; e2 += d2->d[i];
857: }
858: f1 = d1->td - e1; f2 = d2->td - e2;
859: if ( e1 > e2 )
860: return 1;
861: else if ( e1 < e2 )
862: return -1;
863: else {
864: c = cmpdl_lex(dp_nelim,d1,d2);
865: if ( c )
866: return c;
867: else if ( f1 > f2 )
868: return 1;
869: else if ( f1 < f2 )
870: return -1;
871: else {
872: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
873: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
874: }
875: }
876: }
877:
878: int cmpdl_blexrev(n,d1,d2)
879: int n;
880: DL d1,d2;
881: {
882: int e1,e2,f1,f2,c,i;
883:
884: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
885: e1 += d1->d[i]; e2 += d2->d[i];
886: }
887: f1 = d1->td - e1; f2 = d2->td - e2;
888: c = cmpdl_lex(dp_nelim,d1,d2);
889: if ( c )
890: return c;
891: else if ( f1 > f2 )
892: return 1;
893: else if ( f1 < f2 )
894: return -1;
895: else {
896: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
897: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
898: }
899: }
900:
901: int cmpdl_elim(n,d1,d2)
902: int n;
903: DL d1,d2;
904: {
905: int e1,e2,i;
906:
907: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
908: e1 += d1->d[i]; e2 += d2->d[i];
909: }
910: if ( e1 > e2 )
911: return 1;
912: else if ( e1 < e2 )
913: return -1;
914: else
915: return cmpdl_revgradlex(n,d1,d2);
916: }
917:
918: int cmpdl_order_pair(n,d1,d2)
919: int n;
920: DL d1,d2;
921: {
922: int e1,e2,i,j,l;
923: int *t1,*t2;
924: int len;
925: struct order_pair *pair;
926:
927: len = dp_current_spec.ord.block.length;
928: pair = dp_current_spec.ord.block.order_pair;
929:
930: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
931: l = pair[i].length;
932: switch ( pair[i].order ) {
933: case 0:
934: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
935: e1 += t1[j]; e2 += t2[j];
936: }
937: if ( e1 > e2 )
938: return 1;
939: else if ( e1 < e2 )
940: return -1;
941: else {
942: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
943: if ( j >= 0 )
944: return t1[j] < t2[j] ? 1 : -1;
945: }
946: break;
947: case 1:
948: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
949: e1 += t1[j]; e2 += t2[j];
950: }
951: if ( e1 > e2 )
952: return 1;
953: else if ( e1 < e2 )
954: return -1;
955: else {
956: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
957: if ( j < l )
958: return t1[j] > t2[j] ? 1 : -1;
959: }
960: break;
961: case 2:
962: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
963: if ( j < l )
964: return t1[j] > t2[j] ? 1 : -1;
965: break;
966: default:
967: error("cmpdl_order_pair : invalid order"); break;
968: }
969: t1 += l; t2 += l;
970: }
971: return 0;
972: }
973:
974: int cmpdl_matrix(n,d1,d2)
975: int n;
976: DL d1,d2;
977: {
978: int *v,*w,*t1,*t2;
979: int s,i,j,len;
980: int **matrix;
981:
982: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
983: w[i] = t1[i]-t2[i];
984: len = dp_current_spec.ord.matrix.row;
985: matrix = dp_current_spec.ord.matrix.matrix;
986: for ( j = 0; j < len; j++ ) {
987: v = matrix[j];
988: for ( i = 0, s = 0; i < n; i++ )
989: s += v[i]*w[i];
990: if ( s > 0 )
991: return 1;
992: else if ( s < 0 )
993: return -1;
994: }
995: return 0;
996: }