![[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.10 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.9 2001/09/11 03:13:43 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;
1.10 ! noro 266: int c,n,r,i,c1,c2;
1.1 noro 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);
1.10 ! noro 275: c1 = ITOS(C(m));
! 276: DMAR(c1,c,0,mod,c2);
! 277: C(mr) = (P)STOI(c2);
1.1 noro 278: _NEWDL_NOINIT(dt,n); mr->dl = dt;
279: dm = m->dl;
280: dt->td = d->td + dm->td;
281: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
282: *pt++ = *p1++ + *p2++;
283: }
284: NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
285: if ( *pr )
286: (*pr)->sugar = p->sugar + m0->dl->td;
287: }
288: }
289:
1.4 noro 290: void _weyl_mulmmm_dup_bug();
291:
292: void _weyl_mulmdm_dup(mod,m0,p,pr)
1.1 noro 293: int mod;
1.4 noro 294: MP m0;
1.1 noro 295: DP p;
296: DP *pr;
297: {
298: DP r,t,t1;
299: MP m;
1.4 noro 300: DL d0;
301: int n,n2,l,i,j,tlen;
302: static MP *w,*psum;
303: static struct cdlm *tab;
1.1 noro 304: static int wlen;
1.4 noro 305: static int rtlen;
1.1 noro 306:
307: if ( !p )
308: *pr = 0;
309: else {
310: for ( m = BDY(p), 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(p), i = 0; i < l; m = NEXT(m), i++ )
317: w[i] = m;
1.4 noro 318: n = NV(p); n2 = n>>1;
319: d0 = m0->dl;
320:
321: for ( i = 0, tlen = 1; i < n2; i++ )
322: tlen *= d0->d[n2+i]+1;
323: if ( tlen > rtlen ) {
324: if ( tab ) GC_free(tab);
325: if ( psum ) GC_free(psum);
326: rtlen = tlen;
327: tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
328: psum = (MP *)MALLOC(rtlen*sizeof(MP));
1.1 noro 329: }
1.4 noro 330: bzero(psum,tlen*sizeof(MP));
331: for ( i = l-1; i >= 0; i-- ) {
332: bzero(tab,tlen*sizeof(struct cdlm));
333: _weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
334: for ( j = 0; j < tlen; j++ ) {
335: if ( tab[j].c ) {
336: _NEWMP(m); m->dl = tab[j].d;
337: C(m) = STOI(tab[j].c); NEXT(m) = psum[j];
338: psum[j] = m;
339: }
340: }
341: }
342: for ( j = tlen-1, r = 0; j >= 0; j-- )
343: if ( psum[j] ) {
344: _MKDP(n,psum[j],t); _addmd_destructive(mod,r,t,&t1); r = t1;
345: }
1.1 noro 346: if ( r )
347: r->sugar = p->sugar + m0->dl->td;
348: *pr = r;
349: }
350: }
1.4 noro 351:
1.1 noro 352: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
353:
1.4 noro 354: void _weyl_mulmmm_dup(mod,m0,m1,n,rtab,rtablen)
1.1 noro 355: int mod;
356: MP m0,m1;
357: int n;
1.4 noro 358: struct cdlm *rtab;
359: int rtablen;
1.1 noro 360: {
361: MP m,mr,mr0;
362: DP r,t,t1;
363: int c,c0,c1,cc;
1.4 noro 364: DL d,d0,d1,dt;
365: int i,j,a,b,k,l,n2,s,min,h,curlen;
366: struct cdlm *p;
367: static int *ctab;
368: static struct cdlm *tab;
1.1 noro 369: static int tablen;
1.4 noro 370: static struct cdlm *tmptab;
371: static int tmptablen;
1.1 noro 372:
1.4 noro 373: if ( !m0 || !m1 ) {
374: rtab[0].c = 0;
375: rtab[0].d = 0;
376: return;
377: }
378: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
379: c = dmar(c0,c1,0,mod);
380: d0 = m0->dl; d1 = m1->dl;
381: n2 = n>>1;
382: curlen = 1;
383:
384: _NEWDL(d,n);
385: if ( n & 1 )
386: /* offset of h-degree */
387: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
388: else
389: d->td = 0;
390: rtab[0].c = c;
391: rtab[0].d = d;
392:
393: if ( rtablen > tmptablen ) {
394: if ( tmptab ) GC_free(tmptab);
395: tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
396: tmptablen = rtablen;
397: }
1.1 noro 398:
1.4 noro 399: for ( i = 0; i < n2; i++ ) {
400: a = d0->d[i]; b = d1->d[n2+i];
401: k = d0->d[n2+i]; l = d1->d[i];
402: if ( !k || !l ) {
403: a += l;
404: b += k;
405: s = a+b;
406: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
407: if ( p->c ) {
408: dt = p->d;
409: dt->d[i] = a;
410: dt->d[n2+i] = b;
411: dt->td += s;
412: }
413: }
414: curlen *= k+1;
415: continue;
416: }
417: if ( k+1 > tablen ) {
418: if ( tab ) GC_free(tab);
419: if ( ctab ) GC_free(ctab);
420: tablen = k+1;
421: tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
422: ctab = (int *)MALLOC(tablen*sizeof(int));
423: }
424: /* degree of xi^a*(Di^k*xi^l)*Di^b */
425: s = a+k+l+b;
426: /* compute xi^a*(Di^k*xi^l)*Di^b */
427: min = MIN(k,l);
428: mkwcm(k,l,mod,ctab);
429: bzero(tab,(k+1)*sizeof(struct cdlm));
430: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
1.1 noro 431: if ( n & 1 )
1.4 noro 432: for ( j = 0; j <= min; j++ ) {
433: _NEWDL(d,n);
434: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
435: d->td = s;
436: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
437: tab[j].d = d;
438: tab[j].c = ctab[j];
439: }
1.1 noro 440: else
1.4 noro 441: for ( j = 0; j <= min; j++ ) {
442: _NEWDL(d,n);
443: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
444: d->td = d->d[i]+d->d[n2+i]; /* XXX */
445: tab[j].d = d;
446: tab[j].c = ctab[j];
447: }
448: #if 0
449: _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
450: for ( j = 0; j < curlen; j++ )
1.5 noro 451: if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
1.4 noro 452: for ( j = 0; j <= min; j++ )
1.5 noro 453: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4 noro 454: curlen *= k+1;
455: bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
456: #else
457: _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
458: for ( j = 0; j <= min; j++ )
1.5 noro 459: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4 noro 460: curlen *= k+1;
461: #endif
462: }
463: }
464:
465: /* direct product of two cdlm tables
466: rt[] = [
467: t[0]*t1[0],...,t[n-1]*t1[0],
468: t[0]*t1[1],...,t[n-1]*t1[1],
469: ...
470: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
471: ]
472: */
473:
474: void _comm_mulmd_tab(mod,nv,t,n,t1,n1,rt)
475: int mod;
476: int nv;
477: struct cdlm *t;
478: int n;
479: struct cdlm *t1;
480: int n1;
481: struct cdlm *rt;
482: {
483: int i,j;
484: struct cdlm *p;
485: int c;
486: DL d;
487:
488: bzero(rt,n*n1*sizeof(struct cdlm));
489: for ( j = 0, p = rt; j < n1; j++ ) {
490: c = t1[j].c;
491: d = t1[j].d;
492: if ( !c )
493: break;
494: for ( i = 0; i < n; i++, p++ ) {
495: if ( t[i].c ) {
496: p->c = dmar(t[i].c,c,0,mod);
1.5 noro 497: _adddl_dup(nv,t[i].d,d,&p->d);
1.4 noro 498: }
499: }
500: }
501: }
502:
503: void _comm_mulmd_tab_destructive(mod,nv,t,n,t1,n1)
504: int mod;
505: int nv;
506: struct cdlm *t;
507: int n;
508: struct cdlm *t1;
509: int n1;
510: {
511: int i,j;
512: struct cdlm *p;
513: int c;
514: DL d;
515:
516: for ( j = 1, p = t+n; j < n1; j++ ) {
517: c = t1[j].c;
518: d = t1[j].d;
519: if ( !c )
520: break;
521: for ( i = 0; i < n; i++, p++ ) {
522: if ( t[i].c ) {
523: p->c = dmar(t[i].c,c,0,mod);
1.5 noro 524: _adddl_dup(nv,t[i].d,d,&p->d);
1.1 noro 525: }
526: }
527: }
1.4 noro 528: c = t1[0].c;
529: d = t1[0].d;
530: for ( i = 0, p = t; i < n; i++, p++ )
531: if ( t[i].c ) {
532: p->c = dmar(t[i].c,c,0,mod);
533: /* t[i].d += d */
534: adddl_destructive(nv,t[i].d,d);
535: }
1.1 noro 536: }
537:
1.5 noro 538: void dlto_dl(d,dr)
539: DL d;
540: DL *dr;
1.1 noro 541: {
542: int i,n;
1.5 noro 543: DL t;
1.1 noro 544:
1.5 noro 545: n = current_dl_length;
546: _NEWDL(t,n); *dr = t;
547: t->td = d->td;
1.1 noro 548: for ( i = 0; i < n; i++ )
1.5 noro 549: t->d[i] = d->d[i];
1.1 noro 550: }
551:
552: void _dltodl(d,dr)
553: DL d;
554: DL *dr;
555: {
556: int i,n;
557: DL t;
558:
559: n = current_dl_length;
560: NEWDL(t,n); *dr = t;
561: t->td = d->td;
562: for ( i = 0; i < n; i++ )
563: t->d[i] = d->d[i];
564: }
565:
1.5 noro 566: void _adddl_dup(n,d1,d2,dr)
1.1 noro 567: int n;
568: DL d1,d2;
569: DL *dr;
570: {
571: DL dt;
572: int i;
573:
1.4 noro 574: _NEWDL(dt,n);
575: *dr = dt;
1.1 noro 576: dt->td = d1->td + d2->td;
577: for ( i = 0; i < n; i++ )
578: dt->d[i] = d1->d[i]+d2->d[i];
1.4 noro 579: }
580:
581: void _free_dlarray(a,n)
582: DL *a;
583: int n;
584: {
585: int i;
586:
1.5 noro 587: for ( i = 0; i < n; i++ ) { _FREEDL(a[i]); }
1.1 noro 588: }
589:
590: void _free_dp(f)
591: DP f;
592: {
593: MP m,s;
594:
595: if ( !f )
596: return;
597: m = f->body;
598: while ( m ) {
1.5 noro 599: s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
600: }
601: _FREEDP(f);
602: }
603:
604: void dpto_dp(p,r)
605: DP p;
606: DP *r;
607: {
608: MP m,mr0,mr;
1.7 noro 609: DL t;
1.5 noro 610:
611: if ( !p )
612: *r = 0;
613: else {
1.7 noro 614: /* XXX : dummy call to set current_dl_length */
615: _NEWDL_NOINIT(t,NV(p));
616:
1.5 noro 617: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
618: _NEXTMP(mr0,mr);
619: dlto_dl(m->dl,&mr->dl);
620: mr->c = m->c;
621: }
622: NEXT(mr) = 0;
623: _MKDP(p->nv,mr0,*r);
624: (*r)->sugar = p->sugar;
1.1 noro 625: }
626: }
627:
628: void _dptodp(p,r)
629: DP p;
630: DP *r;
631: {
632: MP m,mr0,mr;
633:
634: if ( !p )
635: *r = 0;
636: else {
637: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
638: NEXTMP(mr0,mr);
639: _dltodl(m->dl,&mr->dl);
640: mr->c = m->c;
641: }
642: NEXT(mr) = 0;
643: MKDP(p->nv,mr0,*r);
644: (*r)->sugar = p->sugar;
645: }
646: }
647:
1.5 noro 648: /*
649: * destructive merge of two list
650: *
651: * p1, p2 : list of DL
652: * return : a merged list
653: */
654:
655: NODE _symb_merge(m1,m2,n)
656: NODE m1,m2;
657: int n;
658: {
659: NODE top,prev,cur,m,t;
660:
661: if ( !m1 )
662: return m2;
663: else if ( !m2 )
664: return m1;
665: else {
666: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
667: case 0:
668: top = m1; _FREEDL((DL)BDY(m2)); m = NEXT(m2);
669: break;
670: case 1:
671: top = m1; m = m2;
672: break;
673: case -1:
674: top = m2; m = m1;
1.1 noro 675: break;
676: }
1.5 noro 677: prev = top; cur = NEXT(top);
678: /* BDY(prev) > BDY(m) always holds */
679: while ( cur && m ) {
680: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
681: case 0:
682: _FREEDL(BDY(m)); m = NEXT(m);
683: prev = cur; cur = NEXT(cur);
684: break;
685: case 1:
686: t = NEXT(cur); NEXT(cur) = m; m = t;
687: prev = cur; cur = NEXT(cur);
688: break;
689: case -1:
690: NEXT(prev) = m; m = cur;
691: prev = NEXT(prev); cur = NEXT(prev);
692: break;
1.1 noro 693: }
694: }
1.5 noro 695: if ( !cur )
696: NEXT(prev) = m;
697: return top;
1.1 noro 698: }
1.9 noro 699: }
700:
701: /* XXX : bellow should be placed in another file */
702:
703: void dpto_dp();
704: void _dptodp();
705: void _adddl_dup();
706: void adddl_destructive();
707: void _muldm_dup();
708: void _free_dp();
709: void _comm_muld_dup();
710: void _weyl_muld_dup();
711: void _weyl_mulmm_dup();
712: void _weyl_muldm_dup();
713: void _comm_muld_tab();
714: void _comm_muld_tab_destructive();
715:
716: /* merge p1 and p2 into pr */
717:
718: void _addd_destructive(vl,p1,p2,pr)
719: VL vl;
720: DP p1,p2,*pr;
721: {
722: int n;
723: MP m1,m2,mr,mr0,s;
724: P t;
725:
726: if ( !p1 )
727: *pr = p2;
728: else if ( !p2 )
729: *pr = p1;
730: else {
731: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
732: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
733: case 0:
734: addp(vl,C(m1),C(m2),&t);
735: s = m1; m1 = NEXT(m1);
736: if ( t ) {
737: _NEXTMP2(mr0,mr,s); C(mr) = t;
738: } else {
739: _FREEDL(s->dl); _FREEMP(s);
740: }
741: s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
742: break;
743: case 1:
744: s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
745: break;
746: case -1:
747: s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
748: break;
749: }
750: if ( !mr0 )
751: if ( m1 )
752: mr0 = m1;
753: else if ( m2 )
754: mr0 = m2;
755: else {
756: *pr = 0;
757: return;
758: }
759: else if ( m1 )
760: NEXT(mr) = m1;
761: else if ( m2 )
762: NEXT(mr) = m2;
763: else
764: NEXT(mr) = 0;
765: _MKDP(NV(p1),mr0,*pr);
766: if ( *pr )
767: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
768: _FREEDP(p1); _FREEDP(p2);
769: }
770: }
771:
772: void _muld_dup(vl,p1,p2,pr)
773: VL vl;
774: DP p1,p2,*pr;
775: {
776: if ( !do_weyl )
777: _comm_muld_dup(vl,p1,p2,pr);
778: else
779: _weyl_muld_dup(vl,p1,p2,pr);
780: }
781:
782: void _comm_muld_dup(vl,p1,p2,pr)
783: VL vl;
784: DP p1,p2,*pr;
785: {
786: MP m;
787: DP s,t,u;
788: int i,l,l1;
789: static MP *w;
790: static int wlen;
791:
792: if ( !p1 || !p2 )
793: *pr = 0;
794: else {
795: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
796: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
797: if ( l1 < l ) {
798: t = p1; p1 = p2; p2 = t;
799: l = l1;
800: }
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(p2), i = 0; i < l; m = NEXT(m), i++ )
807: w[i] = m;
808: for ( s = 0, i = l-1; i >= 0; i-- ) {
809: _muldm_dup(vl,p1,w[i],&t); _addd_destructive(vl,s,t,&u); s = u;
810: }
811: bzero(w,l*sizeof(MP));
812: *pr = s;
813: }
814: }
815:
816: void _weyl_muld_dup(vl,p1,p2,pr)
817: VL vl;
818: DP p1,p2,*pr;
819: {
820: MP m;
821: DP s,t,u;
822: int i,l,l1;
823: static MP *w;
824: static int wlen;
825:
826: if ( !p1 || !p2 )
827: *pr = 0;
828: else {
829: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
830: if ( l > wlen ) {
831: if ( w ) GC_free(w);
832: w = (MP *)MALLOC(l*sizeof(MP));
833: wlen = l;
834: }
835: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
836: w[i] = m;
837: for ( s = 0, i = l-1; i >= 0; i-- ) {
838: _weyl_muldm_dup(vl,w[i],p2,&t); _addd_destructive(vl,s,t,&u); s = u;
839: }
840: bzero(w,l*sizeof(MP));
841: *pr = s;
842: }
843: }
844:
845: void _muldm_dup(vl,p,m0,pr)
846: VL vl;
847: DP p;
848: MP m0;
849: DP *pr;
850: {
851: MP m,mr,mr0;
852: DL d,dt,dm;
853: P c;
854: int n,r,i;
855: int *pt,*p1,*p2;
856:
857: if ( !p )
858: *pr = 0;
859: else {
860: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
861: m; m = NEXT(m) ) {
862: _NEXTMP(mr0,mr);
863: mulp(vl,C(m),c,&C(mr));
864: _NEWDL_NOINIT(dt,n); mr->dl = dt;
865: dm = m->dl;
866: dt->td = d->td + dm->td;
867: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
868: *pt++ = *p1++ + *p2++;
869: }
870: NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
871: if ( *pr )
872: (*pr)->sugar = p->sugar + m0->dl->td;
873: }
874: }
875:
876: void _weyl_muldm_dup(vl,m0,p,pr)
877: VL vl;
878: MP m0;
879: DP p;
880: DP *pr;
881: {
882: DP r,t,t1;
883: MP m;
884: DL d0;
885: int n,n2,l,i,j,tlen;
886: static MP *w,*psum;
887: static struct cdl *tab;
888: static int wlen;
889: static int rtlen;
890:
891: if ( !p )
892: *pr = 0;
893: else {
894: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
895: if ( l > wlen ) {
896: if ( w ) GC_free(w);
897: w = (MP *)MALLOC(l*sizeof(MP));
898: wlen = l;
899: }
900: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
901: w[i] = m;
902: n = NV(p); n2 = n>>1;
903: d0 = m0->dl;
904:
905: for ( i = 0, tlen = 1; i < n2; i++ )
906: tlen *= d0->d[n2+i]+1;
907: if ( tlen > rtlen ) {
908: if ( tab ) GC_free(tab);
909: if ( psum ) GC_free(psum);
910: rtlen = tlen;
911: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
912: psum = (MP *)MALLOC(rtlen*sizeof(MP));
913: }
914: bzero(psum,tlen*sizeof(MP));
915: for ( i = l-1; i >= 0; i-- ) {
916: bzero(tab,tlen*sizeof(struct cdl));
917: _weyl_mulmm_dup(vl,m0,w[i],n,tab,tlen);
918: for ( j = 0; j < tlen; j++ ) {
919: if ( tab[j].c ) {
920: _NEWMP(m); m->dl = tab[j].d;
921: C(m) = tab[j].c; NEXT(m) = psum[j];
922: psum[j] = m;
923: }
924: }
925: }
926: for ( j = tlen-1, r = 0; j >= 0; j-- )
927: if ( psum[j] ) {
928: _MKDP(n,psum[j],t); _addd_destructive(vl,r,t,&t1); r = t1;
929: }
930: if ( r )
931: r->sugar = p->sugar + m0->dl->td;
932: *pr = r;
933: }
934: }
935:
936: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
937:
938: void _weyl_mulmm_dup(vl,m0,m1,n,rtab,rtablen)
939: VL vl;
940: MP m0,m1;
941: int n;
942: struct cdl *rtab;
943: int rtablen;
944: {
945: MP m,mr,mr0;
946: DP r,t,t1;
947: P c;
948: int c0,c1,cc;
949: DL d,d0,d1,dt;
950: int i,j,a,b,k,l,n2,s,min,h,curlen;
951: struct cdl *p;
952: static P *ctab;
953: static struct cdl *tab;
954: static int tablen;
955: static struct cdl *tmptab;
956: static int tmptablen;
957:
958: if ( !m0 || !m1 ) {
959: rtab[0].c = 0;
960: rtab[0].d = 0;
961: return;
962: }
963: mulp(vl,C(m0),C(m1),&c);
964: d0 = m0->dl; d1 = m1->dl;
965: n2 = n>>1;
966: curlen = 1;
967:
968: _NEWDL(d,n);
969: if ( n & 1 )
970: /* offset of h-degree */
971: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
972: else
973: d->td = 0;
974: rtab[0].c = c;
975: rtab[0].d = d;
976:
977: if ( rtablen > tmptablen ) {
978: if ( tmptab ) GC_free(tmptab);
979: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
980: tmptablen = rtablen;
981: }
982:
983: for ( i = 0; i < n2; i++ ) {
984: a = d0->d[i]; b = d1->d[n2+i];
985: k = d0->d[n2+i]; l = d1->d[i];
986: if ( !k || !l ) {
987: a += l;
988: b += k;
989: s = a+b;
990: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
991: if ( p->c ) {
992: dt = p->d;
993: dt->d[i] = a;
994: dt->d[n2+i] = b;
995: dt->td += s;
996: }
997: }
998: curlen *= k+1;
999: continue;
1000: }
1001: if ( k+1 > tablen ) {
1002: if ( tab ) GC_free(tab);
1003: if ( ctab ) GC_free(ctab);
1004: tablen = k+1;
1005: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
1006: ctab = (P *)MALLOC(tablen*sizeof(P));
1007: }
1008: /* degree of xi^a*(Di^k*xi^l)*Di^b */
1009: s = a+k+l+b;
1010: /* compute xi^a*(Di^k*xi^l)*Di^b */
1011: min = MIN(k,l);
1012: mkwc(k,l,ctab);
1013: bzero(tab,(k+1)*sizeof(struct cdl));
1014: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
1015: if ( n & 1 )
1016: for ( j = 0; j <= min; j++ ) {
1017: _NEWDL(d,n);
1018: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1019: d->td = s;
1020: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
1021: tab[j].d = d;
1022: tab[j].c = ctab[j];
1023: }
1024: else
1025: for ( j = 0; j <= min; j++ ) {
1026: _NEWDL(d,n);
1027: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1028: d->td = d->d[i]+d->d[n2+i]; /* XXX */
1029: tab[j].d = d;
1030: tab[j].c = ctab[j];
1031: }
1032: #if 0
1033: _comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
1034: for ( j = 0; j < curlen; j++ )
1035: if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
1036: for ( j = 0; j <= min; j++ )
1037: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1038: curlen *= k+1;
1039: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
1040: #else
1041: _comm_muld_tab_destructive(vl,n,rtab,curlen,tab,k+1);
1042: for ( j = 0; j <= min; j++ )
1043: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1044: curlen *= k+1;
1045: #endif
1046: }
1047: }
1048:
1049: /* direct product of two cdl tables
1050: rt[] = [
1051: t[0]*t1[0],...,t[n-1]*t1[0],
1052: t[0]*t1[1],...,t[n-1]*t1[1],
1053: ...
1054: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
1055: ]
1056: */
1057:
1058: void _comm_muld_tab(vl,nv,t,n,t1,n1,rt)
1059: VL vl;
1060: int nv;
1061: struct cdl *t;
1062: int n;
1063: struct cdl *t1;
1064: int n1;
1065: struct cdl *rt;
1066: {
1067: int i,j;
1068: struct cdl *p;
1069: P c;
1070: DL d;
1071:
1072: bzero(rt,n*n1*sizeof(struct cdl));
1073: for ( j = 0, p = rt; j < n1; j++ ) {
1074: c = t1[j].c;
1075: d = t1[j].d;
1076: if ( !c )
1077: break;
1078: for ( i = 0; i < n; i++, p++ ) {
1079: if ( t[i].c ) {
1080: mulp(vl,t[i].c,c,&p->c);
1081: _adddl_dup(nv,t[i].d,d,&p->d);
1082: }
1083: }
1084: }
1085: }
1086:
1087: void _comm_muld_tab_destructive(vl,nv,t,n,t1,n1)
1088: VL vl;
1089: int nv;
1090: struct cdl *t;
1091: int n;
1092: struct cdl *t1;
1093: int n1;
1094: {
1095: int i,j;
1096: struct cdl *p;
1097: P c;
1098: DL d;
1099:
1100: for ( j = 1, p = t+n; j < n1; j++ ) {
1101: c = t1[j].c;
1102: d = t1[j].d;
1103: if ( !c )
1104: break;
1105: for ( i = 0; i < n; i++, p++ ) {
1106: if ( t[i].c ) {
1107: mulp(vl,t[i].c,c,&p->c);
1108: _adddl_dup(nv,t[i].d,d,&p->d);
1109: }
1110: }
1111: }
1112: c = t1[0].c;
1113: d = t1[0].d;
1114: for ( i = 0, p = t; i < n; i++, p++ )
1115: if ( t[i].c ) {
1116: mulp(vl,t[i].c,c,&p->c);
1117: /* t[i].d += d */
1118: adddl_destructive(nv,t[i].d,d);
1119: }
1.1 noro 1120: }