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