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