![[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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.6 2000/08/22 05:04:05 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();
60: void _comm_mulmd();
61: void _weyl_mulmd();
62: void _weyl_mulmmm();
63: void _weyl_mulmdm();
1.1 noro 64:
65: void ptomd(vl,mod,dvl,p,pr)
66: VL vl,dvl;
67: int mod;
68: P p;
69: DP *pr;
70: {
71: P t;
72:
73: ptomp(mod,p,&t);
74: mptomd(vl,mod,dvl,t,pr);
75: }
76:
77: void mptomd(vl,mod,dvl,p,pr)
78: VL vl,dvl;
79: int mod;
80: P p;
81: DP *pr;
82: {
83: int n,i;
84: VL tvl;
85: V v;
86: DL d;
87: MP m;
88: DCP dc;
89: DP r,s,t,u;
90: P x,c;
91:
92: if ( !p )
93: *pr = 0;
94: else {
95: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
96: if ( NUM(p) ) {
97: NEWDL(d,n);
98: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr);
99: } else {
100: for ( i = 0, tvl = dvl, v = VR(p);
101: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
102: if ( !tvl ) {
103: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
104: mptomd(vl,mod,dvl,COEF(dc),&t); pwrmp(vl,mod,x,DEG(dc),&c);
105: mulmdc(vl,mod,t,c,&r); addmd(vl,mod,r,s,&t); s = t;
106: }
107: *pr = s;
108: } else {
109: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
110: mptomd(vl,mod,dvl,COEF(dc),&t);
111: NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
112: NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u);
1.2 noro 113: comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t;
1.1 noro 114: }
115: *pr = s;
116: }
117: }
118: }
119: }
120:
121: void mdtop(vl,mod,dvl,p,pr)
122: VL vl,dvl;
123: int mod;
124: DP p;
125: P *pr;
126: {
127: int n,i;
128: DL d;
129: MP m;
130: P r,s,t,u,w;
131: Q q;
132: VL tvl;
133:
134: if ( !p )
135: *pr = 0;
136: else {
137: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
138: for ( i = 0, t = C(m), d = m->dl, tvl = dvl;
139: i < n; tvl = NEXT(tvl), i++ ) {
140: MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
141: mulmp(vl,mod,t,u,&w); t = w;
142: }
143: addmp(vl,mod,s,t,&u); s = u;
144: }
145: mptop(s,pr);
146: }
147: }
148:
149: void addmd(vl,mod,p1,p2,pr)
150: VL vl;
151: int mod;
152: DP p1,p2,*pr;
153: {
154: int n;
155: MP m1,m2,mr,mr0;
156: P t;
157: int tmp;
158: MQ q;
159:
160: if ( !p1 )
161: *pr = p2;
162: else if ( !p2 )
163: *pr = p1;
164: else {
165: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
166: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
167: case 0:
168: if ( NUM(C(m1)) && NUM(C(m2)) ) {
169: tmp = (CONT((MQ)C(m1))+CONT((MQ)C(m2))) - mod;
170: if ( tmp < 0 )
171: tmp += mod;
172: if ( tmp ) {
173: STOMQ(tmp,q); t = (P)q;
174: } else
175: t = 0;
176: } else
177: addmp(vl,mod,C(m1),C(m2),&t);
178: if ( t ) {
179: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
180: }
181: m1 = NEXT(m1); m2 = NEXT(m2); break;
182: case 1:
183: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
184: m1 = NEXT(m1); break;
185: case -1:
186: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
187: m2 = NEXT(m2); break;
188: }
189: if ( !mr0 )
190: if ( m1 )
191: mr0 = m1;
192: else if ( m2 )
193: mr0 = m2;
194: else {
195: *pr = 0;
196: return;
197: }
198: else if ( m1 )
199: NEXT(mr) = m1;
200: else if ( m2 )
201: NEXT(mr) = m2;
202: else
203: NEXT(mr) = 0;
204: MKDP(NV(p1),mr0,*pr);
205: if ( *pr )
206: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
207: }
208: }
209:
210: void submd(vl,mod,p1,p2,pr)
211: VL vl;
212: int mod;
213: DP p1,p2,*pr;
214: {
215: DP t;
216:
217: if ( !p2 )
218: *pr = p1;
219: else {
220: chsgnmd(mod,p2,&t); addmd(vl,mod,p1,t,pr);
221: }
222: }
223:
224: void chsgnmd(mod,p,pr)
225: int mod;
226: DP p,*pr;
227: {
228: MP m,mr,mr0;
229:
230: if ( !p )
231: *pr = 0;
232: else {
233: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
234: NEXTMP(mr0,mr); chsgnmp(mod,C(m),&C(mr)); mr->dl = m->dl;
235: }
236: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
237: if ( *pr )
238: (*pr)->sugar = p->sugar;
239: }
240: }
241:
242: void mulmd(vl,mod,p1,p2,pr)
243: VL vl;
244: int mod;
245: DP p1,p2,*pr;
246: {
1.2 noro 247: if ( !do_weyl )
248: comm_mulmd(vl,mod,p1,p2,pr);
249: else
250: weyl_mulmd(vl,mod,p1,p2,pr);
251: }
252:
253: void comm_mulmd(vl,mod,p1,p2,pr)
254: VL vl;
255: int mod;
256: DP p1,p2,*pr;
257: {
258: MP m;
259: DP s,t,u;
260: int i,l,l1;
261: static MP *w;
262: static int wlen;
263:
264: if ( !p1 || !p2 )
265: *pr = 0;
266: else if ( OID(p1) <= O_P )
267: mulmdc(vl,mod,p2,(P)p1,pr);
268: else if ( OID(p2) <= O_P )
269: mulmdc(vl,mod,p1,(P)p2,pr);
270: else {
271: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
272: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
273: if ( l1 < l ) {
274: t = p1; p1 = p2; p2 = t;
275: l = l1;
276: }
277: if ( l > wlen ) {
278: if ( w ) GC_free(w);
279: w = (MP *)MALLOC(l*sizeof(MP));
280: wlen = l;
281: }
282: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
283: w[i] = m;
284: for ( s = 0, i = l-1; i >= 0; i-- ) {
285: mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
286: }
287: bzero(w,l*sizeof(MP));
288: *pr = s;
289: }
290: }
291:
292: void weyl_mulmd(vl,mod,p1,p2,pr)
293: VL vl;
294: int mod;
295: DP p1,p2,*pr;
296: {
1.1 noro 297: MP m;
298: DP s,t,u;
1.2 noro 299: int i,l,l1;
300: static MP *w;
301: static int wlen;
1.1 noro 302:
303: if ( !p1 || !p2 )
304: *pr = 0;
305: else if ( OID(p1) <= O_P )
306: mulmdc(vl,mod,p2,(P)p1,pr);
307: else if ( OID(p2) <= O_P )
308: mulmdc(vl,mod,p1,(P)p2,pr);
309: else {
1.2 noro 310: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
311: if ( l > wlen ) {
312: if ( w ) GC_free(w);
313: w = (MP *)MALLOC(l*sizeof(MP));
314: wlen = l;
315: }
316: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
317: w[i] = m;
318: for ( s = 0, i = l-1; i >= 0; i-- ) {
319: weyl_mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
1.1 noro 320: }
1.2 noro 321: bzero(w,l*sizeof(MP));
1.1 noro 322: *pr = s;
323: }
324: }
325:
326: void mulmdm(vl,mod,p,m0,pr)
327: VL vl;
328: int mod;
329: DP p;
330: MP m0;
331: DP *pr;
332: {
333: MP m,mr,mr0;
334: P c;
335: DL d;
336: int n,t;
337: MQ q;
338:
339: if ( !p )
340: *pr = 0;
341: else {
342: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
343: m; m = NEXT(m) ) {
344: NEXTMP(mr0,mr);
345: if ( NUM(C(m)) && NUM(c) ) {
346: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
347: STOMQ(t,q); C(mr) = (P)q;
348: } else
349: mulmp(vl,mod,C(m),c,&C(mr));
350: adddl(n,m->dl,d,&mr->dl);
351: }
352: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
353: if ( *pr )
354: (*pr)->sugar = p->sugar + m0->dl->td;
355: }
356: }
357:
1.2 noro 358: void weyl_mulmdm(vl,mod,p,m0,pr)
359: VL vl;
360: int mod;
361: DP p;
362: MP m0;
363: DP *pr;
364: {
365: DP r,t,t1;
366: MP m;
367: int n,l,i;
368: static MP *w;
369: static int wlen;
370:
371: if ( !p )
372: *pr = 0;
373: else {
374: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
375: if ( l > wlen ) {
376: if ( w ) GC_free(w);
377: w = (MP *)MALLOC(l*sizeof(MP));
378: wlen = l;
379: }
380: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
381: w[i] = m;
382: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
383: weyl_mulmmm(vl,mod,w[i],m0,n,&t);
384: addmd(vl,mod,r,t,&t1); r = t1;
385: }
386: bzero(w,l*sizeof(MP));
387: if ( r )
388: r->sugar = p->sugar + m0->dl->td;
389: *pr = r;
390: }
391: }
392:
393: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
394:
395: void weyl_mulmmm(vl,mod,m0,m1,n,pr)
396: VL vl;
397: int mod;
398: MP m0,m1;
399: int n;
400: DP *pr;
401: {
402: MP m,mr,mr0;
403: MQ mq;
404: DP r,t,t1;
405: P c,c0,c1,cc;
406: DL d,d0,d1;
1.3 noro 407: int i,j,a,b,k,l,n2,s,min;
1.2 noro 408: static int *tab;
409: static int tablen;
410:
411: if ( !m0 || !m1 )
412: *pr = 0;
413: else {
414: c0 = C(m0); c1 = C(m1);
415: mulmp(vl,mod,c0,c1,&c);
416: d0 = m0->dl; d1 = m1->dl;
417: n2 = n>>1;
418: if ( n & 1 ) {
419: /* homogenized computation; dx-xd=h^2 */
420: /* offset of h-degree */
1.3 noro 421: NEWDL(d,n);
422: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1.4 noro 423: NEWMP(mr); mr->c = (P)ONEM; mr->dl = d; NEXT(mr) = 0;
1.3 noro 424: MKDP(n,mr,r); r->sugar = d->td;
425: } else
426: r = (DP)ONEM;
427: for ( i = 0; i < n2; i++ ) {
428: a = d0->d[i]; b = d1->d[n2+i];
429: k = d0->d[n2+i]; l = d1->d[i];
430: /* degree of xi^a*(Di^k*xi^l)*Di^b */
431: s = a+k+l+b;
432: /* compute xi^a*(Di^k*xi^l)*Di^b */
433: min = MIN(k,l);
434:
435: if ( min+1 > tablen ) {
436: if ( tab ) GC_free(tab);
437: tab = (int *)MALLOC((min+1)*sizeof(int));
438: tablen = min+1;
439: }
440: mkwcm(k,l,mod,tab);
441: if ( n & 1 )
1.2 noro 442: for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3 noro 443: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 444: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.3 noro 445: d->td = s;
446: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
447: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 448: }
1.3 noro 449: else
1.2 noro 450: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3 noro 451: NEXTMP(mr0,mr); NEWDL(d,n);
1.2 noro 452: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
453: d->td = d->d[i]+d->d[n2+i]; /* XXX */
454: s = MAX(s,d->td); /* XXX */
1.3 noro 455: STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2 noro 456: }
1.3 noro 457: bzero(tab,(min+1)*sizeof(int));
458: if ( mr0 )
459: NEXT(mr) = 0;
460: MKDP(n,mr0,t);
461: if ( t )
462: t->sugar = s;
463: comm_mulmd(vl,mod,r,t,&t1); r = t1;
464: }
1.2 noro 465: mulmdc(vl,mod,r,c,pr);
466: }
467: }
468:
1.1 noro 469: void mulmdc(vl,mod,p,c,pr)
470: VL vl;
471: int mod;
472: DP p;
473: P c;
474: DP *pr;
475: {
476: MP m,mr,mr0;
477: int t;
478: MQ q;
479:
480: if ( !p || !c )
481: *pr = 0;
482: else {
483: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
484: NEXTMP(mr0,mr);
485: if ( NUM(C(m)) && NUM(c) ) {
486: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
487: STOMQ(t,q); C(mr) = (P)q;
488: } else
489: mulmp(vl,mod,C(m),c,&C(mr));
490: mr->dl = m->dl;
491: }
492: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
493: if ( *pr )
494: (*pr)->sugar = p->sugar;
495: }
496: }
497:
498: void divsmdc(vl,mod,p,c,pr)
499: VL vl;
500: int mod;
501: DP p;
502: P c;
503: DP *pr;
504: {
505: MP m,mr,mr0;
506:
507: if ( !c )
508: error("disvsdc : division by 0");
509: else if ( !p )
510: *pr = 0;
511: else {
512: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
513: NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl;
514: }
515: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
516: if ( *pr )
517: (*pr)->sugar = p->sugar;
518: }
519: }
520:
521: void _dtop_mod(vl,dvl,p,pr)
522: VL vl,dvl;
523: DP p;
524: P *pr;
525: {
526: int n,i;
527: DL d;
528: MP m;
529: P r,s,t,u,w;
530: Q q;
531: VL tvl;
532:
533: if ( !p )
534: *pr = 0;
535: else {
536: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
537: i = ITOS(m->c); STOQ(i,q); t = (P)q;
538: for ( i = 0, d = m->dl, tvl = dvl;
539: i < n; tvl = NEXT(tvl), i++ ) {
540: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
541: mulp(vl,t,u,&w); t = w;
542: }
543: addp(vl,s,t,&u); s = u;
544: }
545: *pr = s;
546: }
547: }
548:
549: void _dp_mod(p,mod,subst,rp)
550: DP p;
551: int mod;
552: NODE subst;
553: DP *rp;
554: {
555: MP m,mr,mr0;
556: P t,s;
557: NODE tn;
558:
559: if ( !p )
560: *rp = 0;
561: else {
562: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
563: for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) {
564: substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t);
565: s = t;
566: }
567: ptomp(mod,s,&t);
568: if ( t ) {
569: NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl;
570: }
571: }
572: if ( mr0 ) {
573: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
574: } else
575: *rp = 0;
576: }
577: }
578:
1.2 noro 579: void _dp_monic(p,mod,rp)
580: DP p;
581: int mod;
582: DP *rp;
583: {
584: MP m,mr,mr0;
585: int c,c1;
586: NODE tn;
587:
588: if ( !p )
589: *rp = 0;
590: else {
591: c = invm(ITOS(BDY(p)->c),mod);
592: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
593: c1 = dmar(ITOS(m->c),c,0,mod);
594: NEXTMP(mr0,mr); mr->c = STOI(c1); mr->dl = m->dl;
595: }
596: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
597: }
1.1 noro 598: }
599:
600: void _printdp(d)
601: DP d;
602: {
603: int n,i;
604: MP m;
605: DL dl;
606:
607: if ( !d ) {
608: printf("0"); return;
609: }
610: for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) {
611: printf("%d*<<",ITOS(m->c));
612: for ( i = 0, dl = m->dl; i < n-1; i++ )
613: printf("%d,",dl->d[i]);
614: printf("%d",dl->d[i]);
615: printf(">>");
616: if ( NEXT(m) )
617: printf("+");
618: }
619: }