Annotation of OpenXM_contrib2/asir2000/engine/dist.c, Revision 1.4
1.4 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.3 2000/04/13 06:01:02 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.4 ! noro 417: int i,l;
! 418: MP *w;
1.1 noro 419:
420: if ( !p1 || !p2 )
421: *pr = 0;
422: else if ( OID(p1) <= O_P )
423: muldc(vl,p2,(P)p1,pr);
424: else if ( OID(p2) <= O_P )
425: muldc(vl,p1,(P)p2,pr);
426: else {
1.4 ! noro 427: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 428: w = ALLOCA(l*sizeof(MP));
! 429: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 430: w[i] = m;
! 431: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 432: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1 noro 433: }
434: *pr = s;
435: }
436: }
437:
438: void muldm(vl,p,m0,pr)
439: VL vl;
440: DP p;
441: MP m0;
442: DP *pr;
443: {
444: MP m,mr,mr0;
445: P c;
446: DL d;
447: int n;
448:
449: if ( !p )
450: *pr = 0;
451: else {
452: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
453: m; m = NEXT(m) ) {
454: NEXTMP(mr0,mr);
455: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
456: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
457: else
458: mulp(vl,C(m),c,&C(mr));
459: adddl(n,m->dl,d,&mr->dl);
460: }
461: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
462: if ( *pr )
463: (*pr)->sugar = p->sugar + m0->dl->td;
1.2 noro 464: }
465: }
466:
467: void weyl_muld(vl,p1,p2,pr)
468: VL vl;
469: DP p1,p2,*pr;
470: {
471: MP m;
472: DP s,t,u;
1.4 ! noro 473: int i,l;
! 474: MP *w;
1.2 noro 475:
476: if ( !p1 || !p2 )
477: *pr = 0;
478: else if ( OID(p1) <= O_P )
479: muldc(vl,p2,(P)p1,pr);
480: else if ( OID(p2) <= O_P )
481: muldc(vl,p1,(P)p2,pr);
482: else {
1.4 ! noro 483: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 484: w = ALLOCA(l*sizeof(MP));
! 485: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 486: w[i] = m;
! 487: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 488: weyl_muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.2 noro 489: }
490: *pr = s;
491: }
492: }
493:
494: void weyl_muldm(vl,p,m0,pr)
495: VL vl;
496: DP p;
497: MP m0;
498: DP *pr;
499: {
500: DP r,t,t1;
501: MP m;
1.4 ! noro 502: int n,l,i;
! 503: MP *w;
1.2 noro 504:
505: if ( !p )
506: *pr = 0;
507: else {
1.4 ! noro 508: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 509: w = ALLOCA(l*sizeof(MP));
! 510: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 511: w[i] = m;
! 512: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
! 513: weyl_mulmm(vl,w[i],m0,n,&t);
1.2 noro 514: addd(vl,r,t,&t1); r = t1;
515: }
516: if ( r )
517: r->sugar = p->sugar + m0->dl->td;
518: *pr = r;
519: }
520: }
521:
522: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
523:
524: void weyl_mulmm(vl,m0,m1,n,pr)
525: VL vl;
526: MP m0,m1;
527: int n;
528: DP *pr;
529: {
530: MP m,mr,mr0;
531: DP r,t,t1;
532: P c,c0,c1,cc;
533: DL d,d0,d1;
534: int i,j,a,b,k,l,n2,s,min;
535: Q *tab;
536:
537: if ( !m0 || !m1 )
538: *pr = 0;
539: else {
540: c0 = C(m0); c1 = C(m1);
541: mulp(vl,c0,c1,&c);
542: d0 = m0->dl; d1 = m1->dl;
543: n2 = n>>1;
544: for ( i = 0, r = (DP)ONE; i < n2; i++ ) {
545: a = d0->d[i]; b = d1->d[n2+i];
546: k = d0->d[n2+i]; l = d1->d[i];
547: /* compute xi^a*(Di^k*xi^l)*Di^b */
548: min = MIN(k,l);
549: tab = (Q *)ALLOCA((min+1)*sizeof(Q));
550: mkwc(k,l,tab);
551: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
552: NEXTMP(mr0,mr);
553: NEWDL(d,n);
554: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
555: d->td = d->d[i]+d->d[n2+i]; /* XXX */
556: s = MAX(s,d->td); /* XXX */
557: mr->c = (P)tab[j];
558: mr->dl = d;
559: }
560: if ( mr0 )
561: NEXT(mr) = 0;
562: MKDP(n,mr0,t);
563: if ( t )
564: t->sugar = s;
565: comm_muld(vl,r,t,&t1); r = t1;
566: }
567: muldc(vl,r,c,pr);
1.1 noro 568: }
569: }
570:
571: void muldc(vl,p,c,pr)
572: VL vl;
573: DP p;
574: P c;
575: DP *pr;
576: {
577: MP m,mr,mr0;
578:
579: if ( !p || !c )
580: *pr = 0;
581: else if ( NUM(c) && UNIQ((Q)c) )
582: *pr = p;
583: else if ( NUM(c) && MUNIQ((Q)c) )
584: chsgnd(p,pr);
585: else {
586: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
587: NEXTMP(mr0,mr);
588: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
589: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
590: else
591: mulp(vl,C(m),c,&C(mr));
592: mr->dl = m->dl;
593: }
594: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
595: if ( *pr )
596: (*pr)->sugar = p->sugar;
597: }
598: }
599:
600: void divsdc(vl,p,c,pr)
601: VL vl;
602: DP p;
603: P c;
604: DP *pr;
605: {
606: MP m,mr,mr0;
607:
608: if ( !c )
609: error("disvsdc : division by 0");
610: else if ( !p )
611: *pr = 0;
612: else {
613: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
614: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
615: }
616: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
617: if ( *pr )
618: (*pr)->sugar = p->sugar;
619: }
620: }
621:
622: void adddl(n,d1,d2,dr)
623: int n;
624: DL d1,d2;
625: DL *dr;
626: {
627: DL dt;
628: int i;
629:
630: if ( !d1->td )
631: *dr = d2;
632: else if ( !d2->td )
633: *dr = d1;
634: else {
635: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
636: dt->td = d1->td + d2->td;
637: for ( i = 0; i < n; i++ )
638: dt->d[i] = d1->d[i]+d2->d[i];
639: }
640: }
641:
642: int compd(vl,p1,p2)
643: VL vl;
644: DP p1,p2;
645: {
646: int n,t;
647: MP m1,m2;
648:
649: if ( !p1 )
650: return p2 ? -1 : 0;
651: else if ( !p2 )
652: return 1;
653: else {
654: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
655: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
656: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
657: (t = compp(vl,C(m1),C(m2)) ) )
658: return t;
659: if ( m1 )
660: return 1;
661: else if ( m2 )
662: return -1;
663: else
664: return 0;
665: }
666: }
667:
668: int cmpdl_lex(n,d1,d2)
669: int n;
670: DL d1,d2;
671: {
672: int i;
673:
674: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
675: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
676: }
677:
678: int cmpdl_revlex(n,d1,d2)
679: int n;
680: DL d1,d2;
681: {
682: int i;
683:
684: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
685: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
686: }
687:
688: int cmpdl_gradlex(n,d1,d2)
689: int n;
690: DL d1,d2;
691: {
692: if ( d1->td > d2->td )
693: return 1;
694: else if ( d1->td < d2->td )
695: return -1;
696: else
697: return cmpdl_lex(n,d1,d2);
698: }
699:
700: int cmpdl_revgradlex(n,d1,d2)
701: int n;
702: DL d1,d2;
703: {
704: if ( d1->td > d2->td )
705: return 1;
706: else if ( d1->td < d2->td )
707: return -1;
708: else
709: return cmpdl_revlex(n,d1,d2);
710: }
711:
712: int cmpdl_blex(n,d1,d2)
713: int n;
714: DL d1,d2;
715: {
716: int c;
717:
718: if ( c = cmpdl_lex(n-1,d1,d2) )
719: return c;
720: else {
721: c = d1->d[n-1] - d2->d[n-1];
722: return c > 0 ? 1 : c < 0 ? -1 : 0;
723: }
724: }
725:
726: int cmpdl_bgradlex(n,d1,d2)
727: int n;
728: DL d1,d2;
729: {
730: int e1,e2,c;
731:
732: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
733: if ( e1 > e2 )
734: return 1;
735: else if ( e1 < e2 )
736: return -1;
737: else {
738: c = cmpdl_lex(n-1,d1,d2);
739: if ( c )
740: return c;
741: else
742: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
743: }
744: }
745:
746: int cmpdl_brevgradlex(n,d1,d2)
747: int n;
748: DL d1,d2;
749: {
750: int e1,e2,c;
751:
752: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
753: if ( e1 > e2 )
754: return 1;
755: else if ( e1 < e2 )
756: return -1;
757: else {
758: c = cmpdl_revlex(n-1,d1,d2);
759: if ( c )
760: return c;
761: else
762: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
763: }
764: }
765:
766: int cmpdl_brevrev(n,d1,d2)
767: int n;
768: DL d1,d2;
769: {
770: int e1,e2,f1,f2,c,i;
771:
772: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
773: e1 += d1->d[i]; e2 += d2->d[i];
774: }
775: f1 = d1->td - e1; f2 = d2->td - e2;
776: if ( e1 > e2 )
777: return 1;
778: else if ( e1 < e2 )
779: return -1;
780: else {
781: c = cmpdl_revlex(dp_nelim,d1,d2);
782: if ( c )
783: return c;
784: else if ( f1 > f2 )
785: return 1;
786: else if ( f1 < f2 )
787: return -1;
788: else {
789: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
790: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
791: }
792: }
793: }
794:
795: int cmpdl_bgradrev(n,d1,d2)
796: int n;
797: DL d1,d2;
798: {
799: int e1,e2,f1,f2,c,i;
800:
801: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
802: e1 += d1->d[i]; e2 += d2->d[i];
803: }
804: f1 = d1->td - e1; f2 = d2->td - e2;
805: if ( e1 > e2 )
806: return 1;
807: else if ( e1 < e2 )
808: return -1;
809: else {
810: c = cmpdl_lex(dp_nelim,d1,d2);
811: if ( c )
812: return c;
813: else if ( f1 > f2 )
814: return 1;
815: else if ( f1 < f2 )
816: return -1;
817: else {
818: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
819: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
820: }
821: }
822: }
823:
824: int cmpdl_blexrev(n,d1,d2)
825: int n;
826: DL d1,d2;
827: {
828: int e1,e2,f1,f2,c,i;
829:
830: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
831: e1 += d1->d[i]; e2 += d2->d[i];
832: }
833: f1 = d1->td - e1; f2 = d2->td - e2;
834: c = cmpdl_lex(dp_nelim,d1,d2);
835: if ( c )
836: return c;
837: else if ( f1 > f2 )
838: return 1;
839: else if ( f1 < f2 )
840: return -1;
841: else {
842: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
843: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
844: }
845: }
846:
847: int cmpdl_elim(n,d1,d2)
848: int n;
849: DL d1,d2;
850: {
851: int e1,e2,i;
852:
853: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
854: e1 += d1->d[i]; e2 += d2->d[i];
855: }
856: if ( e1 > e2 )
857: return 1;
858: else if ( e1 < e2 )
859: return -1;
860: else
861: return cmpdl_revgradlex(n,d1,d2);
862: }
863:
864: int cmpdl_order_pair(n,d1,d2)
865: int n;
866: DL d1,d2;
867: {
868: int e1,e2,i,j,l;
869: int *t1,*t2;
870: int len;
871: struct order_pair *pair;
872:
873: len = dp_current_spec.ord.block.length;
874: pair = dp_current_spec.ord.block.order_pair;
875:
876: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
877: l = pair[i].length;
878: switch ( pair[i].order ) {
879: case 0:
880: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
881: e1 += t1[j]; e2 += t2[j];
882: }
883: if ( e1 > e2 )
884: return 1;
885: else if ( e1 < e2 )
886: return -1;
887: else {
888: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
889: if ( j >= 0 )
890: return t1[j] < t2[j] ? 1 : -1;
891: }
892: break;
893: case 1:
894: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
895: e1 += t1[j]; e2 += t2[j];
896: }
897: if ( e1 > e2 )
898: return 1;
899: else if ( e1 < e2 )
900: return -1;
901: else {
902: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
903: if ( j < l )
904: return t1[j] > t2[j] ? 1 : -1;
905: }
906: break;
907: case 2:
908: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
909: if ( j < l )
910: return t1[j] > t2[j] ? 1 : -1;
911: break;
912: default:
913: error("cmpdl_order_pair : invalid order"); break;
914: }
915: t1 += l; t2 += l;
916: }
917: return 0;
918: }
919:
920: int cmpdl_matrix(n,d1,d2)
921: int n;
922: DL d1,d2;
923: {
924: int *v,*w,*t1,*t2;
925: int s,i,j,len;
926: int **matrix;
927:
928: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
929: w[i] = t1[i]-t2[i];
930: len = dp_current_spec.ord.matrix.row;
931: matrix = dp_current_spec.ord.matrix.matrix;
932: for ( j = 0; j < len; j++ ) {
933: v = matrix[j];
934: for ( i = 0, s = 0; i < n; i++ )
935: s += v[i]*w[i];
936: if ( s > 0 )
937: return 1;
938: else if ( s < 0 )
939: return -1;
940: }
941: return 0;
942: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>