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