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