![[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.5 ! noro 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.4 2000/07/13 05:09:01 noro Exp $
! 49: */
1.1 noro 50: #include "ca.h"
51: #include "inline.h"
52:
53: #define NV(p) ((p)->nv)
54: #define C(p) ((p)->c)
55: #if 0
56: #define ITOS(p) (((unsigned int)(p))>>1)
57: #define STOI(i) ((P)((((unsigned int)(i))<<1)|1))
58: #else
59: #define ITOS(p) (((unsigned int)(p))&0x7fffffff)
60: #define STOI(i) ((P)((unsigned int)(i)|0x80000000))
61: #endif
62:
63: extern int (*cmpdl)();
1.2 noro 64: extern int do_weyl;
65:
66: void comm_mulmd();
67: void weyl_mulmd();
68: void weyl_mulmdm();
69: void weyl_mulmmm();
70: void _comm_mulmd();
71: void _weyl_mulmd();
72: void _weyl_mulmmm();
73: void _weyl_mulmdm();
1.1 noro 74:
75: void ptomd(vl,mod,dvl,p,pr)
76: VL vl,dvl;
77: int mod;
78: P p;
79: DP *pr;
80: {
81: P t;
82:
83: ptomp(mod,p,&t);
84: mptomd(vl,mod,dvl,t,pr);
85: }
86:
87: void mptomd(vl,mod,dvl,p,pr)
88: VL vl,dvl;
89: int mod;
90: P p;
91: DP *pr;
92: {
93: int n,i;
94: VL tvl;
95: V v;
96: DL d;
97: MP m;
98: DCP dc;
99: DP r,s,t,u;
100: P x,c;
101:
102: if ( !p )
103: *pr = 0;
104: else {
105: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
106: if ( NUM(p) ) {
107: NEWDL(d,n);
108: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr);
109: } else {
110: for ( i = 0, tvl = dvl, v = VR(p);
111: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
112: if ( !tvl ) {
113: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
114: mptomd(vl,mod,dvl,COEF(dc),&t); pwrmp(vl,mod,x,DEG(dc),&c);
115: mulmdc(vl,mod,t,c,&r); addmd(vl,mod,r,s,&t); s = t;
116: }
117: *pr = s;
118: } else {
119: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
120: mptomd(vl,mod,dvl,COEF(dc),&t);
121: NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
122: NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u);
1.2 noro 123: comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t;
1.1 noro 124: }
125: *pr = s;
126: }
127: }
128: }
129: }
130:
131: void mdtop(vl,mod,dvl,p,pr)
132: VL vl,dvl;
133: int mod;
134: DP p;
135: P *pr;
136: {
137: int n,i;
138: DL d;
139: MP m;
140: P r,s,t,u,w;
141: Q q;
142: VL tvl;
143:
144: if ( !p )
145: *pr = 0;
146: else {
147: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
148: for ( i = 0, t = C(m), d = m->dl, tvl = dvl;
149: i < n; tvl = NEXT(tvl), i++ ) {
150: MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
151: mulmp(vl,mod,t,u,&w); t = w;
152: }
153: addmp(vl,mod,s,t,&u); s = u;
154: }
155: mptop(s,pr);
156: }
157: }
158:
159: void addmd(vl,mod,p1,p2,pr)
160: VL vl;
161: int mod;
162: DP p1,p2,*pr;
163: {
164: int n;
165: MP m1,m2,mr,mr0;
166: P t;
167: int tmp;
168: MQ q;
169:
170: if ( !p1 )
171: *pr = p2;
172: else if ( !p2 )
173: *pr = p1;
174: else {
175: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
176: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
177: case 0:
178: if ( NUM(C(m1)) && NUM(C(m2)) ) {
179: tmp = (CONT((MQ)C(m1))+CONT((MQ)C(m2))) - mod;
180: if ( tmp < 0 )
181: tmp += mod;
182: if ( tmp ) {
183: STOMQ(tmp,q); t = (P)q;
184: } else
185: t = 0;
186: } else
187: addmp(vl,mod,C(m1),C(m2),&t);
188: if ( t ) {
189: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
190: }
191: m1 = NEXT(m1); m2 = NEXT(m2); break;
192: case 1:
193: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
194: m1 = NEXT(m1); break;
195: case -1:
196: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
197: m2 = NEXT(m2); break;
198: }
199: if ( !mr0 )
200: if ( m1 )
201: mr0 = m1;
202: else if ( m2 )
203: mr0 = m2;
204: else {
205: *pr = 0;
206: return;
207: }
208: else if ( m1 )
209: NEXT(mr) = m1;
210: else if ( m2 )
211: NEXT(mr) = m2;
212: else
213: NEXT(mr) = 0;
214: MKDP(NV(p1),mr0,*pr);
215: if ( *pr )
216: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
217: }
218: }
219:
220: void submd(vl,mod,p1,p2,pr)
221: VL vl;
222: int mod;
223: DP p1,p2,*pr;
224: {
225: DP t;
226:
227: if ( !p2 )
228: *pr = p1;
229: else {
230: chsgnmd(mod,p2,&t); addmd(vl,mod,p1,t,pr);
231: }
232: }
233:
234: void chsgnmd(mod,p,pr)
235: int mod;
236: DP p,*pr;
237: {
238: MP m,mr,mr0;
239:
240: if ( !p )
241: *pr = 0;
242: else {
243: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
244: NEXTMP(mr0,mr); chsgnmp(mod,C(m),&C(mr)); mr->dl = m->dl;
245: }
246: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
247: if ( *pr )
248: (*pr)->sugar = p->sugar;
249: }
250: }
251:
252: void mulmd(vl,mod,p1,p2,pr)
253: VL vl;
254: int mod;
255: DP p1,p2,*pr;
256: {
1.2 noro 257: if ( !do_weyl )
258: comm_mulmd(vl,mod,p1,p2,pr);
259: else
260: weyl_mulmd(vl,mod,p1,p2,pr);
261: }
262:
263: void comm_mulmd(vl,mod,p1,p2,pr)
264: VL vl;
265: int mod;
266: DP p1,p2,*pr;
267: {
268: MP m;
269: DP s,t,u;
270: int i,l,l1;
271: static MP *w;
272: static int wlen;
273:
274: if ( !p1 || !p2 )
275: *pr = 0;
276: else if ( OID(p1) <= O_P )
277: mulmdc(vl,mod,p2,(P)p1,pr);
278: else if ( OID(p2) <= O_P )
279: mulmdc(vl,mod,p1,(P)p2,pr);
280: else {
281: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
282: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
283: if ( l1 < l ) {
284: t = p1; p1 = p2; p2 = t;
285: l = l1;
286: }
287: if ( l > wlen ) {
288: if ( w ) GC_free(w);
289: w = (MP *)MALLOC(l*sizeof(MP));
290: wlen = l;
291: }
292: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
293: w[i] = m;
294: for ( s = 0, i = l-1; i >= 0; i-- ) {
295: mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
296: }
297: bzero(w,l*sizeof(MP));
298: *pr = s;
299: }
300: }
301:
302: void weyl_mulmd(vl,mod,p1,p2,pr)
303: VL vl;
304: int mod;
305: DP p1,p2,*pr;
306: {
1.1 noro 307: MP m;
308: DP s,t,u;
1.2 noro 309: int i,l,l1;
310: static MP *w;
311: static int wlen;
1.1 noro 312:
313: if ( !p1 || !p2 )
314: *pr = 0;
315: else if ( OID(p1) <= O_P )
316: mulmdc(vl,mod,p2,(P)p1,pr);
317: else if ( OID(p2) <= O_P )
318: mulmdc(vl,mod,p1,(P)p2,pr);
319: else {
1.2 noro 320: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
321: if ( l > wlen ) {
322: if ( w ) GC_free(w);
323: w = (MP *)MALLOC(l*sizeof(MP));
324: wlen = l;
325: }
326: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
327: w[i] = m;
328: for ( s = 0, i = l-1; i >= 0; i-- ) {
329: weyl_mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
1.1 noro 330: }
1.2 noro 331: bzero(w,l*sizeof(MP));
1.1 noro 332: *pr = s;
333: }
334: }
335:
336: void mulmdm(vl,mod,p,m0,pr)
337: VL vl;
338: int mod;
339: DP p;
340: MP m0;
341: DP *pr;
342: {
343: MP m,mr,mr0;
344: P c;
345: DL d;
346: int n,t;
347: MQ q;
348:
349: if ( !p )
350: *pr = 0;
351: else {
352: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
353: m; m = NEXT(m) ) {
354: NEXTMP(mr0,mr);
355: if ( NUM(C(m)) && NUM(c) ) {
356: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
357: STOMQ(t,q); C(mr) = (P)q;
358: } else
359: mulmp(vl,mod,C(m),c,&C(mr));
360: adddl(n,m->dl,d,&mr->dl);
361: }
362: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
363: if ( *pr )
364: (*pr)->sugar = p->sugar + m0->dl->td;
365: }
366: }
367:
1.2 noro 368: void weyl_mulmdm(vl,mod,p,m0,pr)
369: VL vl;
370: int mod;
371: DP p;
372: MP m0;
373: DP *pr;
374: {
375: DP r,t,t1;
376: MP m;
377: int n,l,i;
378: static MP *w;
379: static int wlen;
380:
381: if ( !p )
382: *pr = 0;
383: else {
384: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
385: if ( l > wlen ) {
386: if ( w ) GC_free(w);
387: w = (MP *)MALLOC(l*sizeof(MP));
388: wlen = l;
389: }
390: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
391: w[i] = m;
392: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
393: weyl_mulmmm(vl,mod,w[i],m0,n,&t);
394: addmd(vl,mod,r,t,&t1); r = t1;
395: }
396: bzero(w,l*sizeof(MP));
397: if ( r )
398: r->sugar = p->sugar + m0->dl->td;
399: *pr = r;
400: }
401: }
402:
403: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
404:
405: void weyl_mulmmm(vl,mod,m0,m1,n,pr)
406: VL vl;
407: int mod;
408: MP m0,m1;
409: int n;
410: DP *pr;
411: {
412: MP m,mr,mr0;
413: MQ mq;
414: DP r,t,t1;
415: P c,c0,c1,cc;
416: DL d,d0,d1;
1.3 noro 417: int i,j,a,b,k,l,n2,s,min;
1.2 noro 418: static int *tab;
419: static int tablen;
420:
421: if ( !m0 || !m1 )
422: *pr = 0;
423: else {
424: c0 = C(m0); c1 = C(m1);
425: mulmp(vl,mod,c0,c1,&c);
426: d0 = m0->dl; d1 = m1->dl;
427: n2 = n>>1;
428: if ( n & 1 ) {
429: /* homogenized computation; dx-xd=h^2 */
430: /* offset of h-degree */
1.3 noro 431: NEWDL(d,n);
432: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1.4 noro 433: NEWMP(mr); mr->c = (P)ONEM; mr->dl = d; NEXT(mr) = 0;
1.3 noro 434: MKDP(n,mr,r); r->sugar = d->td;
435: } else
436: r = (DP)ONEM;
437: for ( i = 0; i < n2; i++ ) {
438: a = d0->d[i]; b = d1->d[n2+i];
439: k = d0->d[n2+i]; l = d1->d[i];
440: /* degree of xi^a*(Di^k*xi^l)*Di^b */
441: s = a+k+l+b;
442: /* compute xi^a*(Di^k*xi^l)*Di^b */
443: min = MIN(k,l);
444:
445: if ( min+1 > tablen ) {
446: if ( tab ) GC_free(tab);
447: tab = (int *)MALLOC((min+1)*sizeof(int));
448: tablen = min+1;
449: }
450: mkwcm(k,l,mod,tab);
451: if ( n & 1 )
1.2 noro 452: for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3 noro 453: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 454: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.3 noro 455: d->td = s;
456: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
457: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 458: }
1.3 noro 459: else
1.2 noro 460: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3 noro 461: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 462: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
463: d->td = d->d[i]+d->d[n2+i]; /* XXX */
464: s = MAX(s,d->td); /* XXX */
1.3 noro 465: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 466: }
1.3 noro 467: bzero(tab,(min+1)*sizeof(int));
468: if ( mr0 )
469: NEXT(mr) = 0;
470: MKDP(n,mr0,t);
471: if ( t )
472: t->sugar = s;
473: comm_mulmd(vl,mod,r,t,&t1); r = t1;
474: }
1.2 noro 475: mulmdc(vl,mod,r,c,pr);
476: }
477: }
478:
1.1 noro 479: void mulmdc(vl,mod,p,c,pr)
480: VL vl;
481: int mod;
482: DP p;
483: P c;
484: DP *pr;
485: {
486: MP m,mr,mr0;
487: int t;
488: MQ q;
489:
490: if ( !p || !c )
491: *pr = 0;
492: else {
493: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
494: NEXTMP(mr0,mr);
495: if ( NUM(C(m)) && NUM(c) ) {
496: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
497: STOMQ(t,q); C(mr) = (P)q;
498: } else
499: mulmp(vl,mod,C(m),c,&C(mr));
500: mr->dl = m->dl;
501: }
502: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
503: if ( *pr )
504: (*pr)->sugar = p->sugar;
505: }
506: }
507:
508: void divsmdc(vl,mod,p,c,pr)
509: VL vl;
510: int mod;
511: DP p;
512: P c;
513: DP *pr;
514: {
515: MP m,mr,mr0;
516:
517: if ( !c )
518: error("disvsdc : division by 0");
519: else if ( !p )
520: *pr = 0;
521: else {
522: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
523: NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl;
524: }
525: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
526: if ( *pr )
527: (*pr)->sugar = p->sugar;
528: }
529: }
530:
531: void _mdtop(vl,mod,dvl,p,pr)
532: VL vl,dvl;
533: int mod;
534: DP p;
535: P *pr;
536: {
537: int n,i;
538: DL d;
539: MP m;
540: P r,s,t,u,w;
541: Q q;
542: MQ tq;
543: VL tvl;
544:
545: if ( !p )
546: *pr = 0;
547: else {
548: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
549: STOMQ(ITOS(C(m)),tq); t = (P)tq;
550: for ( i = 0, d = m->dl, tvl = dvl;
551: i < n; tvl = NEXT(tvl), i++ ) {
552: MKV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
553: mulmp(vl,mod,t,u,&w); t = w;
554: }
555: addmp(vl,mod,s,t,&u); s = u;
556: }
557: mptop(s,pr);
558: }
559: }
560:
561: void _addmd(vl,mod,p1,p2,pr)
562: VL vl;
563: int mod;
564: DP p1,p2,*pr;
565: {
566: int n;
567: MP m1,m2,mr,mr0;
568: int t;
569:
570: if ( !p1 )
571: *pr = p2;
572: else if ( !p2 )
573: *pr = p1;
574: else {
575: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
576: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
577: case 0:
578: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
579: if ( t < 0 )
580: t += mod;
581: if ( t ) {
582: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = STOI(t);
583: }
584: m1 = NEXT(m1); m2 = NEXT(m2); break;
585: case 1:
586: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
587: m1 = NEXT(m1); break;
588: case -1:
589: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
590: m2 = NEXT(m2); break;
591: }
592: if ( !mr0 )
593: if ( m1 )
594: mr0 = m1;
595: else if ( m2 )
596: mr0 = m2;
597: else {
598: *pr = 0;
599: return;
600: }
601: else if ( m1 )
602: NEXT(mr) = m1;
603: else if ( m2 )
604: NEXT(mr) = m2;
605: else
606: NEXT(mr) = 0;
1.3 noro 607: MKDP(NV(p1),mr0,*pr);
1.1 noro 608: if ( *pr )
609: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
610: }
611: }
612:
613: void _submd(vl,mod,p1,p2,pr)
614: VL vl;
615: int mod;
616: DP p1,p2,*pr;
617: {
618: DP t;
619:
620: if ( !p2 )
621: *pr = p1;
622: else {
623: _chsgnmd(mod,p2,&t); _addmd(vl,mod,p1,t,pr);
624: }
625: }
626:
627: void _chsgnmd(mod,p,pr)
628: int mod;
629: DP p,*pr;
630: {
631: MP m,mr,mr0;
632:
633: if ( !p )
634: *pr = 0;
635: else {
636: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
637: NEXTMP(mr0,mr); C(mr) = STOI(mod - ITOS(C(m))); mr->dl = m->dl;
638: }
1.3 noro 639: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1.1 noro 640: if ( *pr )
641: (*pr)->sugar = p->sugar;
642: }
643: }
644:
645: void _mulmd(vl,mod,p1,p2,pr)
646: VL vl;
647: int mod;
648: DP p1,p2,*pr;
649: {
1.2 noro 650: if ( !do_weyl )
651: _comm_mulmd(vl,mod,p1,p2,pr);
652: else
653: _weyl_mulmd(vl,mod,p1,p2,pr);
654: }
655:
656: void _comm_mulmd(vl,mod,p1,p2,pr)
657: VL vl;
658: int mod;
659: DP p1,p2,*pr;
660: {
1.1 noro 661: MP m;
662: DP s,t,u;
1.2 noro 663: int i,l,l1;
664: static MP *w;
665: static int wlen;
1.1 noro 666:
667: if ( !p1 || !p2 )
668: *pr = 0;
669: else {
1.2 noro 670: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
671: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
672: if ( l1 < l ) {
673: t = p1; p1 = p2; p2 = t;
674: l = l1;
675: }
676: if ( l > wlen ) {
677: if ( w ) GC_free(w);
678: w = (MP *)MALLOC(l*sizeof(MP));
679: wlen = l;
680: }
681: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
682: w[i] = m;
683: for ( s = 0, i = l-1; i >= 0; i-- ) {
684: _mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
1.1 noro 685: }
1.2 noro 686: bzero(w,l*sizeof(MP));
687: *pr = s;
688: }
689: }
690:
691: void _weyl_mulmd(vl,mod,p1,p2,pr)
692: VL vl;
693: int mod;
694: DP p1,p2,*pr;
695: {
696: MP m;
697: DP s,t,u;
698: int i,l,l1;
699: static MP *w;
700: static int wlen;
701:
702: if ( !p1 || !p2 )
703: *pr = 0;
704: else {
705: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
706: if ( l > wlen ) {
707: if ( w ) GC_free(w);
708: w = (MP *)MALLOC(l*sizeof(MP));
709: wlen = l;
710: }
711: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
712: w[i] = m;
713: for ( s = 0, i = l-1; i >= 0; i-- ) {
714: _weyl_mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
715: }
716: bzero(w,l*sizeof(MP));
1.1 noro 717: *pr = s;
718: }
719: }
720:
721: void _mulmdm(vl,mod,p,m0,pr)
722: VL vl;
723: int mod;
724: DP p;
725: MP m0;
726: DP *pr;
727: {
728: MP m,mr,mr0;
729: DL d;
730: int c,n,r;
731:
732: if ( !p )
733: *pr = 0;
734: else {
735: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
736: m; m = NEXT(m) ) {
737: NEXTMP(mr0,mr);
738: C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
739: adddl(n,m->dl,d,&mr->dl);
740: }
1.3 noro 741: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1.1 noro 742: if ( *pr )
743: (*pr)->sugar = p->sugar + m0->dl->td;
744: }
745: }
746:
1.2 noro 747: void _weyl_mulmdm(vl,mod,p,m0,pr)
748: VL vl;
749: int mod;
750: DP p;
751: MP m0;
752: DP *pr;
753: {
754: DP r,t,t1;
755: MP m;
756: int n,l,i;
757: static MP *w;
758: static int wlen;
759:
760: if ( !p )
761: *pr = 0;
762: else {
763: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
764: if ( l > wlen ) {
765: if ( w ) GC_free(w);
766: w = (MP *)MALLOC(l*sizeof(MP));
767: wlen = l;
768: }
769: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
770: w[i] = m;
771: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
772: _weyl_mulmmm(vl,mod,w[i],m0,n,&t);
773: _addmd(vl,mod,r,t,&t1); r = t1;
774: }
775: bzero(w,l*sizeof(MP));
776: if ( r )
777: r->sugar = p->sugar + m0->dl->td;
778: *pr = r;
779: }
780: }
781:
782: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
783:
784: void _weyl_mulmmm(vl,mod,m0,m1,n,pr)
785: VL vl;
786: int mod;
787: MP m0,m1;
788: int n;
789: DP *pr;
790: {
791: MP m,mr,mr0;
792: DP r,t,t1;
793: int c,c0,c1,cc;
794: DL d,d0,d1;
795: int i,j,a,b,k,l,n2,s,min,h;
796: static int *tab;
797: static int tablen;
798:
799: if ( !m0 || !m1 )
800: *pr = 0;
801: else {
802: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
803: c = dmar(c0,c1,0,mod);
804: d0 = m0->dl; d1 = m1->dl;
805: n2 = n>>1;
806:
1.3 noro 807: NEWDL(d,n);
808: if ( n & 1 )
1.2 noro 809: /* offset of h-degree */
1.3 noro 810: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
811: else
812: d->td = 0;
1.4 noro 813: NEWMP(mr); mr->c = STOI(c); mr->dl = d; NEXT(mr) = 0;
1.3 noro 814: MKDP(n,mr,r); r->sugar = d->td;
815:
816: /* homogenized computation; dx-xd=h^2 */
817: for ( i = 0; i < n2; i++ ) {
818: a = d0->d[i]; b = d1->d[n2+i];
819: k = d0->d[n2+i]; l = d1->d[i];
820: /* degree of xi^a*(Di^k*xi^l)*Di^b */
821: s = a+k+l+b;
822: /* compute xi^a*(Di^k*xi^l)*Di^b */
823: min = MIN(k,l);
824:
825: if ( min+1 > tablen ) {
826: if ( tab ) GC_free(tab);
827: tab = (int *)MALLOC((min+1)*sizeof(int));
828: tablen = min+1;
829: }
830: mkwcm(k,l,mod,tab);
831: if ( n & 1 )
1.2 noro 832: for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3 noro 833: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 834: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.3 noro 835: d->td = s;
836: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
837: mr->c = STOI(tab[j]); mr->dl = d;
1.2 noro 838: }
1.3 noro 839: else
1.2 noro 840: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3 noro 841: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 842: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
843: d->td = d->d[i]+d->d[n2+i]; /* XXX */
844: s = MAX(s,d->td); /* XXX */
1.3 noro 845: mr->c = STOI(tab[j]); mr->dl = d;
1.2 noro 846: }
1.3 noro 847: bzero(tab,(min+1)*sizeof(int));
848: if ( mr0 )
849: NEXT(mr) = 0;
850: MKDP(n,mr0,t);
851: if ( t )
852: t->sugar = s;
853: _comm_mulmd(vl,mod,r,t,&t1); r = t1;
1.2 noro 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: }