![[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.4 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.3 2000/05/30 01:35:12 noro Exp $ */
1.1 noro 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)();
1.2 noro 16: extern int do_weyl;
17:
18: void comm_mulmd();
19: void weyl_mulmd();
20: void weyl_mulmdm();
21: void weyl_mulmmm();
22: void _comm_mulmd();
23: void _weyl_mulmd();
24: void _weyl_mulmmm();
25: void _weyl_mulmdm();
1.1 noro 26:
27: void ptomd(vl,mod,dvl,p,pr)
28: VL vl,dvl;
29: int mod;
30: P p;
31: DP *pr;
32: {
33: P t;
34:
35: ptomp(mod,p,&t);
36: mptomd(vl,mod,dvl,t,pr);
37: }
38:
39: void mptomd(vl,mod,dvl,p,pr)
40: VL vl,dvl;
41: int mod;
42: P p;
43: DP *pr;
44: {
45: int n,i;
46: VL tvl;
47: V v;
48: DL d;
49: MP m;
50: DCP dc;
51: DP r,s,t,u;
52: P x,c;
53:
54: if ( !p )
55: *pr = 0;
56: else {
57: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
58: if ( NUM(p) ) {
59: NEWDL(d,n);
60: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr);
61: } else {
62: for ( i = 0, tvl = dvl, v = VR(p);
63: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
64: if ( !tvl ) {
65: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
66: mptomd(vl,mod,dvl,COEF(dc),&t); pwrmp(vl,mod,x,DEG(dc),&c);
67: mulmdc(vl,mod,t,c,&r); addmd(vl,mod,r,s,&t); s = t;
68: }
69: *pr = s;
70: } else {
71: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
72: mptomd(vl,mod,dvl,COEF(dc),&t);
73: NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
74: NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u);
1.2 noro 75: comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t;
1.1 noro 76: }
77: *pr = s;
78: }
79: }
80: }
81: }
82:
83: void mdtop(vl,mod,dvl,p,pr)
84: VL vl,dvl;
85: int mod;
86: DP p;
87: P *pr;
88: {
89: int n,i;
90: DL d;
91: MP m;
92: P r,s,t,u,w;
93: Q q;
94: VL tvl;
95:
96: if ( !p )
97: *pr = 0;
98: else {
99: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
100: for ( i = 0, t = C(m), d = m->dl, tvl = dvl;
101: i < n; tvl = NEXT(tvl), i++ ) {
102: MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
103: mulmp(vl,mod,t,u,&w); t = w;
104: }
105: addmp(vl,mod,s,t,&u); s = u;
106: }
107: mptop(s,pr);
108: }
109: }
110:
111: void addmd(vl,mod,p1,p2,pr)
112: VL vl;
113: int mod;
114: DP p1,p2,*pr;
115: {
116: int n;
117: MP m1,m2,mr,mr0;
118: P t;
119: int tmp;
120: MQ q;
121:
122: if ( !p1 )
123: *pr = p2;
124: else if ( !p2 )
125: *pr = p1;
126: else {
127: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
128: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
129: case 0:
130: if ( NUM(C(m1)) && NUM(C(m2)) ) {
131: tmp = (CONT((MQ)C(m1))+CONT((MQ)C(m2))) - mod;
132: if ( tmp < 0 )
133: tmp += mod;
134: if ( tmp ) {
135: STOMQ(tmp,q); t = (P)q;
136: } else
137: t = 0;
138: } else
139: addmp(vl,mod,C(m1),C(m2),&t);
140: if ( t ) {
141: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
142: }
143: m1 = NEXT(m1); m2 = NEXT(m2); break;
144: case 1:
145: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
146: m1 = NEXT(m1); break;
147: case -1:
148: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
149: m2 = NEXT(m2); break;
150: }
151: if ( !mr0 )
152: if ( m1 )
153: mr0 = m1;
154: else if ( m2 )
155: mr0 = m2;
156: else {
157: *pr = 0;
158: return;
159: }
160: else if ( m1 )
161: NEXT(mr) = m1;
162: else if ( m2 )
163: NEXT(mr) = m2;
164: else
165: NEXT(mr) = 0;
166: MKDP(NV(p1),mr0,*pr);
167: if ( *pr )
168: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
169: }
170: }
171:
172: void submd(vl,mod,p1,p2,pr)
173: VL vl;
174: int mod;
175: DP p1,p2,*pr;
176: {
177: DP t;
178:
179: if ( !p2 )
180: *pr = p1;
181: else {
182: chsgnmd(mod,p2,&t); addmd(vl,mod,p1,t,pr);
183: }
184: }
185:
186: void chsgnmd(mod,p,pr)
187: int mod;
188: DP p,*pr;
189: {
190: MP m,mr,mr0;
191:
192: if ( !p )
193: *pr = 0;
194: else {
195: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
196: NEXTMP(mr0,mr); chsgnmp(mod,C(m),&C(mr)); mr->dl = m->dl;
197: }
198: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
199: if ( *pr )
200: (*pr)->sugar = p->sugar;
201: }
202: }
203:
204: void mulmd(vl,mod,p1,p2,pr)
205: VL vl;
206: int mod;
207: DP p1,p2,*pr;
208: {
1.2 noro 209: if ( !do_weyl )
210: comm_mulmd(vl,mod,p1,p2,pr);
211: else
212: weyl_mulmd(vl,mod,p1,p2,pr);
213: }
214:
215: void comm_mulmd(vl,mod,p1,p2,pr)
216: VL vl;
217: int mod;
218: DP p1,p2,*pr;
219: {
220: MP m;
221: DP s,t,u;
222: int i,l,l1;
223: static MP *w;
224: static int wlen;
225:
226: if ( !p1 || !p2 )
227: *pr = 0;
228: else if ( OID(p1) <= O_P )
229: mulmdc(vl,mod,p2,(P)p1,pr);
230: else if ( OID(p2) <= O_P )
231: mulmdc(vl,mod,p1,(P)p2,pr);
232: else {
233: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
234: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
235: if ( l1 < l ) {
236: t = p1; p1 = p2; p2 = t;
237: l = l1;
238: }
239: if ( l > wlen ) {
240: if ( w ) GC_free(w);
241: w = (MP *)MALLOC(l*sizeof(MP));
242: wlen = l;
243: }
244: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
245: w[i] = m;
246: for ( s = 0, i = l-1; i >= 0; i-- ) {
247: mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
248: }
249: bzero(w,l*sizeof(MP));
250: *pr = s;
251: }
252: }
253:
254: void weyl_mulmd(vl,mod,p1,p2,pr)
255: VL vl;
256: int mod;
257: DP p1,p2,*pr;
258: {
1.1 noro 259: MP m;
260: DP s,t,u;
1.2 noro 261: int i,l,l1;
262: static MP *w;
263: static int wlen;
1.1 noro 264:
265: if ( !p1 || !p2 )
266: *pr = 0;
267: else if ( OID(p1) <= O_P )
268: mulmdc(vl,mod,p2,(P)p1,pr);
269: else if ( OID(p2) <= O_P )
270: mulmdc(vl,mod,p1,(P)p2,pr);
271: else {
1.2 noro 272: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
273: if ( l > wlen ) {
274: if ( w ) GC_free(w);
275: w = (MP *)MALLOC(l*sizeof(MP));
276: wlen = l;
277: }
278: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
279: w[i] = m;
280: for ( s = 0, i = l-1; i >= 0; i-- ) {
281: weyl_mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
1.1 noro 282: }
1.2 noro 283: bzero(w,l*sizeof(MP));
1.1 noro 284: *pr = s;
285: }
286: }
287:
288: void mulmdm(vl,mod,p,m0,pr)
289: VL vl;
290: int mod;
291: DP p;
292: MP m0;
293: DP *pr;
294: {
295: MP m,mr,mr0;
296: P c;
297: DL d;
298: int n,t;
299: MQ q;
300:
301: if ( !p )
302: *pr = 0;
303: else {
304: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
305: m; m = NEXT(m) ) {
306: NEXTMP(mr0,mr);
307: if ( NUM(C(m)) && NUM(c) ) {
308: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
309: STOMQ(t,q); C(mr) = (P)q;
310: } else
311: mulmp(vl,mod,C(m),c,&C(mr));
312: adddl(n,m->dl,d,&mr->dl);
313: }
314: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
315: if ( *pr )
316: (*pr)->sugar = p->sugar + m0->dl->td;
317: }
318: }
319:
1.2 noro 320: void weyl_mulmdm(vl,mod,p,m0,pr)
321: VL vl;
322: int mod;
323: DP p;
324: MP m0;
325: DP *pr;
326: {
327: DP r,t,t1;
328: MP m;
329: int n,l,i;
330: static MP *w;
331: static int wlen;
332:
333: if ( !p )
334: *pr = 0;
335: else {
336: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
337: if ( l > wlen ) {
338: if ( w ) GC_free(w);
339: w = (MP *)MALLOC(l*sizeof(MP));
340: wlen = l;
341: }
342: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
343: w[i] = m;
344: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
345: weyl_mulmmm(vl,mod,w[i],m0,n,&t);
346: addmd(vl,mod,r,t,&t1); r = t1;
347: }
348: bzero(w,l*sizeof(MP));
349: if ( r )
350: r->sugar = p->sugar + m0->dl->td;
351: *pr = r;
352: }
353: }
354:
355: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
356:
357: void weyl_mulmmm(vl,mod,m0,m1,n,pr)
358: VL vl;
359: int mod;
360: MP m0,m1;
361: int n;
362: DP *pr;
363: {
364: MP m,mr,mr0;
365: MQ mq;
366: DP r,t,t1;
367: P c,c0,c1,cc;
368: DL d,d0,d1;
1.3 noro 369: int i,j,a,b,k,l,n2,s,min;
1.2 noro 370: static int *tab;
371: static int tablen;
372:
373: if ( !m0 || !m1 )
374: *pr = 0;
375: else {
376: c0 = C(m0); c1 = C(m1);
377: mulmp(vl,mod,c0,c1,&c);
378: d0 = m0->dl; d1 = m1->dl;
379: n2 = n>>1;
380: if ( n & 1 ) {
381: /* homogenized computation; dx-xd=h^2 */
382: /* offset of h-degree */
1.3 noro 383: NEWDL(d,n);
384: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1.4 ! noro 385: NEWMP(mr); mr->c = (P)ONEM; mr->dl = d; NEXT(mr) = 0;
1.3 noro 386: MKDP(n,mr,r); r->sugar = d->td;
387: } else
388: r = (DP)ONEM;
389: for ( i = 0; i < n2; i++ ) {
390: a = d0->d[i]; b = d1->d[n2+i];
391: k = d0->d[n2+i]; l = d1->d[i];
392: /* degree of xi^a*(Di^k*xi^l)*Di^b */
393: s = a+k+l+b;
394: /* compute xi^a*(Di^k*xi^l)*Di^b */
395: min = MIN(k,l);
396:
397: if ( min+1 > tablen ) {
398: if ( tab ) GC_free(tab);
399: tab = (int *)MALLOC((min+1)*sizeof(int));
400: tablen = min+1;
401: }
402: mkwcm(k,l,mod,tab);
403: if ( n & 1 )
1.2 noro 404: for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3 noro 405: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 406: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.3 noro 407: d->td = s;
408: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
409: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 410: }
1.3 noro 411: else
1.2 noro 412: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3 noro 413: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 414: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
415: d->td = d->d[i]+d->d[n2+i]; /* XXX */
416: s = MAX(s,d->td); /* XXX */
1.3 noro 417: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 418: }
1.3 noro 419: bzero(tab,(min+1)*sizeof(int));
420: if ( mr0 )
421: NEXT(mr) = 0;
422: MKDP(n,mr0,t);
423: if ( t )
424: t->sugar = s;
425: comm_mulmd(vl,mod,r,t,&t1); r = t1;
426: }
1.2 noro 427: mulmdc(vl,mod,r,c,pr);
428: }
429: }
430:
1.1 noro 431: void mulmdc(vl,mod,p,c,pr)
432: VL vl;
433: int mod;
434: DP p;
435: P c;
436: DP *pr;
437: {
438: MP m,mr,mr0;
439: int t;
440: MQ q;
441:
442: if ( !p || !c )
443: *pr = 0;
444: else {
445: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
446: NEXTMP(mr0,mr);
447: if ( NUM(C(m)) && NUM(c) ) {
448: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
449: STOMQ(t,q); C(mr) = (P)q;
450: } else
451: mulmp(vl,mod,C(m),c,&C(mr));
452: mr->dl = m->dl;
453: }
454: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
455: if ( *pr )
456: (*pr)->sugar = p->sugar;
457: }
458: }
459:
460: void divsmdc(vl,mod,p,c,pr)
461: VL vl;
462: int mod;
463: DP p;
464: P c;
465: DP *pr;
466: {
467: MP m,mr,mr0;
468:
469: if ( !c )
470: error("disvsdc : division by 0");
471: else if ( !p )
472: *pr = 0;
473: else {
474: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
475: NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl;
476: }
477: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
478: if ( *pr )
479: (*pr)->sugar = p->sugar;
480: }
481: }
482:
483: void _mdtop(vl,mod,dvl,p,pr)
484: VL vl,dvl;
485: int mod;
486: DP p;
487: P *pr;
488: {
489: int n,i;
490: DL d;
491: MP m;
492: P r,s,t,u,w;
493: Q q;
494: MQ tq;
495: VL tvl;
496:
497: if ( !p )
498: *pr = 0;
499: else {
500: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
501: STOMQ(ITOS(C(m)),tq); t = (P)tq;
502: for ( i = 0, d = m->dl, tvl = dvl;
503: i < n; tvl = NEXT(tvl), i++ ) {
504: MKV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
505: mulmp(vl,mod,t,u,&w); t = w;
506: }
507: addmp(vl,mod,s,t,&u); s = u;
508: }
509: mptop(s,pr);
510: }
511: }
512:
513: void _addmd(vl,mod,p1,p2,pr)
514: VL vl;
515: int mod;
516: DP p1,p2,*pr;
517: {
518: int n;
519: MP m1,m2,mr,mr0;
520: int t;
521:
522: if ( !p1 )
523: *pr = p2;
524: else if ( !p2 )
525: *pr = p1;
526: else {
527: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
528: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
529: case 0:
530: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
531: if ( t < 0 )
532: t += mod;
533: if ( t ) {
534: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = STOI(t);
535: }
536: m1 = NEXT(m1); m2 = NEXT(m2); break;
537: case 1:
538: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
539: m1 = NEXT(m1); break;
540: case -1:
541: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
542: m2 = NEXT(m2); break;
543: }
544: if ( !mr0 )
545: if ( m1 )
546: mr0 = m1;
547: else if ( m2 )
548: mr0 = m2;
549: else {
550: *pr = 0;
551: return;
552: }
553: else if ( m1 )
554: NEXT(mr) = m1;
555: else if ( m2 )
556: NEXT(mr) = m2;
557: else
558: NEXT(mr) = 0;
1.3 noro 559: MKDP(NV(p1),mr0,*pr);
1.1 noro 560: if ( *pr )
561: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
562: }
563: }
564:
565: void _submd(vl,mod,p1,p2,pr)
566: VL vl;
567: int mod;
568: DP p1,p2,*pr;
569: {
570: DP t;
571:
572: if ( !p2 )
573: *pr = p1;
574: else {
575: _chsgnmd(mod,p2,&t); _addmd(vl,mod,p1,t,pr);
576: }
577: }
578:
579: void _chsgnmd(mod,p,pr)
580: int mod;
581: DP p,*pr;
582: {
583: MP m,mr,mr0;
584:
585: if ( !p )
586: *pr = 0;
587: else {
588: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
589: NEXTMP(mr0,mr); C(mr) = STOI(mod - ITOS(C(m))); mr->dl = m->dl;
590: }
1.3 noro 591: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1.1 noro 592: if ( *pr )
593: (*pr)->sugar = p->sugar;
594: }
595: }
596:
597: void _mulmd(vl,mod,p1,p2,pr)
598: VL vl;
599: int mod;
600: DP p1,p2,*pr;
601: {
1.2 noro 602: if ( !do_weyl )
603: _comm_mulmd(vl,mod,p1,p2,pr);
604: else
605: _weyl_mulmd(vl,mod,p1,p2,pr);
606: }
607:
608: void _comm_mulmd(vl,mod,p1,p2,pr)
609: VL vl;
610: int mod;
611: DP p1,p2,*pr;
612: {
1.1 noro 613: MP m;
614: DP s,t,u;
1.2 noro 615: int i,l,l1;
616: static MP *w;
617: static int wlen;
1.1 noro 618:
619: if ( !p1 || !p2 )
620: *pr = 0;
621: else {
1.2 noro 622: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
623: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
624: if ( l1 < l ) {
625: t = p1; p1 = p2; p2 = t;
626: l = l1;
627: }
628: if ( l > wlen ) {
629: if ( w ) GC_free(w);
630: w = (MP *)MALLOC(l*sizeof(MP));
631: wlen = l;
632: }
633: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
634: w[i] = m;
635: for ( s = 0, i = l-1; i >= 0; i-- ) {
636: _mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
1.1 noro 637: }
1.2 noro 638: bzero(w,l*sizeof(MP));
639: *pr = s;
640: }
641: }
642:
643: void _weyl_mulmd(vl,mod,p1,p2,pr)
644: VL vl;
645: int mod;
646: DP p1,p2,*pr;
647: {
648: MP m;
649: DP s,t,u;
650: int i,l,l1;
651: static MP *w;
652: static int wlen;
653:
654: if ( !p1 || !p2 )
655: *pr = 0;
656: else {
657: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
658: if ( l > wlen ) {
659: if ( w ) GC_free(w);
660: w = (MP *)MALLOC(l*sizeof(MP));
661: wlen = l;
662: }
663: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
664: w[i] = m;
665: for ( s = 0, i = l-1; i >= 0; i-- ) {
666: _weyl_mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
667: }
668: bzero(w,l*sizeof(MP));
1.1 noro 669: *pr = s;
670: }
671: }
672:
673: void _mulmdm(vl,mod,p,m0,pr)
674: VL vl;
675: int mod;
676: DP p;
677: MP m0;
678: DP *pr;
679: {
680: MP m,mr,mr0;
681: DL d;
682: int c,n,r;
683:
684: if ( !p )
685: *pr = 0;
686: else {
687: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
688: m; m = NEXT(m) ) {
689: NEXTMP(mr0,mr);
690: C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
691: adddl(n,m->dl,d,&mr->dl);
692: }
1.3 noro 693: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1.1 noro 694: if ( *pr )
695: (*pr)->sugar = p->sugar + m0->dl->td;
696: }
697: }
698:
1.2 noro 699: void _weyl_mulmdm(vl,mod,p,m0,pr)
700: VL vl;
701: int mod;
702: DP p;
703: MP m0;
704: DP *pr;
705: {
706: DP r,t,t1;
707: MP m;
708: int n,l,i;
709: static MP *w;
710: static int wlen;
711:
712: if ( !p )
713: *pr = 0;
714: else {
715: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
716: if ( l > wlen ) {
717: if ( w ) GC_free(w);
718: w = (MP *)MALLOC(l*sizeof(MP));
719: wlen = l;
720: }
721: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
722: w[i] = m;
723: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
724: _weyl_mulmmm(vl,mod,w[i],m0,n,&t);
725: _addmd(vl,mod,r,t,&t1); r = t1;
726: }
727: bzero(w,l*sizeof(MP));
728: if ( r )
729: r->sugar = p->sugar + m0->dl->td;
730: *pr = r;
731: }
732: }
733:
734: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
735:
736: void _weyl_mulmmm(vl,mod,m0,m1,n,pr)
737: VL vl;
738: int mod;
739: MP m0,m1;
740: int n;
741: DP *pr;
742: {
743: MP m,mr,mr0;
744: DP r,t,t1;
745: int c,c0,c1,cc;
746: DL d,d0,d1;
747: int i,j,a,b,k,l,n2,s,min,h;
748: static int *tab;
749: static int tablen;
750:
751: if ( !m0 || !m1 )
752: *pr = 0;
753: else {
754: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
755: c = dmar(c0,c1,0,mod);
756: d0 = m0->dl; d1 = m1->dl;
757: n2 = n>>1;
758:
1.3 noro 759: NEWDL(d,n);
760: if ( n & 1 )
1.2 noro 761: /* offset of h-degree */
1.3 noro 762: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
763: else
764: d->td = 0;
1.4 ! noro 765: NEWMP(mr); mr->c = STOI(c); mr->dl = d; NEXT(mr) = 0;
1.3 noro 766: MKDP(n,mr,r); r->sugar = d->td;
767:
768: /* homogenized computation; dx-xd=h^2 */
769: for ( i = 0; i < n2; i++ ) {
770: a = d0->d[i]; b = d1->d[n2+i];
771: k = d0->d[n2+i]; l = d1->d[i];
772: /* degree of xi^a*(Di^k*xi^l)*Di^b */
773: s = a+k+l+b;
774: /* compute xi^a*(Di^k*xi^l)*Di^b */
775: min = MIN(k,l);
776:
777: if ( min+1 > tablen ) {
778: if ( tab ) GC_free(tab);
779: tab = (int *)MALLOC((min+1)*sizeof(int));
780: tablen = min+1;
781: }
782: mkwcm(k,l,mod,tab);
783: if ( n & 1 )
1.2 noro 784: for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3 noro 785: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 786: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.3 noro 787: d->td = s;
788: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
789: mr->c = STOI(tab[j]); mr->dl = d;
1.2 noro 790: }
1.3 noro 791: else
1.2 noro 792: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3 noro 793: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 794: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
795: d->td = d->d[i]+d->d[n2+i]; /* XXX */
796: s = MAX(s,d->td); /* XXX */
1.3 noro 797: mr->c = STOI(tab[j]); mr->dl = d;
1.2 noro 798: }
1.3 noro 799: bzero(tab,(min+1)*sizeof(int));
800: if ( mr0 )
801: NEXT(mr) = 0;
802: MKDP(n,mr0,t);
803: if ( t )
804: t->sugar = s;
805: _comm_mulmd(vl,mod,r,t,&t1); r = t1;
1.2 noro 806: }
807: *pr = r;
808: }
809: }
810:
1.1 noro 811: void _dtop_mod(vl,dvl,p,pr)
812: VL vl,dvl;
813: DP p;
814: P *pr;
815: {
816: int n,i;
817: DL d;
818: MP m;
819: P r,s,t,u,w;
820: Q q;
821: VL tvl;
822:
823: if ( !p )
824: *pr = 0;
825: else {
826: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
827: i = ITOS(m->c); STOQ(i,q); t = (P)q;
828: for ( i = 0, d = m->dl, tvl = dvl;
829: i < n; tvl = NEXT(tvl), i++ ) {
830: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
831: mulp(vl,t,u,&w); t = w;
832: }
833: addp(vl,s,t,&u); s = u;
834: }
835: *pr = s;
836: }
837: }
838:
839: void _dp_red_mod(p1,p2,mod,rp)
840: DP p1,p2;
841: int mod;
842: DP *rp;
843: {
844: int i,n;
845: DL d1,d2,d;
846: MP m;
847: DP t,s;
848: int c,c1;
849:
850: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
851: NEWDL(d,n); d->td = d1->td - d2->td;
852: for ( i = 0; i < n; i++ )
853: d->d[i] = d1->d[i]-d2->d[i];
854: c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
855: NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
856: MKDP(n,m,s); s->sugar = d->td;
1.2 noro 857: _mulmd(CO,mod,s,p2,&t); _addmd(CO,mod,p1,t,rp);
1.1 noro 858: }
859:
860: void _dp_mod(p,mod,subst,rp)
861: DP p;
862: int mod;
863: NODE subst;
864: DP *rp;
865: {
866: MP m,mr,mr0;
867: P t,s;
868: NODE tn;
869:
870: if ( !p )
871: *rp = 0;
872: else {
873: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
874: for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) {
875: substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t);
876: s = t;
877: }
878: ptomp(mod,s,&t);
879: if ( t ) {
880: NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl;
881: }
882: }
883: if ( mr0 ) {
884: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
885: } else
886: *rp = 0;
887: }
888: }
889:
1.2 noro 890: void _dp_monic(p,mod,rp)
891: DP p;
892: int mod;
893: DP *rp;
894: {
895: MP m,mr,mr0;
896: int c,c1;
897: NODE tn;
898:
899: if ( !p )
900: *rp = 0;
901: else {
902: c = invm(ITOS(BDY(p)->c),mod);
903: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
904: c1 = dmar(ITOS(m->c),c,0,mod);
905: NEXTMP(mr0,mr); mr->c = STOI(c1); mr->dl = m->dl;
906: }
907: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
908: }
909: }
910:
1.1 noro 911: void _dp_sp_mod(p1,p2,mod,rp)
912: DP p1,p2;
913: int mod;
914: DP *rp;
915: {
916: int i,n,td;
917: int *w;
918: DL d1,d2,d;
919: MP m;
920: DP t,s,u;
921:
922: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
923: w = (int *)ALLOCA(n*sizeof(int));
924: for ( i = 0, td = 0; i < n; i++ ) {
925: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
926: }
927: NEWDL(d,n); d->td = td - d1->td;
928: for ( i = 0; i < n; i++ )
929: d->d[i] = w[i] - d1->d[i];
930: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
1.2 noro 931: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p1,&t);
1.1 noro 932: NEWDL(d,n); d->td = td - d2->td;
933: for ( i = 0; i < n; i++ )
934: d->d[i] = w[i] - d2->d[i];
935: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
1.2 noro 936: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p2,&u);
1.1 noro 937: _addmd(CO,mod,t,u,rp);
938: }
939:
940: void _dp_sp_component_mod(p1,p2,mod,f1,f2)
941: DP p1,p2;
942: int mod;
943: DP *f1,*f2;
944: {
945: int i,n,td;
946: int *w;
947: DL d1,d2,d;
948: MP m;
949: DP t,s,u;
950:
951: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
952: w = (int *)ALLOCA(n*sizeof(int));
953: for ( i = 0, td = 0; i < n; i++ ) {
954: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
955: }
956: NEWDL(d,n); d->td = td - d1->td;
957: for ( i = 0; i < n; i++ )
958: d->d[i] = w[i] - d1->d[i];
959: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
1.2 noro 960: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p1,f1);
1.1 noro 961: NEWDL(d,n); d->td = td - d2->td;
962: for ( i = 0; i < n; i++ )
963: d->d[i] = w[i] - d2->d[i];
964: NEWMP(m); m->dl = d; m->c = BDY(p1)->c; NEXT(m) = 0;
1.2 noro 965: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p2,f2);
1.1 noro 966: }
967:
968: void _printdp(d)
969: DP d;
970: {
971: int n,i;
972: MP m;
973: DL dl;
974:
975: if ( !d ) {
976: printf("0"); return;
977: }
978: for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) {
979: printf("%d*<<",ITOS(m->c));
980: for ( i = 0, dl = m->dl; i < n-1; i++ )
981: printf("%d,",dl->d[i]);
982: printf("%d",dl->d[i]);
983: printf(">>");
984: if ( NEXT(m) )
985: printf("+");
986: }
987: }