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