![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to distm.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 distm.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/distm.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include "inline.h"
4:
5: #define NV(p) ((p)->nv)
6: #define C(p) ((p)->c)
7: #if 0
8: #define ITOS(p) (((unsigned int)(p))>>1)
9: #define STOI(i) ((P)((((unsigned int)(i))<<1)|1))
10: #else
11: #define ITOS(p) (((unsigned int)(p))&0x7fffffff)
12: #define STOI(i) ((P)((unsigned int)(i)|0x80000000))
13: #endif
14:
15: extern int (*cmpdl)();
16:
17: void ptomd(vl,mod,dvl,p,pr)
18: VL vl,dvl;
19: int mod;
20: P p;
21: DP *pr;
22: {
23: P t;
24:
25: ptomp(mod,p,&t);
26: mptomd(vl,mod,dvl,t,pr);
27: }
28:
29: void mptomd(vl,mod,dvl,p,pr)
30: VL vl,dvl;
31: int mod;
32: P p;
33: DP *pr;
34: {
35: int n,i;
36: VL tvl;
37: V v;
38: DL d;
39: MP m;
40: DCP dc;
41: DP r,s,t,u;
42: P x,c;
43:
44: if ( !p )
45: *pr = 0;
46: else {
47: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
48: if ( NUM(p) ) {
49: NEWDL(d,n);
50: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr);
51: } else {
52: for ( i = 0, tvl = dvl, v = VR(p);
53: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
54: if ( !tvl ) {
55: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
56: mptomd(vl,mod,dvl,COEF(dc),&t); pwrmp(vl,mod,x,DEG(dc),&c);
57: mulmdc(vl,mod,t,c,&r); addmd(vl,mod,r,s,&t); s = t;
58: }
59: *pr = s;
60: } else {
61: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
62: mptomd(vl,mod,dvl,COEF(dc),&t);
63: NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
64: NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u);
65: mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t;
66: }
67: *pr = s;
68: }
69: }
70: }
71: }
72:
73: void mdtop(vl,mod,dvl,p,pr)
74: VL vl,dvl;
75: int mod;
76: DP p;
77: P *pr;
78: {
79: int n,i;
80: DL d;
81: MP m;
82: P r,s,t,u,w;
83: Q q;
84: VL tvl;
85:
86: if ( !p )
87: *pr = 0;
88: else {
89: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
90: for ( i = 0, t = C(m), d = m->dl, tvl = dvl;
91: i < n; tvl = NEXT(tvl), i++ ) {
92: MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
93: mulmp(vl,mod,t,u,&w); t = w;
94: }
95: addmp(vl,mod,s,t,&u); s = u;
96: }
97: mptop(s,pr);
98: }
99: }
100:
101: void addmd(vl,mod,p1,p2,pr)
102: VL vl;
103: int mod;
104: DP p1,p2,*pr;
105: {
106: int n;
107: MP m1,m2,mr,mr0;
108: P t;
109: int tmp;
110: MQ q;
111:
112: if ( !p1 )
113: *pr = p2;
114: else if ( !p2 )
115: *pr = p1;
116: else {
117: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
118: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
119: case 0:
120: if ( NUM(C(m1)) && NUM(C(m2)) ) {
121: tmp = (CONT((MQ)C(m1))+CONT((MQ)C(m2))) - mod;
122: if ( tmp < 0 )
123: tmp += mod;
124: if ( tmp ) {
125: STOMQ(tmp,q); t = (P)q;
126: } else
127: t = 0;
128: } else
129: addmp(vl,mod,C(m1),C(m2),&t);
130: if ( t ) {
131: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
132: }
133: m1 = NEXT(m1); m2 = NEXT(m2); break;
134: case 1:
135: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
136: m1 = NEXT(m1); break;
137: case -1:
138: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
139: m2 = NEXT(m2); break;
140: }
141: if ( !mr0 )
142: if ( m1 )
143: mr0 = m1;
144: else if ( m2 )
145: mr0 = m2;
146: else {
147: *pr = 0;
148: return;
149: }
150: else if ( m1 )
151: NEXT(mr) = m1;
152: else if ( m2 )
153: NEXT(mr) = m2;
154: else
155: NEXT(mr) = 0;
156: MKDP(NV(p1),mr0,*pr);
157: if ( *pr )
158: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
159: }
160: }
161:
162: void submd(vl,mod,p1,p2,pr)
163: VL vl;
164: int mod;
165: DP p1,p2,*pr;
166: {
167: DP t;
168:
169: if ( !p2 )
170: *pr = p1;
171: else {
172: chsgnmd(mod,p2,&t); addmd(vl,mod,p1,t,pr);
173: }
174: }
175:
176: void chsgnmd(mod,p,pr)
177: int mod;
178: DP p,*pr;
179: {
180: MP m,mr,mr0;
181:
182: if ( !p )
183: *pr = 0;
184: else {
185: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
186: NEXTMP(mr0,mr); chsgnmp(mod,C(m),&C(mr)); mr->dl = m->dl;
187: }
188: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
189: if ( *pr )
190: (*pr)->sugar = p->sugar;
191: }
192: }
193:
194: void mulmd(vl,mod,p1,p2,pr)
195: VL vl;
196: int mod;
197: DP p1,p2,*pr;
198: {
199: MP m;
200: DP s,t,u;
201:
202: if ( !p1 || !p2 )
203: *pr = 0;
204: else if ( OID(p1) <= O_P )
205: mulmdc(vl,mod,p2,(P)p1,pr);
206: else if ( OID(p2) <= O_P )
207: mulmdc(vl,mod,p1,(P)p2,pr);
208: else {
209: for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) {
210: mulmdm(vl,mod,p1,m,&t); addmd(vl,mod,s,t,&u); s = u;
211: }
212: *pr = s;
213: }
214: }
215:
216: void mulmdm(vl,mod,p,m0,pr)
217: VL vl;
218: int mod;
219: DP p;
220: MP m0;
221: DP *pr;
222: {
223: MP m,mr,mr0;
224: P c;
225: DL d;
226: int n,t;
227: MQ q;
228:
229: if ( !p )
230: *pr = 0;
231: else {
232: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
233: m; m = NEXT(m) ) {
234: NEXTMP(mr0,mr);
235: if ( NUM(C(m)) && NUM(c) ) {
236: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
237: STOMQ(t,q); C(mr) = (P)q;
238: } else
239: mulmp(vl,mod,C(m),c,&C(mr));
240: adddl(n,m->dl,d,&mr->dl);
241: }
242: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
243: if ( *pr )
244: (*pr)->sugar = p->sugar + m0->dl->td;
245: }
246: }
247:
248: void mulmdc(vl,mod,p,c,pr)
249: VL vl;
250: int mod;
251: DP p;
252: P c;
253: DP *pr;
254: {
255: MP m,mr,mr0;
256: int t;
257: MQ q;
258:
259: if ( !p || !c )
260: *pr = 0;
261: else {
262: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
263: NEXTMP(mr0,mr);
264: if ( NUM(C(m)) && NUM(c) ) {
265: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
266: STOMQ(t,q); C(mr) = (P)q;
267: } else
268: mulmp(vl,mod,C(m),c,&C(mr));
269: mr->dl = m->dl;
270: }
271: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
272: if ( *pr )
273: (*pr)->sugar = p->sugar;
274: }
275: }
276:
277: void divsmdc(vl,mod,p,c,pr)
278: VL vl;
279: int mod;
280: DP p;
281: P c;
282: DP *pr;
283: {
284: MP m,mr,mr0;
285:
286: if ( !c )
287: error("disvsdc : division by 0");
288: else if ( !p )
289: *pr = 0;
290: else {
291: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
292: NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl;
293: }
294: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
295: if ( *pr )
296: (*pr)->sugar = p->sugar;
297: }
298: }
299:
300: #define MKDPM(n,m,d) (NEWDP(d),(d)->nv=(n),BDY(d)=(m))
301:
302: void _mdtop(vl,mod,dvl,p,pr)
303: VL vl,dvl;
304: int mod;
305: DP p;
306: P *pr;
307: {
308: int n,i;
309: DL d;
310: MP m;
311: P r,s,t,u,w;
312: Q q;
313: MQ tq;
314: VL tvl;
315:
316: if ( !p )
317: *pr = 0;
318: else {
319: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
320: STOMQ(ITOS(C(m)),tq); t = (P)tq;
321: for ( i = 0, d = m->dl, tvl = dvl;
322: i < n; tvl = NEXT(tvl), i++ ) {
323: MKV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
324: mulmp(vl,mod,t,u,&w); t = w;
325: }
326: addmp(vl,mod,s,t,&u); s = u;
327: }
328: mptop(s,pr);
329: }
330: }
331:
332: void _addmd(vl,mod,p1,p2,pr)
333: VL vl;
334: int mod;
335: DP p1,p2,*pr;
336: {
337: int n;
338: MP m1,m2,mr,mr0;
339: int t;
340:
341: if ( !p1 )
342: *pr = p2;
343: else if ( !p2 )
344: *pr = p1;
345: else {
346: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
347: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
348: case 0:
349: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
350: if ( t < 0 )
351: t += mod;
352: if ( t ) {
353: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = STOI(t);
354: }
355: m1 = NEXT(m1); m2 = NEXT(m2); break;
356: case 1:
357: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
358: m1 = NEXT(m1); break;
359: case -1:
360: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
361: m2 = NEXT(m2); break;
362: }
363: if ( !mr0 )
364: if ( m1 )
365: mr0 = m1;
366: else if ( m2 )
367: mr0 = m2;
368: else {
369: *pr = 0;
370: return;
371: }
372: else if ( m1 )
373: NEXT(mr) = m1;
374: else if ( m2 )
375: NEXT(mr) = m2;
376: else
377: NEXT(mr) = 0;
378: MKDPM(NV(p1),mr0,*pr);
379: if ( *pr )
380: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
381: }
382: }
383:
384: void _submd(vl,mod,p1,p2,pr)
385: VL vl;
386: int mod;
387: DP p1,p2,*pr;
388: {
389: DP t;
390:
391: if ( !p2 )
392: *pr = p1;
393: else {
394: _chsgnmd(mod,p2,&t); _addmd(vl,mod,p1,t,pr);
395: }
396: }
397:
398: void _chsgnmd(mod,p,pr)
399: int mod;
400: DP p,*pr;
401: {
402: MP m,mr,mr0;
403:
404: if ( !p )
405: *pr = 0;
406: else {
407: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
408: NEXTMP(mr0,mr); C(mr) = STOI(mod - ITOS(C(m))); mr->dl = m->dl;
409: }
410: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
411: if ( *pr )
412: (*pr)->sugar = p->sugar;
413: }
414: }
415:
416: void _mulmd(vl,mod,p1,p2,pr)
417: VL vl;
418: int mod;
419: DP p1,p2,*pr;
420: {
421: MP m;
422: DP s,t,u;
423:
424: if ( !p1 || !p2 )
425: *pr = 0;
426: else {
427: for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) {
428: _mulmdm(vl,mod,p1,m,&t); _addmd(vl,mod,s,t,&u); s = u;
429: }
430: *pr = s;
431: }
432: }
433:
434: void _mulmdm(vl,mod,p,m0,pr)
435: VL vl;
436: int mod;
437: DP p;
438: MP m0;
439: DP *pr;
440: {
441: MP m,mr,mr0;
442: DL d;
443: int c,n,r;
444:
445: if ( !p )
446: *pr = 0;
447: else {
448: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
449: m; m = NEXT(m) ) {
450: NEXTMP(mr0,mr);
451: C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
452: adddl(n,m->dl,d,&mr->dl);
453: }
454: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
455: if ( *pr )
456: (*pr)->sugar = p->sugar + m0->dl->td;
457: }
458: }
459:
460: void _dtop_mod(vl,dvl,p,pr)
461: VL vl,dvl;
462: DP p;
463: P *pr;
464: {
465: int n,i;
466: DL d;
467: MP m;
468: P r,s,t,u,w;
469: Q q;
470: VL tvl;
471:
472: if ( !p )
473: *pr = 0;
474: else {
475: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
476: i = ITOS(m->c); STOQ(i,q); t = (P)q;
477: for ( i = 0, d = m->dl, tvl = dvl;
478: i < n; tvl = NEXT(tvl), i++ ) {
479: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
480: mulp(vl,t,u,&w); t = w;
481: }
482: addp(vl,s,t,&u); s = u;
483: }
484: *pr = s;
485: }
486: }
487:
488: void _dp_red_mod(p1,p2,mod,rp)
489: DP p1,p2;
490: int mod;
491: DP *rp;
492: {
493: int i,n;
494: DL d1,d2,d;
495: MP m;
496: DP t,s;
497: int c,c1;
498:
499: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
500: NEWDL(d,n); d->td = d1->td - d2->td;
501: for ( i = 0; i < n; i++ )
502: d->d[i] = d1->d[i]-d2->d[i];
503: c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
504: NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
505: MKDP(n,m,s); s->sugar = d->td;
506: _mulmd(CO,mod,p2,s,&t); _addmd(CO,mod,p1,t,rp);
507: }
508:
509: void _dp_mod(p,mod,subst,rp)
510: DP p;
511: int mod;
512: NODE subst;
513: DP *rp;
514: {
515: MP m,mr,mr0;
516: P t,s;
517: NODE tn;
518:
519: if ( !p )
520: *rp = 0;
521: else {
522: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
523: for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) {
524: substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t);
525: s = t;
526: }
527: ptomp(mod,s,&t);
528: if ( t ) {
529: NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl;
530: }
531: }
532: if ( mr0 ) {
533: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
534: } else
535: *rp = 0;
536: }
537: }
538:
539: void _dp_sp_mod(p1,p2,mod,rp)
540: DP p1,p2;
541: int mod;
542: DP *rp;
543: {
544: int i,n,td;
545: int *w;
546: DL d1,d2,d;
547: MP m;
548: DP t,s,u;
549:
550: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
551: w = (int *)ALLOCA(n*sizeof(int));
552: for ( i = 0, td = 0; i < n; i++ ) {
553: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
554: }
555: NEWDL(d,n); d->td = td - d1->td;
556: for ( i = 0; i < n; i++ )
557: d->d[i] = w[i] - d1->d[i];
558: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
559: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p1,s,&t);
560: NEWDL(d,n); d->td = td - d2->td;
561: for ( i = 0; i < n; i++ )
562: d->d[i] = w[i] - d2->d[i];
563: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
564: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p2,s,&u);
565: _addmd(CO,mod,t,u,rp);
566: }
567:
568: void _dp_sp_component_mod(p1,p2,mod,f1,f2)
569: DP p1,p2;
570: int mod;
571: DP *f1,*f2;
572: {
573: int i,n,td;
574: int *w;
575: DL d1,d2,d;
576: MP m;
577: DP t,s,u;
578:
579: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
580: w = (int *)ALLOCA(n*sizeof(int));
581: for ( i = 0, td = 0; i < n; i++ ) {
582: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
583: }
584: NEWDL(d,n); d->td = td - d1->td;
585: for ( i = 0; i < n; i++ )
586: d->d[i] = w[i] - d1->d[i];
587: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
588: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p1,s,f1);
589: NEWDL(d,n); d->td = td - d2->td;
590: for ( i = 0; i < n; i++ )
591: d->d[i] = w[i] - d2->d[i];
592: NEWMP(m); m->dl = d; m->c = BDY(p1)->c; NEXT(m) = 0;
593: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p2,s,f2);
594: }
595:
596: void _printdp(d)
597: DP d;
598: {
599: int n,i;
600: MP m;
601: DL dl;
602:
603: if ( !d ) {
604: printf("0"); return;
605: }
606: for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) {
607: printf("%d*<<",ITOS(m->c));
608: for ( i = 0, dl = m->dl; i < n-1; i++ )
609: printf("%d,",dl->d[i]);
610: printf("%d",dl->d[i]);
611: printf(">>");
612: if ( NEXT(m) )
613: printf("+");
614: }
615: }