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