![[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.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.1.1.1 1999/12/03 07:39:08 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;
! 369: int i,j,a,b,k,l,n2,s,min,h;
! 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 */
! 383: NEWDL(d,n); d->d[n-1] = d0->d[n-1]+d1->d[n-1]; d->td = d->d[n-1];
! 384: NEWMP(mr); mr->c = (P)ONEM; mr->dl = d;
! 385: MKDP(n,mr,r); r->sugar = d->d[n-1];
! 386:
! 387: for ( i = 0; i < n2; i++ ) {
! 388: a = d0->d[i]; b = d1->d[n2+i];
! 389: k = d0->d[n2+i]; l = d1->d[i];
! 390: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 391: h = a+k+l+b;
! 392: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 393: min = MIN(k,l);
! 394:
! 395: if ( min+1 > tablen ) {
! 396: if ( tab ) GC_free(tab);
! 397: tab = (int *)MALLOC((min+1)*sizeof(int));
! 398: tablen = min+1;
! 399: }
! 400: mkwcm(k,l,mod,tab);
! 401: for ( mr0 = 0, j = 0; j <= min; j++ ) {
! 402: NEXTMP(mr0,mr);
! 403: NEWDL(d,n);
! 404: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 405: d->td = h;
! 406: d->d[n-1] = h-(d->d[i]+d->d[n2+i]);
! 407: STOMQ(tab[j],mq); mr->c = (P)mq;
! 408: mr->dl = d;
! 409: }
! 410: bzero(tab,(min+1)*sizeof(int));
! 411: if ( mr0 )
! 412: NEXT(mr) = 0;
! 413: MKDP(n,mr0,t);
! 414: if ( t )
! 415: t->sugar = h;
! 416: comm_mulmd(vl,mod,r,t,&t1); r = t1;
! 417: }
! 418: } else
! 419: for ( i = 0, r = (DP)ONEM; i < n2; i++ ) {
! 420: a = d0->d[i]; b = d1->d[n2+i];
! 421: k = d0->d[n2+i]; l = d1->d[i];
! 422: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 423: min = MIN(k,l);
! 424:
! 425: if ( min+1 > tablen ) {
! 426: if ( tab ) GC_free(tab);
! 427: tab = (int *)MALLOC((min+1)*sizeof(int));
! 428: tablen = min+1;
! 429: }
! 430: mkwcm(k,l,mod,tab);
! 431: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
! 432: NEXTMP(mr0,mr);
! 433: NEWDL(d,n);
! 434: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 435: d->td = d->d[i]+d->d[n2+i]; /* XXX */
! 436: s = MAX(s,d->td); /* XXX */
! 437: STOMQ(tab[j],mq); mr->c = (P)mq;
! 438: mr->dl = d;
! 439: }
! 440: bzero(tab,(min+1)*sizeof(int));
! 441: if ( mr0 )
! 442: NEXT(mr) = 0;
! 443: MKDP(n,mr0,t);
! 444: if ( t )
! 445: t->sugar = s;
! 446: comm_mulmd(vl,mod,r,t,&t1); r = t1;
! 447: }
! 448: mulmdc(vl,mod,r,c,pr);
! 449: }
! 450: }
! 451:
1.1 noro 452: void mulmdc(vl,mod,p,c,pr)
453: VL vl;
454: int mod;
455: DP p;
456: P c;
457: DP *pr;
458: {
459: MP m,mr,mr0;
460: int t;
461: MQ q;
462:
463: if ( !p || !c )
464: *pr = 0;
465: else {
466: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
467: NEXTMP(mr0,mr);
468: if ( NUM(C(m)) && NUM(c) ) {
469: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
470: STOMQ(t,q); C(mr) = (P)q;
471: } else
472: mulmp(vl,mod,C(m),c,&C(mr));
473: mr->dl = m->dl;
474: }
475: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
476: if ( *pr )
477: (*pr)->sugar = p->sugar;
478: }
479: }
480:
481: void divsmdc(vl,mod,p,c,pr)
482: VL vl;
483: int mod;
484: DP p;
485: P c;
486: DP *pr;
487: {
488: MP m,mr,mr0;
489:
490: if ( !c )
491: error("disvsdc : division by 0");
492: else if ( !p )
493: *pr = 0;
494: else {
495: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
496: NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl;
497: }
498: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
499: if ( *pr )
500: (*pr)->sugar = p->sugar;
501: }
502: }
503:
504: #define MKDPM(n,m,d) (NEWDP(d),(d)->nv=(n),BDY(d)=(m))
505:
506: void _mdtop(vl,mod,dvl,p,pr)
507: VL vl,dvl;
508: int mod;
509: DP p;
510: P *pr;
511: {
512: int n,i;
513: DL d;
514: MP m;
515: P r,s,t,u,w;
516: Q q;
517: MQ tq;
518: VL tvl;
519:
520: if ( !p )
521: *pr = 0;
522: else {
523: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
524: STOMQ(ITOS(C(m)),tq); t = (P)tq;
525: for ( i = 0, d = m->dl, tvl = dvl;
526: i < n; tvl = NEXT(tvl), i++ ) {
527: MKV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
528: mulmp(vl,mod,t,u,&w); t = w;
529: }
530: addmp(vl,mod,s,t,&u); s = u;
531: }
532: mptop(s,pr);
533: }
534: }
535:
536: void _addmd(vl,mod,p1,p2,pr)
537: VL vl;
538: int mod;
539: DP p1,p2,*pr;
540: {
541: int n;
542: MP m1,m2,mr,mr0;
543: int t;
544:
545: if ( !p1 )
546: *pr = p2;
547: else if ( !p2 )
548: *pr = p1;
549: else {
550: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
551: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
552: case 0:
553: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
554: if ( t < 0 )
555: t += mod;
556: if ( t ) {
557: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = STOI(t);
558: }
559: m1 = NEXT(m1); m2 = NEXT(m2); break;
560: case 1:
561: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
562: m1 = NEXT(m1); break;
563: case -1:
564: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
565: m2 = NEXT(m2); break;
566: }
567: if ( !mr0 )
568: if ( m1 )
569: mr0 = m1;
570: else if ( m2 )
571: mr0 = m2;
572: else {
573: *pr = 0;
574: return;
575: }
576: else if ( m1 )
577: NEXT(mr) = m1;
578: else if ( m2 )
579: NEXT(mr) = m2;
580: else
581: NEXT(mr) = 0;
582: MKDPM(NV(p1),mr0,*pr);
583: if ( *pr )
584: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
585: }
586: }
587:
588: void _submd(vl,mod,p1,p2,pr)
589: VL vl;
590: int mod;
591: DP p1,p2,*pr;
592: {
593: DP t;
594:
595: if ( !p2 )
596: *pr = p1;
597: else {
598: _chsgnmd(mod,p2,&t); _addmd(vl,mod,p1,t,pr);
599: }
600: }
601:
602: void _chsgnmd(mod,p,pr)
603: int mod;
604: DP p,*pr;
605: {
606: MP m,mr,mr0;
607:
608: if ( !p )
609: *pr = 0;
610: else {
611: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
612: NEXTMP(mr0,mr); C(mr) = STOI(mod - ITOS(C(m))); mr->dl = m->dl;
613: }
614: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
615: if ( *pr )
616: (*pr)->sugar = p->sugar;
617: }
618: }
619:
620: void _mulmd(vl,mod,p1,p2,pr)
621: VL vl;
622: int mod;
623: DP p1,p2,*pr;
624: {
1.2 ! noro 625: if ( !do_weyl )
! 626: _comm_mulmd(vl,mod,p1,p2,pr);
! 627: else
! 628: _weyl_mulmd(vl,mod,p1,p2,pr);
! 629: }
! 630:
! 631: void _comm_mulmd(vl,mod,p1,p2,pr)
! 632: VL vl;
! 633: int mod;
! 634: DP p1,p2,*pr;
! 635: {
1.1 noro 636: MP m;
637: DP s,t,u;
1.2 ! noro 638: int i,l,l1;
! 639: static MP *w;
! 640: static int wlen;
1.1 noro 641:
642: if ( !p1 || !p2 )
643: *pr = 0;
644: else {
1.2 ! noro 645: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
! 646: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 647: if ( l1 < l ) {
! 648: t = p1; p1 = p2; p2 = t;
! 649: l = l1;
! 650: }
! 651: if ( l > wlen ) {
! 652: if ( w ) GC_free(w);
! 653: w = (MP *)MALLOC(l*sizeof(MP));
! 654: wlen = l;
! 655: }
! 656: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 657: w[i] = m;
! 658: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 659: _mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
1.1 noro 660: }
1.2 ! noro 661: bzero(w,l*sizeof(MP));
! 662: *pr = s;
! 663: }
! 664: }
! 665:
! 666: void _weyl_mulmd(vl,mod,p1,p2,pr)
! 667: VL vl;
! 668: int mod;
! 669: DP p1,p2,*pr;
! 670: {
! 671: MP m;
! 672: DP s,t,u;
! 673: int i,l,l1;
! 674: static MP *w;
! 675: static int wlen;
! 676:
! 677: if ( !p1 || !p2 )
! 678: *pr = 0;
! 679: else {
! 680: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 681: if ( l > wlen ) {
! 682: if ( w ) GC_free(w);
! 683: w = (MP *)MALLOC(l*sizeof(MP));
! 684: wlen = l;
! 685: }
! 686: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 687: w[i] = m;
! 688: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 689: _weyl_mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
! 690: }
! 691: bzero(w,l*sizeof(MP));
1.1 noro 692: *pr = s;
693: }
694: }
695:
696: void _mulmdm(vl,mod,p,m0,pr)
697: VL vl;
698: int mod;
699: DP p;
700: MP m0;
701: DP *pr;
702: {
703: MP m,mr,mr0;
704: DL d;
705: int c,n,r;
706:
707: if ( !p )
708: *pr = 0;
709: else {
710: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
711: m; m = NEXT(m) ) {
712: NEXTMP(mr0,mr);
713: C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
714: adddl(n,m->dl,d,&mr->dl);
715: }
716: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
717: if ( *pr )
718: (*pr)->sugar = p->sugar + m0->dl->td;
719: }
720: }
721:
1.2 ! noro 722: void _weyl_mulmdm(vl,mod,p,m0,pr)
! 723: VL vl;
! 724: int mod;
! 725: DP p;
! 726: MP m0;
! 727: DP *pr;
! 728: {
! 729: DP r,t,t1;
! 730: MP m;
! 731: int n,l,i;
! 732: static MP *w;
! 733: static int wlen;
! 734:
! 735: if ( !p )
! 736: *pr = 0;
! 737: else {
! 738: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 739: if ( l > wlen ) {
! 740: if ( w ) GC_free(w);
! 741: w = (MP *)MALLOC(l*sizeof(MP));
! 742: wlen = l;
! 743: }
! 744: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 745: w[i] = m;
! 746: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
! 747: _weyl_mulmmm(vl,mod,w[i],m0,n,&t);
! 748: _addmd(vl,mod,r,t,&t1); r = t1;
! 749: }
! 750: bzero(w,l*sizeof(MP));
! 751: if ( r )
! 752: r->sugar = p->sugar + m0->dl->td;
! 753: *pr = r;
! 754: }
! 755: }
! 756:
! 757: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
! 758:
! 759: void _weyl_mulmmm(vl,mod,m0,m1,n,pr)
! 760: VL vl;
! 761: int mod;
! 762: MP m0,m1;
! 763: int n;
! 764: DP *pr;
! 765: {
! 766: MP m,mr,mr0;
! 767: DP r,t,t1;
! 768: int c,c0,c1,cc;
! 769: DL d,d0,d1;
! 770: int i,j,a,b,k,l,n2,s,min,h;
! 771: static int *tab;
! 772: static int tablen;
! 773:
! 774: if ( !m0 || !m1 )
! 775: *pr = 0;
! 776: else {
! 777: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
! 778: c = dmar(c0,c1,0,mod);
! 779: d0 = m0->dl; d1 = m1->dl;
! 780: n2 = n>>1;
! 781:
! 782: if ( n & 1 ) {
! 783: /* offset of h-degree */
! 784: NEWDL(d,n); d->d[n-1] = d0->d[n-1]+d1->d[n-1]; d->td = d->d[n-1];
! 785: NEWMP(mr); mr->c = STOI(c); mr->dl = d;
! 786: MKDPM(n,mr,r); r->sugar = d->d[n-1];
! 787:
! 788: /* homogenized computation; dx-xd=h^2 */
! 789: for ( i = 0; i < n2; i++ ) {
! 790: a = d0->d[i]; b = d1->d[n2+i];
! 791: k = d0->d[n2+i]; l = d1->d[i];
! 792: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 793: h = a+k+l+b;
! 794: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 795: min = MIN(k,l);
! 796:
! 797: if ( min+1 > tablen ) {
! 798: if ( tab ) GC_free(tab);
! 799: tab = (int *)MALLOC((min+1)*sizeof(int));
! 800: tablen = min+1;
! 801: }
! 802: mkwcm(k,l,mod,tab);
! 803: for ( mr0 = 0, j = 0; j <= min; j++ ) {
! 804: NEXTMP(mr0,mr);
! 805: NEWDL(d,n);
! 806: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 807: d->td = h;
! 808: d->d[n-1] = h-(d->d[i]+d->d[n2+i]);
! 809: mr->c = STOI(tab[j]);
! 810: mr->dl = d;
! 811: }
! 812: bzero(tab,(min+1)*sizeof(int));
! 813: if ( mr0 )
! 814: NEXT(mr) = 0;
! 815: MKDP(n,mr0,t);
! 816: if ( t )
! 817: t->sugar = h;
! 818: _comm_mulmd(vl,mod,r,t,&t1); r = t1;
! 819: }
! 820: } else {
! 821: NEWDL(d,n); d->td = 0;
! 822: NEWMP(mr); mr->c = STOI(c); mr->dl = d;
! 823: MKDPM(n,mr,r); r->sugar = 0;
! 824:
! 825: for ( i = 0; i < n2; i++ ) {
! 826: a = d0->d[i]; b = d1->d[n2+i];
! 827: k = d0->d[n2+i]; l = d1->d[i];
! 828: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 829: min = MIN(k,l);
! 830:
! 831: if ( min+1 > tablen ) {
! 832: if ( tab ) GC_free(tab);
! 833: tab = (int *)MALLOC((min+1)*sizeof(int));
! 834: tablen = min+1;
! 835: }
! 836: mkwcm(k,l,mod,tab);
! 837: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
! 838: NEXTMP(mr0,mr);
! 839: NEWDL(d,n);
! 840: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 841: d->td = d->d[i]+d->d[n2+i]; /* XXX */
! 842: s = MAX(s,d->td); /* XXX */
! 843: mr->c = STOI(tab[j]);
! 844: mr->dl = d;
! 845: }
! 846: bzero(tab,(min+1)*sizeof(int));
! 847: if ( mr0 )
! 848: NEXT(mr) = 0;
! 849: MKDP(n,mr0,t);
! 850: if ( t )
! 851: t->sugar = s;
! 852: _comm_mulmd(vl,mod,r,t,&t1); r = t1;
! 853: }
! 854: }
! 855: *pr = r;
! 856: }
! 857: }
! 858:
1.1 noro 859: void _dtop_mod(vl,dvl,p,pr)
860: VL vl,dvl;
861: DP p;
862: P *pr;
863: {
864: int n,i;
865: DL d;
866: MP m;
867: P r,s,t,u,w;
868: Q q;
869: VL tvl;
870:
871: if ( !p )
872: *pr = 0;
873: else {
874: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
875: i = ITOS(m->c); STOQ(i,q); t = (P)q;
876: for ( i = 0, d = m->dl, tvl = dvl;
877: i < n; tvl = NEXT(tvl), i++ ) {
878: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
879: mulp(vl,t,u,&w); t = w;
880: }
881: addp(vl,s,t,&u); s = u;
882: }
883: *pr = s;
884: }
885: }
886:
887: void _dp_red_mod(p1,p2,mod,rp)
888: DP p1,p2;
889: int mod;
890: DP *rp;
891: {
892: int i,n;
893: DL d1,d2,d;
894: MP m;
895: DP t,s;
896: int c,c1;
897:
898: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
899: NEWDL(d,n); d->td = d1->td - d2->td;
900: for ( i = 0; i < n; i++ )
901: d->d[i] = d1->d[i]-d2->d[i];
902: c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
903: NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
904: MKDP(n,m,s); s->sugar = d->td;
1.2 ! noro 905: _mulmd(CO,mod,s,p2,&t); _addmd(CO,mod,p1,t,rp);
1.1 noro 906: }
907:
908: void _dp_mod(p,mod,subst,rp)
909: DP p;
910: int mod;
911: NODE subst;
912: DP *rp;
913: {
914: MP m,mr,mr0;
915: P t,s;
916: NODE tn;
917:
918: if ( !p )
919: *rp = 0;
920: else {
921: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
922: for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) {
923: substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t);
924: s = t;
925: }
926: ptomp(mod,s,&t);
927: if ( t ) {
928: NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl;
929: }
930: }
931: if ( mr0 ) {
932: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
933: } else
934: *rp = 0;
935: }
936: }
937:
1.2 ! noro 938: void _dp_monic(p,mod,rp)
! 939: DP p;
! 940: int mod;
! 941: DP *rp;
! 942: {
! 943: MP m,mr,mr0;
! 944: int c,c1;
! 945: NODE tn;
! 946:
! 947: if ( !p )
! 948: *rp = 0;
! 949: else {
! 950: c = invm(ITOS(BDY(p)->c),mod);
! 951: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 952: c1 = dmar(ITOS(m->c),c,0,mod);
! 953: NEXTMP(mr0,mr); mr->c = STOI(c1); mr->dl = m->dl;
! 954: }
! 955: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 956: }
! 957: }
! 958:
1.1 noro 959: void _dp_sp_mod(p1,p2,mod,rp)
960: DP p1,p2;
961: int mod;
962: DP *rp;
963: {
964: int i,n,td;
965: int *w;
966: DL d1,d2,d;
967: MP m;
968: DP t,s,u;
969:
970: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
971: w = (int *)ALLOCA(n*sizeof(int));
972: for ( i = 0, td = 0; i < n; i++ ) {
973: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
974: }
975: NEWDL(d,n); d->td = td - d1->td;
976: for ( i = 0; i < n; i++ )
977: d->d[i] = w[i] - d1->d[i];
978: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
1.2 ! noro 979: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p1,&t);
1.1 noro 980: NEWDL(d,n); d->td = td - d2->td;
981: for ( i = 0; i < n; i++ )
982: d->d[i] = w[i] - d2->d[i];
983: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
1.2 ! noro 984: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p2,&u);
1.1 noro 985: _addmd(CO,mod,t,u,rp);
986: }
987:
988: void _dp_sp_component_mod(p1,p2,mod,f1,f2)
989: DP p1,p2;
990: int mod;
991: DP *f1,*f2;
992: {
993: int i,n,td;
994: int *w;
995: DL d1,d2,d;
996: MP m;
997: DP t,s,u;
998:
999: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1000: w = (int *)ALLOCA(n*sizeof(int));
1001: for ( i = 0, td = 0; i < n; i++ ) {
1002: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
1003: }
1004: NEWDL(d,n); d->td = td - d1->td;
1005: for ( i = 0; i < n; i++ )
1006: d->d[i] = w[i] - d1->d[i];
1007: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
1.2 ! noro 1008: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p1,f1);
1.1 noro 1009: NEWDL(d,n); d->td = td - d2->td;
1010: for ( i = 0; i < n; i++ )
1011: d->d[i] = w[i] - d2->d[i];
1012: NEWMP(m); m->dl = d; m->c = BDY(p1)->c; NEXT(m) = 0;
1.2 ! noro 1013: MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p2,f2);
1.1 noro 1014: }
1015:
1016: void _printdp(d)
1017: DP d;
1018: {
1019: int n,i;
1020: MP m;
1021: DL dl;
1022:
1023: if ( !d ) {
1024: printf("0"); return;
1025: }
1026: for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) {
1027: printf("%d*<<",ITOS(m->c));
1028: for ( i = 0, dl = m->dl; i < n-1; i++ )
1029: printf("%d,",dl->d[i]);
1030: printf("%d",dl->d[i]);
1031: printf(">>");
1032: if ( NEXT(m) )
1033: printf("+");
1034: }
1035: }