Annotation of OpenXM_contrib2/asir2000/engine/dist.c, Revision 1.1.1.1
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/dist.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
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:
21: int dp_nelim,dp_fcoeffs;
22: struct order_spec dp_current_spec;
23: int *dp_dl_work;
24:
25: int has_fcoef(DP);
26: int has_fcoef_p(P);
27:
28: int has_fcoef(f)
29: DP f;
30: {
31: MP t;
32:
33: if ( !f )
34: return 0;
35: for ( t = BDY(f); t; t = NEXT(t) )
36: if ( has_fcoef_p(t->c) )
37: break;
38: return t ? 1 : 0;
39: }
40:
41: int has_fcoef_p(f)
42: P f;
43: {
44: DCP dc;
45:
46: if ( !f )
47: return 0;
48: else if ( NUM(f) )
49: return (NID((Num)f) == N_LM || NID((Num)f) == N_GF2N) ? 1 : 0;
50: else {
51: for ( dc = DC(f); dc; dc = NEXT(dc) )
52: if ( has_fcoef_p(COEF(dc)) )
53: return 1;
54: return 0;
55: }
56: }
57:
58: void initd(spec)
59: struct order_spec *spec;
60: {
61: switch ( spec->id ) {
62: case 2:
63: cmpdl = cmpdl_matrix;
64: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
65: break;
66: case 1:
67: cmpdl = cmpdl_order_pair;
68: break;
69: default:
70: switch ( spec->ord.simple ) {
71: case ORD_REVGRADLEX:
72: cmpdl = cmpdl_revgradlex; break;
73: case ORD_GRADLEX:
74: cmpdl = cmpdl_gradlex; break;
75: case ORD_BREVGRADLEX:
76: cmpdl = cmpdl_brevgradlex; break;
77: case ORD_BGRADLEX:
78: cmpdl = cmpdl_bgradlex; break;
79: case ORD_BLEX:
80: cmpdl = cmpdl_blex; break;
81: case ORD_BREVREV:
82: cmpdl = cmpdl_brevrev; break;
83: case ORD_BGRADREV:
84: cmpdl = cmpdl_bgradrev; break;
85: case ORD_BLEXREV:
86: cmpdl = cmpdl_blexrev; break;
87: case ORD_ELIM:
88: cmpdl = cmpdl_elim; break;
89: case ORD_LEX: default:
90: cmpdl = cmpdl_lex; break;
91: }
92: break;
93: }
94: dp_current_spec = *spec;
95: }
96:
97: void ptod(vl,dvl,p,pr)
98: VL vl,dvl;
99: P p;
100: DP *pr;
101: {
102: int isconst = 0;
103: int n,i;
104: VL tvl;
105: V v;
106: DL d;
107: MP m;
108: DCP dc;
109: DP r,s,t,u;
110: P x,c;
111:
112: if ( !p )
113: *pr = 0;
114: else {
115: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
116: if ( NUM(p) ) {
117: NEWDL(d,n);
118: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
119: } else {
120: for ( i = 0, tvl = dvl, v = VR(p);
121: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
122: if ( !tvl ) {
123: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
124: ptod(vl,dvl,COEF(dc),&t); pwrp(vl,x,DEG(dc),&c);
125: muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
126: }
127: *pr = s;
128: } else {
129: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
130: ptod(vl,dvl,COEF(dc),&t);
131: NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
132: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
133: muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
134: }
135: *pr = s;
136: }
137: }
138: }
139: if ( !dp_fcoeffs && has_fcoef(*pr) )
140: dp_fcoeffs = 1;
141: }
142:
143: void dtop(vl,dvl,p,pr)
144: VL vl,dvl;
145: DP p;
146: P *pr;
147: {
148: int n,i;
149: DL d;
150: MP m;
151: P r,s,t,u,w;
152: Q q;
153: VL tvl;
154:
155: if ( !p )
156: *pr = 0;
157: else {
158: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
159: t = C(m);
160: if ( NUM(t) && NID((Num)t) == N_M ) {
161: mptop(t,&u); t = u;
162: }
163: for ( i = 0, d = m->dl, tvl = dvl;
164: i < n; tvl = NEXT(tvl), i++ ) {
165: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
166: mulp(vl,t,u,&w); t = w;
167: }
168: addp(vl,s,t,&u); s = u;
169: }
170: *pr = s;
171: }
172: }
173:
174: void nodetod(node,dp)
175: NODE node;
176: DP *dp;
177: {
178: NODE t;
179: int len,i,td;
180: Q e;
181: DL d;
182: MP m;
183: DP u;
184:
185: for ( t = node, len = 0; t; t = NEXT(t), len++ );
186: NEWDL(d,len);
187: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
188: e = (Q)BDY(t);
189: if ( !e )
190: d->d[i] = 0;
191: else if ( !NUM(e) || !RATN(e) || !INT(e) )
192: error("nodetod : invalid input");
193: else {
194: d->d[i] = QTOS((Q)e); td += d->d[i];
195: }
196: }
197: d->td = td;
198: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0;
199: MKDP(len,m,u); u->sugar = td; *dp = u;
200: }
201:
202: int sugard(m)
203: MP m;
204: {
205: int s;
206:
207: for ( s = 0; m; m = NEXT(m) )
208: s = MAX(s,m->dl->td);
209: return s;
210: }
211:
212: void addd(vl,p1,p2,pr)
213: VL vl;
214: DP p1,p2,*pr;
215: {
216: int n;
217: MP m1,m2,mr,mr0;
218: P t;
219:
220: if ( !p1 )
221: *pr = p2;
222: else if ( !p2 )
223: *pr = p1;
224: else {
225: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
226: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
227: case 0:
228: addp(vl,C(m1),C(m2),&t);
229: if ( t ) {
230: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
231: }
232: m1 = NEXT(m1); m2 = NEXT(m2); break;
233: case 1:
234: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
235: m1 = NEXT(m1); break;
236: case -1:
237: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
238: m2 = NEXT(m2); break;
239: }
240: if ( !mr0 )
241: if ( m1 )
242: mr0 = m1;
243: else if ( m2 )
244: mr0 = m2;
245: else {
246: *pr = 0;
247: return;
248: }
249: else if ( m1 )
250: NEXT(mr) = m1;
251: else if ( m2 )
252: NEXT(mr) = m2;
253: else
254: NEXT(mr) = 0;
255: MKDP(NV(p1),mr0,*pr);
256: if ( *pr )
257: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
258: }
259: }
260:
261: /* for F4 symbolic reduction */
262:
263: void symb_addd(p1,p2,pr)
264: DP p1,p2,*pr;
265: {
266: int n;
267: MP m1,m2,mr,mr0;
268: P t;
269:
270: if ( !p1 )
271: *pr = p2;
272: else if ( !p2 )
273: *pr = p1;
274: else {
275: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
276: NEXTMP(mr0,mr); C(mr) = (P)ONE;
277: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
278: case 0:
279: mr->dl = m1->dl;
280: m1 = NEXT(m1); m2 = NEXT(m2); break;
281: case 1:
282: mr->dl = m1->dl;
283: m1 = NEXT(m1); break;
284: case -1:
285: mr->dl = m2->dl;
286: m2 = NEXT(m2); break;
287: }
288: }
289: if ( !mr0 )
290: if ( m1 )
291: mr0 = m1;
292: else if ( m2 )
293: mr0 = m2;
294: else {
295: *pr = 0;
296: return;
297: }
298: else if ( m1 )
299: NEXT(mr) = m1;
300: else if ( m2 )
301: NEXT(mr) = m2;
302: else
303: NEXT(mr) = 0;
304: MKDP(NV(p1),mr0,*pr);
305: if ( *pr )
306: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
307: }
308: }
309:
310: void subd(vl,p1,p2,pr)
311: VL vl;
312: DP p1,p2,*pr;
313: {
314: DP t;
315:
316: if ( !p2 )
317: *pr = p1;
318: else {
319: chsgnd(p2,&t); addd(vl,p1,t,pr);
320: }
321: }
322:
323: void chsgnd(p,pr)
324: DP p,*pr;
325: {
326: MP m,mr,mr0;
327:
328: if ( !p )
329: *pr = 0;
330: else {
331: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
332: NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
333: }
334: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
335: if ( *pr )
336: (*pr)->sugar = p->sugar;
337: }
338: }
339:
340: void muld(vl,p1,p2,pr)
341: VL vl;
342: DP p1,p2,*pr;
343: {
344: MP m;
345: DP s,t,u;
346:
347: if ( !p1 || !p2 )
348: *pr = 0;
349: else if ( OID(p1) <= O_P )
350: muldc(vl,p2,(P)p1,pr);
351: else if ( OID(p2) <= O_P )
352: muldc(vl,p1,(P)p2,pr);
353: else {
354: for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) {
355: muldm(vl,p1,m,&t); addd(vl,s,t,&u); s = u;
356: }
357: *pr = s;
358: }
359: }
360:
361: void muldm(vl,p,m0,pr)
362: VL vl;
363: DP p;
364: MP m0;
365: DP *pr;
366: {
367: MP m,mr,mr0;
368: P c;
369: DL d;
370: int n;
371:
372: if ( !p )
373: *pr = 0;
374: else {
375: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
376: m; m = NEXT(m) ) {
377: NEXTMP(mr0,mr);
378: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
379: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
380: else
381: mulp(vl,C(m),c,&C(mr));
382: adddl(n,m->dl,d,&mr->dl);
383: }
384: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
385: if ( *pr )
386: (*pr)->sugar = p->sugar + m0->dl->td;
387: }
388: }
389:
390: void muldc(vl,p,c,pr)
391: VL vl;
392: DP p;
393: P c;
394: DP *pr;
395: {
396: MP m,mr,mr0;
397:
398: if ( !p || !c )
399: *pr = 0;
400: else if ( NUM(c) && UNIQ((Q)c) )
401: *pr = p;
402: else if ( NUM(c) && MUNIQ((Q)c) )
403: chsgnd(p,pr);
404: else {
405: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
406: NEXTMP(mr0,mr);
407: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
408: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
409: else
410: mulp(vl,C(m),c,&C(mr));
411: mr->dl = m->dl;
412: }
413: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
414: if ( *pr )
415: (*pr)->sugar = p->sugar;
416: }
417: }
418:
419: void divsdc(vl,p,c,pr)
420: VL vl;
421: DP p;
422: P c;
423: DP *pr;
424: {
425: MP m,mr,mr0;
426:
427: if ( !c )
428: error("disvsdc : division by 0");
429: else if ( !p )
430: *pr = 0;
431: else {
432: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
433: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
434: }
435: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
436: if ( *pr )
437: (*pr)->sugar = p->sugar;
438: }
439: }
440:
441: void adddl(n,d1,d2,dr)
442: int n;
443: DL d1,d2;
444: DL *dr;
445: {
446: DL dt;
447: int i;
448:
449: if ( !d1->td )
450: *dr = d2;
451: else if ( !d2->td )
452: *dr = d1;
453: else {
454: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
455: dt->td = d1->td + d2->td;
456: for ( i = 0; i < n; i++ )
457: dt->d[i] = d1->d[i]+d2->d[i];
458: }
459: }
460:
461: int compd(vl,p1,p2)
462: VL vl;
463: DP p1,p2;
464: {
465: int n,t;
466: MP m1,m2;
467:
468: if ( !p1 )
469: return p2 ? -1 : 0;
470: else if ( !p2 )
471: return 1;
472: else {
473: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
474: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
475: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
476: (t = compp(vl,C(m1),C(m2)) ) )
477: return t;
478: if ( m1 )
479: return 1;
480: else if ( m2 )
481: return -1;
482: else
483: return 0;
484: }
485: }
486:
487: int cmpdl_lex(n,d1,d2)
488: int n;
489: DL d1,d2;
490: {
491: int i;
492:
493: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
494: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
495: }
496:
497: int cmpdl_revlex(n,d1,d2)
498: int n;
499: DL d1,d2;
500: {
501: int i;
502:
503: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
504: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
505: }
506:
507: int cmpdl_gradlex(n,d1,d2)
508: int n;
509: DL d1,d2;
510: {
511: if ( d1->td > d2->td )
512: return 1;
513: else if ( d1->td < d2->td )
514: return -1;
515: else
516: return cmpdl_lex(n,d1,d2);
517: }
518:
519: int cmpdl_revgradlex(n,d1,d2)
520: int n;
521: DL d1,d2;
522: {
523: if ( d1->td > d2->td )
524: return 1;
525: else if ( d1->td < d2->td )
526: return -1;
527: else
528: return cmpdl_revlex(n,d1,d2);
529: }
530:
531: int cmpdl_blex(n,d1,d2)
532: int n;
533: DL d1,d2;
534: {
535: int c;
536:
537: if ( c = cmpdl_lex(n-1,d1,d2) )
538: return c;
539: else {
540: c = d1->d[n-1] - d2->d[n-1];
541: return c > 0 ? 1 : c < 0 ? -1 : 0;
542: }
543: }
544:
545: int cmpdl_bgradlex(n,d1,d2)
546: int n;
547: DL d1,d2;
548: {
549: int e1,e2,c;
550:
551: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
552: if ( e1 > e2 )
553: return 1;
554: else if ( e1 < e2 )
555: return -1;
556: else {
557: c = cmpdl_lex(n-1,d1,d2);
558: if ( c )
559: return c;
560: else
561: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
562: }
563: }
564:
565: int cmpdl_brevgradlex(n,d1,d2)
566: int n;
567: DL d1,d2;
568: {
569: int e1,e2,c;
570:
571: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
572: if ( e1 > e2 )
573: return 1;
574: else if ( e1 < e2 )
575: return -1;
576: else {
577: c = cmpdl_revlex(n-1,d1,d2);
578: if ( c )
579: return c;
580: else
581: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
582: }
583: }
584:
585: int cmpdl_brevrev(n,d1,d2)
586: int n;
587: DL d1,d2;
588: {
589: int e1,e2,f1,f2,c,i;
590:
591: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
592: e1 += d1->d[i]; e2 += d2->d[i];
593: }
594: f1 = d1->td - e1; f2 = d2->td - e2;
595: if ( e1 > e2 )
596: return 1;
597: else if ( e1 < e2 )
598: return -1;
599: else {
600: c = cmpdl_revlex(dp_nelim,d1,d2);
601: if ( c )
602: return c;
603: else if ( f1 > f2 )
604: return 1;
605: else if ( f1 < f2 )
606: return -1;
607: else {
608: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
609: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
610: }
611: }
612: }
613:
614: int cmpdl_bgradrev(n,d1,d2)
615: int n;
616: DL d1,d2;
617: {
618: int e1,e2,f1,f2,c,i;
619:
620: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
621: e1 += d1->d[i]; e2 += d2->d[i];
622: }
623: f1 = d1->td - e1; f2 = d2->td - e2;
624: if ( e1 > e2 )
625: return 1;
626: else if ( e1 < e2 )
627: return -1;
628: else {
629: c = cmpdl_lex(dp_nelim,d1,d2);
630: if ( c )
631: return c;
632: else if ( f1 > f2 )
633: return 1;
634: else if ( f1 < f2 )
635: return -1;
636: else {
637: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
638: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
639: }
640: }
641: }
642:
643: int cmpdl_blexrev(n,d1,d2)
644: int n;
645: DL d1,d2;
646: {
647: int e1,e2,f1,f2,c,i;
648:
649: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
650: e1 += d1->d[i]; e2 += d2->d[i];
651: }
652: f1 = d1->td - e1; f2 = d2->td - e2;
653: c = cmpdl_lex(dp_nelim,d1,d2);
654: if ( c )
655: return c;
656: else if ( f1 > f2 )
657: return 1;
658: else if ( f1 < f2 )
659: return -1;
660: else {
661: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
662: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
663: }
664: }
665:
666: int cmpdl_elim(n,d1,d2)
667: int n;
668: DL d1,d2;
669: {
670: int e1,e2,i;
671:
672: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
673: e1 += d1->d[i]; e2 += d2->d[i];
674: }
675: if ( e1 > e2 )
676: return 1;
677: else if ( e1 < e2 )
678: return -1;
679: else
680: return cmpdl_revgradlex(n,d1,d2);
681: }
682:
683: int cmpdl_order_pair(n,d1,d2)
684: int n;
685: DL d1,d2;
686: {
687: int e1,e2,i,j,l;
688: int *t1,*t2;
689: int len;
690: struct order_pair *pair;
691:
692: len = dp_current_spec.ord.block.length;
693: pair = dp_current_spec.ord.block.order_pair;
694:
695: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
696: l = pair[i].length;
697: switch ( pair[i].order ) {
698: case 0:
699: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
700: e1 += t1[j]; e2 += t2[j];
701: }
702: if ( e1 > e2 )
703: return 1;
704: else if ( e1 < e2 )
705: return -1;
706: else {
707: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
708: if ( j >= 0 )
709: return t1[j] < t2[j] ? 1 : -1;
710: }
711: break;
712: case 1:
713: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
714: e1 += t1[j]; e2 += t2[j];
715: }
716: if ( e1 > e2 )
717: return 1;
718: else if ( e1 < e2 )
719: return -1;
720: else {
721: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
722: if ( j < l )
723: return t1[j] > t2[j] ? 1 : -1;
724: }
725: break;
726: case 2:
727: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
728: if ( j < l )
729: return t1[j] > t2[j] ? 1 : -1;
730: break;
731: default:
732: error("cmpdl_order_pair : invalid order"); break;
733: }
734: t1 += l; t2 += l;
735: }
736: return 0;
737: }
738:
739: int cmpdl_matrix(n,d1,d2)
740: int n;
741: DL d1,d2;
742: {
743: int *v,*w,*t1,*t2;
744: int s,i,j,len;
745: int **matrix;
746:
747: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
748: w[i] = t1[i]-t2[i];
749: len = dp_current_spec.ord.matrix.row;
750: matrix = dp_current_spec.ord.matrix.matrix;
751: for ( j = 0; j < len; j++ ) {
752: v = matrix[j];
753: for ( i = 0, s = 0; i < n; i++ )
754: s += v[i]*w[i];
755: if ( s > 0 )
756: return 1;
757: else if ( s < 0 )
758: return -1;
759: }
760: return 0;
761: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>