![[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.2 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.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 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 2001/03/19 04:02:03 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "inline.h"
52:
53: extern int (*cmpdl)();
54: extern int do_weyl;
55:
1.5 noro 56: void dpto_dp();
1.1 noro 57: void _dptodp();
1.5 noro 58: void _adddl_dup();
1.4 noro 59: void adddl_destructive();
1.1 noro 60: void _mulmdm_dup();
61: void _free_dp();
62: void _comm_mulmd_dup();
63: void _weyl_mulmd_dup();
64: void _weyl_mulmmm_dup();
65: void _weyl_mulmdm_dup();
1.4 noro 66: void _comm_mulmd_tab();
67: void _comm_mulmd_tab_destructive();
1.1 noro 68:
69: MP _mp_free_list;
70: DP _dp_free_list;
71: DL _dl_free_list;
72: int current_dl_length;
1.6 noro 73:
74: void _free_private_storage()
75: {
76: _mp_free_list = 0;
77: _dp_free_list = 0;
78: _dl_free_list = 0;
79: GC_gcollect();
80: }
1.1 noro 81:
82: void _DL_alloc()
83: {
84: int *p;
85: int i,dl_len;
86: static int DL_alloc_count;
87:
88: /* fprintf(stderr,"DL_alloc : %d \n",++DL_alloc_count); */
89: dl_len = (current_dl_length+1);
1.8 noro 90: #if defined(LONG_IS_64BIT)
91: if ( dl_len & 1 )
92: dl_len += 1;
93: #endif
1.1 noro 94: p = (int *)GC_malloc(128*dl_len*sizeof(int));
95: for ( i = 0; i < 128; i++, p += dl_len ) {
96: *(DL *)p = _dl_free_list;
97: _dl_free_list = (DL)p;
98: }
99: }
100:
101: void _MP_alloc()
102: {
103: MP p;
104: int i;
105: static int MP_alloc_count;
106:
107: /* fprintf(stderr,"MP_alloc : %d \n",++MP_alloc_count); */
108: p = (MP)GC_malloc(1024*sizeof(struct oMP));
109: for ( i = 0; i < 1024; i++ ) {
110: p[i].next = _mp_free_list; _mp_free_list = &p[i];
111: }
112: }
113:
114: void _DP_alloc()
115: {
116: DP p;
117: int i;
118: static int DP_alloc_count;
119:
120: /* fprintf(stderr,"DP_alloc : %d \n",++DP_alloc_count); */
121: p = (DP)GC_malloc(1024*sizeof(struct oDP));
122: for ( i = 0; i < 1024; i++ ) {
123: p[i].body = (MP)_dp_free_list; _dp_free_list = &p[i];
124: }
125: }
126:
127: /* merge p1 and p2 into pr */
128:
129: void _addmd_destructive(mod,p1,p2,pr)
130: int mod;
131: DP p1,p2,*pr;
132: {
133: int n;
134: MP m1,m2,mr,mr0,s;
135: int t;
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: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
146: if ( t < 0 )
147: t += mod;
148: s = m1; m1 = NEXT(m1);
149: if ( t ) {
150: _NEXTMP2(mr0,mr,s); C(mr) = STOI(t);
151: } else {
1.5 noro 152: _FREEDL(s->dl); _FREEMP(s);
1.1 noro 153: }
1.5 noro 154: s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
1.1 noro 155: break;
156: case 1:
157: s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
158: break;
159: case -1:
160: s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
161: break;
162: }
163: if ( !mr0 )
164: if ( m1 )
165: mr0 = m1;
166: else if ( m2 )
167: mr0 = m2;
168: else {
169: *pr = 0;
170: return;
171: }
172: else if ( m1 )
173: NEXT(mr) = m1;
174: else if ( m2 )
175: NEXT(mr) = m2;
176: else
177: NEXT(mr) = 0;
178: _MKDP(NV(p1),mr0,*pr);
179: if ( *pr )
180: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.5 noro 181: _FREEDP(p1); _FREEDP(p2);
1.1 noro 182: }
183: }
184:
185: void _mulmd_dup(mod,p1,p2,pr)
186: int mod;
187: DP p1,p2,*pr;
188: {
189: if ( !do_weyl )
190: _comm_mulmd_dup(mod,p1,p2,pr);
191: else
192: _weyl_mulmd_dup(mod,p1,p2,pr);
193: }
194:
195: void _comm_mulmd_dup(mod,p1,p2,pr)
196: int mod;
197: DP p1,p2,*pr;
198: {
199: MP m;
200: DP s,t,u;
201: int i,l,l1;
202: static MP *w;
203: static int wlen;
204:
205: if ( !p1 || !p2 )
206: *pr = 0;
207: else {
208: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
209: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
210: if ( l1 < l ) {
211: t = p1; p1 = p2; p2 = t;
212: l = l1;
213: }
214: if ( l > wlen ) {
215: if ( w ) GC_free(w);
216: w = (MP *)MALLOC(l*sizeof(MP));
217: wlen = l;
218: }
219: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
220: w[i] = m;
221: for ( s = 0, i = l-1; i >= 0; i-- ) {
222: _mulmdm_dup(mod,p1,w[i],&t); _addmd_destructive(mod,s,t,&u); s = u;
223: }
224: bzero(w,l*sizeof(MP));
225: *pr = s;
226: }
227: }
228:
229: void _weyl_mulmd_dup(mod,p1,p2,pr)
230: int mod;
231: DP p1,p2,*pr;
232: {
233: MP m;
234: DP s,t,u;
235: int i,l,l1;
236: static MP *w;
237: static int wlen;
238:
239: if ( !p1 || !p2 )
240: *pr = 0;
241: else {
1.4 noro 242: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.1 noro 243: if ( l > wlen ) {
244: if ( w ) GC_free(w);
245: w = (MP *)MALLOC(l*sizeof(MP));
246: wlen = l;
247: }
1.4 noro 248: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.1 noro 249: w[i] = m;
250: for ( s = 0, i = l-1; i >= 0; i-- ) {
1.4 noro 251: _weyl_mulmdm_dup(mod,w[i],p2,&t); _addmd_destructive(mod,s,t,&u); s = u;
1.1 noro 252: }
253: bzero(w,l*sizeof(MP));
254: *pr = s;
255: }
256: }
257:
258: void _mulmdm_dup(mod,p,m0,pr)
259: int mod;
260: DP p;
261: MP m0;
262: DP *pr;
263: {
264: MP m,mr,mr0;
265: DL d,dt,dm;
266: int c,n,r,i;
267: int *pt,*p1,*p2;
268:
269: if ( !p )
270: *pr = 0;
271: else {
272: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
273: m; m = NEXT(m) ) {
274: _NEXTMP(mr0,mr);
275: C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
276: _NEWDL_NOINIT(dt,n); mr->dl = dt;
277: dm = m->dl;
278: dt->td = d->td + dm->td;
279: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
280: *pt++ = *p1++ + *p2++;
281: }
282: NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
283: if ( *pr )
284: (*pr)->sugar = p->sugar + m0->dl->td;
285: }
286: }
287:
1.4 noro 288: void _weyl_mulmmm_dup_bug();
289:
290: void _weyl_mulmdm_dup(mod,m0,p,pr)
1.1 noro 291: int mod;
1.4 noro 292: MP m0;
1.1 noro 293: DP p;
294: DP *pr;
295: {
296: DP r,t,t1;
297: MP m;
1.4 noro 298: DL d0;
299: int n,n2,l,i,j,tlen;
300: static MP *w,*psum;
301: static struct cdlm *tab;
1.1 noro 302: static int wlen;
1.4 noro 303: static int rtlen;
1.1 noro 304:
305: if ( !p )
306: *pr = 0;
307: else {
308: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
309: if ( l > wlen ) {
310: if ( w ) GC_free(w);
311: w = (MP *)MALLOC(l*sizeof(MP));
312: wlen = l;
313: }
314: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
315: w[i] = m;
1.4 noro 316: n = NV(p); n2 = n>>1;
317: d0 = m0->dl;
318:
319: for ( i = 0, tlen = 1; i < n2; i++ )
320: tlen *= d0->d[n2+i]+1;
321: if ( tlen > rtlen ) {
322: if ( tab ) GC_free(tab);
323: if ( psum ) GC_free(psum);
324: rtlen = tlen;
325: tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
326: psum = (MP *)MALLOC(rtlen*sizeof(MP));
1.1 noro 327: }
1.4 noro 328: bzero(psum,tlen*sizeof(MP));
329: for ( i = l-1; i >= 0; i-- ) {
330: bzero(tab,tlen*sizeof(struct cdlm));
331: _weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
332: for ( j = 0; j < tlen; j++ ) {
333: if ( tab[j].c ) {
334: _NEWMP(m); m->dl = tab[j].d;
335: C(m) = STOI(tab[j].c); NEXT(m) = psum[j];
336: psum[j] = m;
337: }
338: }
339: }
340: for ( j = tlen-1, r = 0; j >= 0; j-- )
341: if ( psum[j] ) {
342: _MKDP(n,psum[j],t); _addmd_destructive(mod,r,t,&t1); r = t1;
343: }
1.1 noro 344: if ( r )
345: r->sugar = p->sugar + m0->dl->td;
346: *pr = r;
347: }
348: }
1.4 noro 349:
1.1 noro 350: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
351:
1.4 noro 352: void _weyl_mulmmm_dup(mod,m0,m1,n,rtab,rtablen)
1.1 noro 353: int mod;
354: MP m0,m1;
355: int n;
1.4 noro 356: struct cdlm *rtab;
357: int rtablen;
1.1 noro 358: {
359: MP m,mr,mr0;
360: DP r,t,t1;
361: int c,c0,c1,cc;
1.4 noro 362: DL d,d0,d1,dt;
363: int i,j,a,b,k,l,n2,s,min,h,curlen;
364: struct cdlm *p;
365: static int *ctab;
366: static struct cdlm *tab;
1.1 noro 367: static int tablen;
1.4 noro 368: static struct cdlm *tmptab;
369: static int tmptablen;
1.1 noro 370:
1.4 noro 371: if ( !m0 || !m1 ) {
372: rtab[0].c = 0;
373: rtab[0].d = 0;
374: return;
375: }
376: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
377: c = dmar(c0,c1,0,mod);
378: d0 = m0->dl; d1 = m1->dl;
379: n2 = n>>1;
380: curlen = 1;
381:
382: _NEWDL(d,n);
383: if ( n & 1 )
384: /* offset of h-degree */
385: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
386: else
387: d->td = 0;
388: rtab[0].c = c;
389: rtab[0].d = d;
390:
391: if ( rtablen > tmptablen ) {
392: if ( tmptab ) GC_free(tmptab);
393: tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
394: tmptablen = rtablen;
395: }
1.1 noro 396:
1.4 noro 397: for ( i = 0; i < n2; i++ ) {
398: a = d0->d[i]; b = d1->d[n2+i];
399: k = d0->d[n2+i]; l = d1->d[i];
400: if ( !k || !l ) {
401: a += l;
402: b += k;
403: s = a+b;
404: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
405: if ( p->c ) {
406: dt = p->d;
407: dt->d[i] = a;
408: dt->d[n2+i] = b;
409: dt->td += s;
410: }
411: }
412: curlen *= k+1;
413: continue;
414: }
415: if ( k+1 > tablen ) {
416: if ( tab ) GC_free(tab);
417: if ( ctab ) GC_free(ctab);
418: tablen = k+1;
419: tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
420: ctab = (int *)MALLOC(tablen*sizeof(int));
421: }
422: /* degree of xi^a*(Di^k*xi^l)*Di^b */
423: s = a+k+l+b;
424: /* compute xi^a*(Di^k*xi^l)*Di^b */
425: min = MIN(k,l);
426: mkwcm(k,l,mod,ctab);
427: bzero(tab,(k+1)*sizeof(struct cdlm));
428: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
1.1 noro 429: if ( n & 1 )
1.4 noro 430: for ( j = 0; j <= min; j++ ) {
431: _NEWDL(d,n);
432: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
433: d->td = s;
434: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
435: tab[j].d = d;
436: tab[j].c = ctab[j];
437: }
1.1 noro 438: else
1.4 noro 439: for ( j = 0; j <= min; j++ ) {
440: _NEWDL(d,n);
441: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
442: d->td = d->d[i]+d->d[n2+i]; /* XXX */
443: tab[j].d = d;
444: tab[j].c = ctab[j];
445: }
446: #if 0
447: _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
448: for ( j = 0; j < curlen; j++ )
1.5 noro 449: if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
1.4 noro 450: for ( j = 0; j <= min; j++ )
1.5 noro 451: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4 noro 452: curlen *= k+1;
453: bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
454: #else
455: _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
456: for ( j = 0; j <= min; j++ )
1.5 noro 457: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4 noro 458: curlen *= k+1;
459: #endif
460: }
461: }
462:
463: /* direct product of two cdlm tables
464: rt[] = [
465: t[0]*t1[0],...,t[n-1]*t1[0],
466: t[0]*t1[1],...,t[n-1]*t1[1],
467: ...
468: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
469: ]
470: */
471:
472: void _comm_mulmd_tab(mod,nv,t,n,t1,n1,rt)
473: int mod;
474: int nv;
475: struct cdlm *t;
476: int n;
477: struct cdlm *t1;
478: int n1;
479: struct cdlm *rt;
480: {
481: int i,j;
482: struct cdlm *p;
483: int c;
484: DL d;
485:
486: bzero(rt,n*n1*sizeof(struct cdlm));
487: for ( j = 0, p = rt; j < n1; j++ ) {
488: c = t1[j].c;
489: d = t1[j].d;
490: if ( !c )
491: break;
492: for ( i = 0; i < n; i++, p++ ) {
493: if ( t[i].c ) {
494: p->c = dmar(t[i].c,c,0,mod);
1.5 noro 495: _adddl_dup(nv,t[i].d,d,&p->d);
1.4 noro 496: }
497: }
498: }
499: }
500:
501: void _comm_mulmd_tab_destructive(mod,nv,t,n,t1,n1)
502: int mod;
503: int nv;
504: struct cdlm *t;
505: int n;
506: struct cdlm *t1;
507: int n1;
508: {
509: int i,j;
510: struct cdlm *p;
511: int c;
512: DL d;
513:
514: for ( j = 1, p = t+n; j < n1; j++ ) {
515: c = t1[j].c;
516: d = t1[j].d;
517: if ( !c )
518: break;
519: for ( i = 0; i < n; i++, p++ ) {
520: if ( t[i].c ) {
521: p->c = dmar(t[i].c,c,0,mod);
1.5 noro 522: _adddl_dup(nv,t[i].d,d,&p->d);
1.1 noro 523: }
524: }
525: }
1.4 noro 526: c = t1[0].c;
527: d = t1[0].d;
528: for ( i = 0, p = t; i < n; i++, p++ )
529: if ( t[i].c ) {
530: p->c = dmar(t[i].c,c,0,mod);
531: /* t[i].d += d */
532: adddl_destructive(nv,t[i].d,d);
533: }
1.1 noro 534: }
535:
1.5 noro 536: void dlto_dl(d,dr)
537: DL d;
538: DL *dr;
1.1 noro 539: {
540: int i,n;
1.5 noro 541: DL t;
1.1 noro 542:
1.5 noro 543: n = current_dl_length;
544: _NEWDL(t,n); *dr = t;
545: t->td = d->td;
1.1 noro 546: for ( i = 0; i < n; i++ )
1.5 noro 547: t->d[i] = d->d[i];
1.1 noro 548: }
549:
550: void _dltodl(d,dr)
551: DL d;
552: DL *dr;
553: {
554: int i,n;
555: DL t;
556:
557: n = current_dl_length;
558: NEWDL(t,n); *dr = t;
559: t->td = d->td;
560: for ( i = 0; i < n; i++ )
561: t->d[i] = d->d[i];
562: }
563:
1.5 noro 564: void _adddl_dup(n,d1,d2,dr)
1.1 noro 565: int n;
566: DL d1,d2;
567: DL *dr;
568: {
569: DL dt;
570: int i;
571:
1.4 noro 572: _NEWDL(dt,n);
573: *dr = dt;
1.1 noro 574: dt->td = d1->td + d2->td;
575: for ( i = 0; i < n; i++ )
576: dt->d[i] = d1->d[i]+d2->d[i];
1.4 noro 577: }
578:
579: void _free_dlarray(a,n)
580: DL *a;
581: int n;
582: {
583: int i;
584:
1.5 noro 585: for ( i = 0; i < n; i++ ) { _FREEDL(a[i]); }
1.1 noro 586: }
587:
588: void _free_dp(f)
589: DP f;
590: {
591: MP m,s;
592:
593: if ( !f )
594: return;
595: m = f->body;
596: while ( m ) {
1.5 noro 597: s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
598: }
599: _FREEDP(f);
600: }
601:
602: void dpto_dp(p,r)
603: DP p;
604: DP *r;
605: {
606: MP m,mr0,mr;
1.7 noro 607: DL t;
1.5 noro 608:
609: if ( !p )
610: *r = 0;
611: else {
1.7 noro 612: /* XXX : dummy call to set current_dl_length */
613: _NEWDL_NOINIT(t,NV(p));
614:
1.5 noro 615: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
616: _NEXTMP(mr0,mr);
617: dlto_dl(m->dl,&mr->dl);
618: mr->c = m->c;
619: }
620: NEXT(mr) = 0;
621: _MKDP(p->nv,mr0,*r);
622: (*r)->sugar = p->sugar;
1.1 noro 623: }
624: }
625:
626: void _dptodp(p,r)
627: DP p;
628: DP *r;
629: {
630: MP m,mr0,mr;
631:
632: if ( !p )
633: *r = 0;
634: else {
635: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
636: NEXTMP(mr0,mr);
637: _dltodl(m->dl,&mr->dl);
638: mr->c = m->c;
639: }
640: NEXT(mr) = 0;
641: MKDP(p->nv,mr0,*r);
642: (*r)->sugar = p->sugar;
643: }
644: }
645:
1.5 noro 646: /*
647: * destructive merge of two list
648: *
649: * p1, p2 : list of DL
650: * return : a merged list
651: */
652:
653: NODE _symb_merge(m1,m2,n)
654: NODE m1,m2;
655: int n;
656: {
657: NODE top,prev,cur,m,t;
658:
659: if ( !m1 )
660: return m2;
661: else if ( !m2 )
662: return m1;
663: else {
664: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
665: case 0:
666: top = m1; _FREEDL((DL)BDY(m2)); m = NEXT(m2);
667: break;
668: case 1:
669: top = m1; m = m2;
670: break;
671: case -1:
672: top = m2; m = m1;
1.1 noro 673: break;
674: }
1.5 noro 675: prev = top; cur = NEXT(top);
676: /* BDY(prev) > BDY(m) always holds */
677: while ( cur && m ) {
678: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
679: case 0:
680: _FREEDL(BDY(m)); m = NEXT(m);
681: prev = cur; cur = NEXT(cur);
682: break;
683: case 1:
684: t = NEXT(cur); NEXT(cur) = m; m = t;
685: prev = cur; cur = NEXT(cur);
686: break;
687: case -1:
688: NEXT(prev) = m; m = cur;
689: prev = NEXT(prev); cur = NEXT(prev);
690: break;
1.1 noro 691: }
692: }
1.5 noro 693: if ( !cur )
694: NEXT(prev) = m;
695: return top;
1.1 noro 696: }
1.9 ! noro 697: }
! 698:
! 699: /* XXX : bellow should be placed in another file */
! 700:
! 701: void dpto_dp();
! 702: void _dptodp();
! 703: void _adddl_dup();
! 704: void adddl_destructive();
! 705: void _muldm_dup();
! 706: void _free_dp();
! 707: void _comm_muld_dup();
! 708: void _weyl_muld_dup();
! 709: void _weyl_mulmm_dup();
! 710: void _weyl_muldm_dup();
! 711: void _comm_muld_tab();
! 712: void _comm_muld_tab_destructive();
! 713:
! 714: /* merge p1 and p2 into pr */
! 715:
! 716: void _addd_destructive(vl,p1,p2,pr)
! 717: VL vl;
! 718: DP p1,p2,*pr;
! 719: {
! 720: int n;
! 721: MP m1,m2,mr,mr0,s;
! 722: P t;
! 723:
! 724: if ( !p1 )
! 725: *pr = p2;
! 726: else if ( !p2 )
! 727: *pr = p1;
! 728: else {
! 729: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
! 730: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
! 731: case 0:
! 732: addp(vl,C(m1),C(m2),&t);
! 733: s = m1; m1 = NEXT(m1);
! 734: if ( t ) {
! 735: _NEXTMP2(mr0,mr,s); C(mr) = t;
! 736: } else {
! 737: _FREEDL(s->dl); _FREEMP(s);
! 738: }
! 739: s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
! 740: break;
! 741: case 1:
! 742: s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
! 743: break;
! 744: case -1:
! 745: s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
! 746: break;
! 747: }
! 748: if ( !mr0 )
! 749: if ( m1 )
! 750: mr0 = m1;
! 751: else if ( m2 )
! 752: mr0 = m2;
! 753: else {
! 754: *pr = 0;
! 755: return;
! 756: }
! 757: else if ( m1 )
! 758: NEXT(mr) = m1;
! 759: else if ( m2 )
! 760: NEXT(mr) = m2;
! 761: else
! 762: NEXT(mr) = 0;
! 763: _MKDP(NV(p1),mr0,*pr);
! 764: if ( *pr )
! 765: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 766: _FREEDP(p1); _FREEDP(p2);
! 767: }
! 768: }
! 769:
! 770: void _muld_dup(vl,p1,p2,pr)
! 771: VL vl;
! 772: DP p1,p2,*pr;
! 773: {
! 774: if ( !do_weyl )
! 775: _comm_muld_dup(vl,p1,p2,pr);
! 776: else
! 777: _weyl_muld_dup(vl,p1,p2,pr);
! 778: }
! 779:
! 780: void _comm_muld_dup(vl,p1,p2,pr)
! 781: VL vl;
! 782: DP p1,p2,*pr;
! 783: {
! 784: MP m;
! 785: DP s,t,u;
! 786: int i,l,l1;
! 787: static MP *w;
! 788: static int wlen;
! 789:
! 790: if ( !p1 || !p2 )
! 791: *pr = 0;
! 792: else {
! 793: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
! 794: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 795: if ( l1 < l ) {
! 796: t = p1; p1 = p2; p2 = t;
! 797: l = l1;
! 798: }
! 799: if ( l > wlen ) {
! 800: if ( w ) GC_free(w);
! 801: w = (MP *)MALLOC(l*sizeof(MP));
! 802: wlen = l;
! 803: }
! 804: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 805: w[i] = m;
! 806: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 807: _muldm_dup(vl,p1,w[i],&t); _addd_destructive(vl,s,t,&u); s = u;
! 808: }
! 809: bzero(w,l*sizeof(MP));
! 810: *pr = s;
! 811: }
! 812: }
! 813:
! 814: void _weyl_muld_dup(vl,p1,p2,pr)
! 815: VL vl;
! 816: DP p1,p2,*pr;
! 817: {
! 818: MP m;
! 819: DP s,t,u;
! 820: int i,l,l1;
! 821: static MP *w;
! 822: static int wlen;
! 823:
! 824: if ( !p1 || !p2 )
! 825: *pr = 0;
! 826: else {
! 827: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
! 828: if ( l > wlen ) {
! 829: if ( w ) GC_free(w);
! 830: w = (MP *)MALLOC(l*sizeof(MP));
! 831: wlen = l;
! 832: }
! 833: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
! 834: w[i] = m;
! 835: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 836: _weyl_muldm_dup(vl,w[i],p2,&t); _addd_destructive(vl,s,t,&u); s = u;
! 837: }
! 838: bzero(w,l*sizeof(MP));
! 839: *pr = s;
! 840: }
! 841: }
! 842:
! 843: void _muldm_dup(vl,p,m0,pr)
! 844: VL vl;
! 845: DP p;
! 846: MP m0;
! 847: DP *pr;
! 848: {
! 849: MP m,mr,mr0;
! 850: DL d,dt,dm;
! 851: P c;
! 852: int n,r,i;
! 853: int *pt,*p1,*p2;
! 854:
! 855: if ( !p )
! 856: *pr = 0;
! 857: else {
! 858: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
! 859: m; m = NEXT(m) ) {
! 860: _NEXTMP(mr0,mr);
! 861: mulp(vl,C(m),c,&C(mr));
! 862: _NEWDL_NOINIT(dt,n); mr->dl = dt;
! 863: dm = m->dl;
! 864: dt->td = d->td + dm->td;
! 865: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
! 866: *pt++ = *p1++ + *p2++;
! 867: }
! 868: NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
! 869: if ( *pr )
! 870: (*pr)->sugar = p->sugar + m0->dl->td;
! 871: }
! 872: }
! 873:
! 874: void _weyl_muldm_dup(vl,m0,p,pr)
! 875: VL vl;
! 876: MP m0;
! 877: DP p;
! 878: DP *pr;
! 879: {
! 880: DP r,t,t1;
! 881: MP m;
! 882: DL d0;
! 883: int n,n2,l,i,j,tlen;
! 884: static MP *w,*psum;
! 885: static struct cdl *tab;
! 886: static int wlen;
! 887: static int rtlen;
! 888:
! 889: if ( !p )
! 890: *pr = 0;
! 891: else {
! 892: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 893: if ( l > wlen ) {
! 894: if ( w ) GC_free(w);
! 895: w = (MP *)MALLOC(l*sizeof(MP));
! 896: wlen = l;
! 897: }
! 898: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 899: w[i] = m;
! 900: n = NV(p); n2 = n>>1;
! 901: d0 = m0->dl;
! 902:
! 903: for ( i = 0, tlen = 1; i < n2; i++ )
! 904: tlen *= d0->d[n2+i]+1;
! 905: if ( tlen > rtlen ) {
! 906: if ( tab ) GC_free(tab);
! 907: if ( psum ) GC_free(psum);
! 908: rtlen = tlen;
! 909: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
! 910: psum = (MP *)MALLOC(rtlen*sizeof(MP));
! 911: }
! 912: bzero(psum,tlen*sizeof(MP));
! 913: for ( i = l-1; i >= 0; i-- ) {
! 914: bzero(tab,tlen*sizeof(struct cdl));
! 915: _weyl_mulmm_dup(vl,m0,w[i],n,tab,tlen);
! 916: for ( j = 0; j < tlen; j++ ) {
! 917: if ( tab[j].c ) {
! 918: _NEWMP(m); m->dl = tab[j].d;
! 919: C(m) = tab[j].c; NEXT(m) = psum[j];
! 920: psum[j] = m;
! 921: }
! 922: }
! 923: }
! 924: for ( j = tlen-1, r = 0; j >= 0; j-- )
! 925: if ( psum[j] ) {
! 926: _MKDP(n,psum[j],t); _addd_destructive(vl,r,t,&t1); r = t1;
! 927: }
! 928: if ( r )
! 929: r->sugar = p->sugar + m0->dl->td;
! 930: *pr = r;
! 931: }
! 932: }
! 933:
! 934: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
! 935:
! 936: void _weyl_mulmm_dup(vl,m0,m1,n,rtab,rtablen)
! 937: VL vl;
! 938: MP m0,m1;
! 939: int n;
! 940: struct cdl *rtab;
! 941: int rtablen;
! 942: {
! 943: MP m,mr,mr0;
! 944: DP r,t,t1;
! 945: P c;
! 946: int c0,c1,cc;
! 947: DL d,d0,d1,dt;
! 948: int i,j,a,b,k,l,n2,s,min,h,curlen;
! 949: struct cdl *p;
! 950: static P *ctab;
! 951: static struct cdl *tab;
! 952: static int tablen;
! 953: static struct cdl *tmptab;
! 954: static int tmptablen;
! 955:
! 956: if ( !m0 || !m1 ) {
! 957: rtab[0].c = 0;
! 958: rtab[0].d = 0;
! 959: return;
! 960: }
! 961: mulp(vl,C(m0),C(m1),&c);
! 962: d0 = m0->dl; d1 = m1->dl;
! 963: n2 = n>>1;
! 964: curlen = 1;
! 965:
! 966: _NEWDL(d,n);
! 967: if ( n & 1 )
! 968: /* offset of h-degree */
! 969: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
! 970: else
! 971: d->td = 0;
! 972: rtab[0].c = c;
! 973: rtab[0].d = d;
! 974:
! 975: if ( rtablen > tmptablen ) {
! 976: if ( tmptab ) GC_free(tmptab);
! 977: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
! 978: tmptablen = rtablen;
! 979: }
! 980:
! 981: for ( i = 0; i < n2; i++ ) {
! 982: a = d0->d[i]; b = d1->d[n2+i];
! 983: k = d0->d[n2+i]; l = d1->d[i];
! 984: if ( !k || !l ) {
! 985: a += l;
! 986: b += k;
! 987: s = a+b;
! 988: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
! 989: if ( p->c ) {
! 990: dt = p->d;
! 991: dt->d[i] = a;
! 992: dt->d[n2+i] = b;
! 993: dt->td += s;
! 994: }
! 995: }
! 996: curlen *= k+1;
! 997: continue;
! 998: }
! 999: if ( k+1 > tablen ) {
! 1000: if ( tab ) GC_free(tab);
! 1001: if ( ctab ) GC_free(ctab);
! 1002: tablen = k+1;
! 1003: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
! 1004: ctab = (P *)MALLOC(tablen*sizeof(P));
! 1005: }
! 1006: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 1007: s = a+k+l+b;
! 1008: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 1009: min = MIN(k,l);
! 1010: mkwc(k,l,ctab);
! 1011: bzero(tab,(k+1)*sizeof(struct cdl));
! 1012: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
! 1013: if ( n & 1 )
! 1014: for ( j = 0; j <= min; j++ ) {
! 1015: _NEWDL(d,n);
! 1016: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 1017: d->td = s;
! 1018: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
! 1019: tab[j].d = d;
! 1020: tab[j].c = ctab[j];
! 1021: }
! 1022: else
! 1023: for ( j = 0; j <= min; j++ ) {
! 1024: _NEWDL(d,n);
! 1025: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 1026: d->td = d->d[i]+d->d[n2+i]; /* XXX */
! 1027: tab[j].d = d;
! 1028: tab[j].c = ctab[j];
! 1029: }
! 1030: #if 0
! 1031: _comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
! 1032: for ( j = 0; j < curlen; j++ )
! 1033: if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
! 1034: for ( j = 0; j <= min; j++ )
! 1035: if ( tab[j].d ) { _FREEDL(tab[j].d); }
! 1036: curlen *= k+1;
! 1037: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
! 1038: #else
! 1039: _comm_muld_tab_destructive(vl,n,rtab,curlen,tab,k+1);
! 1040: for ( j = 0; j <= min; j++ )
! 1041: if ( tab[j].d ) { _FREEDL(tab[j].d); }
! 1042: curlen *= k+1;
! 1043: #endif
! 1044: }
! 1045: }
! 1046:
! 1047: /* direct product of two cdl tables
! 1048: rt[] = [
! 1049: t[0]*t1[0],...,t[n-1]*t1[0],
! 1050: t[0]*t1[1],...,t[n-1]*t1[1],
! 1051: ...
! 1052: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
! 1053: ]
! 1054: */
! 1055:
! 1056: void _comm_muld_tab(vl,nv,t,n,t1,n1,rt)
! 1057: VL vl;
! 1058: int nv;
! 1059: struct cdl *t;
! 1060: int n;
! 1061: struct cdl *t1;
! 1062: int n1;
! 1063: struct cdl *rt;
! 1064: {
! 1065: int i,j;
! 1066: struct cdl *p;
! 1067: P c;
! 1068: DL d;
! 1069:
! 1070: bzero(rt,n*n1*sizeof(struct cdl));
! 1071: for ( j = 0, p = rt; j < n1; j++ ) {
! 1072: c = t1[j].c;
! 1073: d = t1[j].d;
! 1074: if ( !c )
! 1075: break;
! 1076: for ( i = 0; i < n; i++, p++ ) {
! 1077: if ( t[i].c ) {
! 1078: mulp(vl,t[i].c,c,&p->c);
! 1079: _adddl_dup(nv,t[i].d,d,&p->d);
! 1080: }
! 1081: }
! 1082: }
! 1083: }
! 1084:
! 1085: void _comm_muld_tab_destructive(vl,nv,t,n,t1,n1)
! 1086: VL vl;
! 1087: int nv;
! 1088: struct cdl *t;
! 1089: int n;
! 1090: struct cdl *t1;
! 1091: int n1;
! 1092: {
! 1093: int i,j;
! 1094: struct cdl *p;
! 1095: P c;
! 1096: DL d;
! 1097:
! 1098: for ( j = 1, p = t+n; j < n1; j++ ) {
! 1099: c = t1[j].c;
! 1100: d = t1[j].d;
! 1101: if ( !c )
! 1102: break;
! 1103: for ( i = 0; i < n; i++, p++ ) {
! 1104: if ( t[i].c ) {
! 1105: mulp(vl,t[i].c,c,&p->c);
! 1106: _adddl_dup(nv,t[i].d,d,&p->d);
! 1107: }
! 1108: }
! 1109: }
! 1110: c = t1[0].c;
! 1111: d = t1[0].d;
! 1112: for ( i = 0, p = t; i < n; i++, p++ )
! 1113: if ( t[i].c ) {
! 1114: mulp(vl,t[i].c,c,&p->c);
! 1115: /* t[i].d += d */
! 1116: adddl_destructive(nv,t[i].d,d);
! 1117: }
1.1 noro 1118: }