![[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.13 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.12 2002/01/28 00:54:42 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: for ( i = 0; i < 128; i++, p += dl_len ) {
1.13 ! noro 84: p = (int *)GC_malloc(dl_len*sizeof(int));
1.1 noro 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: for ( i = 0; i < 1024; i++ ) {
1.13 ! noro 98: p = (MP)GC_malloc(sizeof(struct oMP));
! 99: p->next = _mp_free_list; _mp_free_list = p;
1.1 noro 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: for ( i = 0; i < 1024; i++ ) {
1.13 ! noro 111: p = (DP)GC_malloc(sizeof(struct oDP));
! 112: p->body = (MP)_dp_free_list; _dp_free_list = p;
1.1 noro 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];
1.12 noro 366:
367: /* degree of xi^a*(Di^k*xi^l)*Di^b */
368: a += l;
369: b += k;
370: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
371:
1.4 noro 372: if ( !k || !l ) {
373: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
374: if ( p->c ) {
375: dt = p->d;
376: dt->d[i] = a;
377: dt->d[n2+i] = b;
378: dt->td += s;
379: }
380: }
381: curlen *= k+1;
382: continue;
383: }
384: if ( k+1 > tablen ) {
385: if ( tab ) GC_free(tab);
386: if ( ctab ) GC_free(ctab);
387: tablen = k+1;
388: tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
389: ctab = (int *)MALLOC(tablen*sizeof(int));
390: }
391: /* compute xi^a*(Di^k*xi^l)*Di^b */
392: min = MIN(k,l);
393: mkwcm(k,l,mod,ctab);
394: bzero(tab,(k+1)*sizeof(struct cdlm));
395: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
1.1 noro 396: if ( n & 1 )
1.4 noro 397: for ( j = 0; j <= min; j++ ) {
398: _NEWDL(d,n);
1.12 noro 399: d->d[i] = a-j; d->d[n2+i] = b-j;
1.4 noro 400: d->td = s;
1.12 noro 401: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.4 noro 402: tab[j].d = d;
403: tab[j].c = ctab[j];
404: }
1.1 noro 405: else
1.4 noro 406: for ( j = 0; j <= min; j++ ) {
407: _NEWDL(d,n);
1.12 noro 408: d->d[i] = a-j; d->d[n2+i] = b-j;
409: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.4 noro 410: tab[j].d = d;
411: tab[j].c = ctab[j];
412: }
413: #if 0
414: _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
415: for ( j = 0; j < curlen; j++ )
1.5 noro 416: if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
1.4 noro 417: for ( j = 0; j <= min; j++ )
1.5 noro 418: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4 noro 419: curlen *= k+1;
420: bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
421: #else
422: _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
423: for ( j = 0; j <= min; j++ )
1.5 noro 424: if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4 noro 425: curlen *= k+1;
426: #endif
427: }
428: }
429:
430: /* direct product of two cdlm tables
431: rt[] = [
432: t[0]*t1[0],...,t[n-1]*t1[0],
433: t[0]*t1[1],...,t[n-1]*t1[1],
434: ...
435: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
436: ]
437: */
438:
1.11 noro 439: void _comm_mulmd_tab(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1,struct cdlm *rt)
1.4 noro 440: {
441: int i,j;
442: struct cdlm *p;
443: int c;
444: DL d;
445:
446: bzero(rt,n*n1*sizeof(struct cdlm));
447: for ( j = 0, p = rt; j < n1; j++ ) {
448: c = t1[j].c;
449: d = t1[j].d;
450: if ( !c )
451: break;
452: for ( i = 0; i < n; i++, p++ ) {
453: if ( t[i].c ) {
454: p->c = dmar(t[i].c,c,0,mod);
1.5 noro 455: _adddl_dup(nv,t[i].d,d,&p->d);
1.4 noro 456: }
457: }
458: }
459: }
460:
1.11 noro 461: void _comm_mulmd_tab_destructive(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1)
1.4 noro 462: {
463: int i,j;
464: struct cdlm *p;
465: int c;
466: DL d;
467:
468: for ( j = 1, p = t+n; j < n1; j++ ) {
469: c = t1[j].c;
470: d = t1[j].d;
471: if ( !c )
472: break;
473: for ( i = 0; i < n; i++, p++ ) {
474: if ( t[i].c ) {
475: p->c = dmar(t[i].c,c,0,mod);
1.5 noro 476: _adddl_dup(nv,t[i].d,d,&p->d);
1.1 noro 477: }
478: }
479: }
1.4 noro 480: c = t1[0].c;
481: d = t1[0].d;
482: for ( i = 0, p = t; i < n; i++, p++ )
483: if ( t[i].c ) {
484: p->c = dmar(t[i].c,c,0,mod);
485: /* t[i].d += d */
486: adddl_destructive(nv,t[i].d,d);
487: }
1.1 noro 488: }
489:
1.11 noro 490: void dlto_dl(DL d,DL *dr)
1.1 noro 491: {
492: int i,n;
1.5 noro 493: DL t;
1.1 noro 494:
1.5 noro 495: n = current_dl_length;
496: _NEWDL(t,n); *dr = t;
497: t->td = d->td;
1.1 noro 498: for ( i = 0; i < n; i++ )
1.5 noro 499: t->d[i] = d->d[i];
1.1 noro 500: }
501:
1.11 noro 502: void _dltodl(DL d,DL *dr)
1.1 noro 503: {
504: int i,n;
505: DL t;
506:
507: n = current_dl_length;
508: NEWDL(t,n); *dr = t;
509: t->td = d->td;
510: for ( i = 0; i < n; i++ )
511: t->d[i] = d->d[i];
512: }
513:
1.11 noro 514: void _adddl_dup(int n,DL d1,DL d2,DL *dr)
1.1 noro 515: {
516: DL dt;
517: int i;
518:
1.4 noro 519: _NEWDL(dt,n);
520: *dr = dt;
1.1 noro 521: dt->td = d1->td + d2->td;
522: for ( i = 0; i < n; i++ )
523: dt->d[i] = d1->d[i]+d2->d[i];
1.4 noro 524: }
525:
1.11 noro 526: void _free_dlarray(DL *a,int n)
1.4 noro 527: {
528: int i;
529:
1.5 noro 530: for ( i = 0; i < n; i++ ) { _FREEDL(a[i]); }
1.1 noro 531: }
532:
1.11 noro 533: void _free_dp(DP f)
1.1 noro 534: {
535: MP m,s;
536:
537: if ( !f )
538: return;
539: m = f->body;
540: while ( m ) {
1.5 noro 541: s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
542: }
543: _FREEDP(f);
544: }
545:
1.11 noro 546: void dpto_dp(DP p,DP *r)
1.5 noro 547: {
548: MP m,mr0,mr;
1.7 noro 549: DL t;
1.5 noro 550:
551: if ( !p )
552: *r = 0;
553: else {
1.7 noro 554: /* XXX : dummy call to set current_dl_length */
555: _NEWDL_NOINIT(t,NV(p));
556:
1.5 noro 557: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
558: _NEXTMP(mr0,mr);
559: dlto_dl(m->dl,&mr->dl);
560: mr->c = m->c;
561: }
562: NEXT(mr) = 0;
563: _MKDP(p->nv,mr0,*r);
564: (*r)->sugar = p->sugar;
1.1 noro 565: }
566: }
567:
1.11 noro 568: void _dptodp(DP p,DP *r)
1.1 noro 569: {
570: MP m,mr0,mr;
571:
572: if ( !p )
573: *r = 0;
574: else {
575: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
576: NEXTMP(mr0,mr);
577: _dltodl(m->dl,&mr->dl);
578: mr->c = m->c;
579: }
580: NEXT(mr) = 0;
581: MKDP(p->nv,mr0,*r);
582: (*r)->sugar = p->sugar;
583: }
584: }
585:
1.5 noro 586: /*
587: * destructive merge of two list
588: *
589: * p1, p2 : list of DL
590: * return : a merged list
591: */
592:
1.11 noro 593: NODE _symb_merge(NODE m1,NODE m2,int n)
1.5 noro 594: {
595: NODE top,prev,cur,m,t;
596:
597: if ( !m1 )
598: return m2;
599: else if ( !m2 )
600: return m1;
601: else {
602: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
603: case 0:
604: top = m1; _FREEDL((DL)BDY(m2)); m = NEXT(m2);
605: break;
606: case 1:
607: top = m1; m = m2;
608: break;
609: case -1:
610: top = m2; m = m1;
1.1 noro 611: break;
612: }
1.5 noro 613: prev = top; cur = NEXT(top);
614: /* BDY(prev) > BDY(m) always holds */
615: while ( cur && m ) {
616: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
617: case 0:
618: _FREEDL(BDY(m)); m = NEXT(m);
619: prev = cur; cur = NEXT(cur);
620: break;
621: case 1:
622: t = NEXT(cur); NEXT(cur) = m; m = t;
623: prev = cur; cur = NEXT(cur);
624: break;
625: case -1:
626: NEXT(prev) = m; m = cur;
627: prev = NEXT(prev); cur = NEXT(prev);
628: break;
1.1 noro 629: }
630: }
1.5 noro 631: if ( !cur )
632: NEXT(prev) = m;
633: return top;
1.1 noro 634: }
1.9 noro 635: }
636:
637: /* merge p1 and p2 into pr */
638:
1.11 noro 639: void _addd_destructive(VL vl,DP p1,DP p2,DP *pr)
1.9 noro 640: {
641: int n;
642: MP m1,m2,mr,mr0,s;
643: P t;
644:
645: if ( !p1 )
646: *pr = p2;
647: else if ( !p2 )
648: *pr = p1;
649: else {
650: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
651: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
652: case 0:
653: addp(vl,C(m1),C(m2),&t);
654: s = m1; m1 = NEXT(m1);
655: if ( t ) {
656: _NEXTMP2(mr0,mr,s); C(mr) = t;
657: } else {
658: _FREEDL(s->dl); _FREEMP(s);
659: }
660: s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
661: break;
662: case 1:
663: s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
664: break;
665: case -1:
666: s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
667: break;
668: }
669: if ( !mr0 )
670: if ( m1 )
671: mr0 = m1;
672: else if ( m2 )
673: mr0 = m2;
674: else {
675: *pr = 0;
676: return;
677: }
678: else if ( m1 )
679: NEXT(mr) = m1;
680: else if ( m2 )
681: NEXT(mr) = m2;
682: else
683: NEXT(mr) = 0;
684: _MKDP(NV(p1),mr0,*pr);
685: if ( *pr )
686: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
687: _FREEDP(p1); _FREEDP(p2);
688: }
689: }
690:
1.11 noro 691: void _muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9 noro 692: {
693: if ( !do_weyl )
694: _comm_muld_dup(vl,p1,p2,pr);
695: else
696: _weyl_muld_dup(vl,p1,p2,pr);
697: }
698:
1.11 noro 699: void _comm_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9 noro 700: {
701: MP m;
702: DP s,t,u;
703: int i,l,l1;
704: static MP *w;
705: static int wlen;
706:
707: if ( !p1 || !p2 )
708: *pr = 0;
709: else {
710: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
711: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
712: if ( l1 < l ) {
713: t = p1; p1 = p2; p2 = t;
714: l = l1;
715: }
716: if ( l > wlen ) {
717: if ( w ) GC_free(w);
718: w = (MP *)MALLOC(l*sizeof(MP));
719: wlen = l;
720: }
721: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
722: w[i] = m;
723: for ( s = 0, i = l-1; i >= 0; i-- ) {
724: _muldm_dup(vl,p1,w[i],&t); _addd_destructive(vl,s,t,&u); s = u;
725: }
726: bzero(w,l*sizeof(MP));
727: *pr = s;
728: }
729: }
730:
1.11 noro 731: void _weyl_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9 noro 732: {
733: MP m;
734: DP s,t,u;
1.11 noro 735: int i,l;
1.9 noro 736: static MP *w;
737: static int wlen;
738:
739: if ( !p1 || !p2 )
740: *pr = 0;
741: else {
742: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
743: if ( l > wlen ) {
744: if ( w ) GC_free(w);
745: w = (MP *)MALLOC(l*sizeof(MP));
746: wlen = l;
747: }
748: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
749: w[i] = m;
750: for ( s = 0, i = l-1; i >= 0; i-- ) {
751: _weyl_muldm_dup(vl,w[i],p2,&t); _addd_destructive(vl,s,t,&u); s = u;
752: }
753: bzero(w,l*sizeof(MP));
754: *pr = s;
755: }
756: }
757:
1.11 noro 758: void _muldm_dup(VL vl,DP p,MP m0,DP *pr)
1.9 noro 759: {
760: MP m,mr,mr0;
761: DL d,dt,dm;
762: P c;
1.11 noro 763: int n,i;
1.9 noro 764: int *pt,*p1,*p2;
765:
766: if ( !p )
767: *pr = 0;
768: else {
769: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
770: m; m = NEXT(m) ) {
771: _NEXTMP(mr0,mr);
772: mulp(vl,C(m),c,&C(mr));
773: _NEWDL_NOINIT(dt,n); mr->dl = dt;
774: dm = m->dl;
775: dt->td = d->td + dm->td;
776: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
777: *pt++ = *p1++ + *p2++;
778: }
779: NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
780: if ( *pr )
781: (*pr)->sugar = p->sugar + m0->dl->td;
782: }
783: }
784:
1.11 noro 785: void _weyl_muldm_dup(VL vl,MP m0,DP p,DP *pr)
1.9 noro 786: {
787: DP r,t,t1;
788: MP m;
789: DL d0;
790: int n,n2,l,i,j,tlen;
791: static MP *w,*psum;
792: static struct cdl *tab;
793: static int wlen;
794: static int rtlen;
795:
796: if ( !p )
797: *pr = 0;
798: else {
799: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
800: if ( l > wlen ) {
801: if ( w ) GC_free(w);
802: w = (MP *)MALLOC(l*sizeof(MP));
803: wlen = l;
804: }
805: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
806: w[i] = m;
807: n = NV(p); n2 = n>>1;
808: d0 = m0->dl;
809:
810: for ( i = 0, tlen = 1; i < n2; i++ )
811: tlen *= d0->d[n2+i]+1;
812: if ( tlen > rtlen ) {
813: if ( tab ) GC_free(tab);
814: if ( psum ) GC_free(psum);
815: rtlen = tlen;
816: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
817: psum = (MP *)MALLOC(rtlen*sizeof(MP));
818: }
819: bzero(psum,tlen*sizeof(MP));
820: for ( i = l-1; i >= 0; i-- ) {
821: bzero(tab,tlen*sizeof(struct cdl));
822: _weyl_mulmm_dup(vl,m0,w[i],n,tab,tlen);
823: for ( j = 0; j < tlen; j++ ) {
824: if ( tab[j].c ) {
825: _NEWMP(m); m->dl = tab[j].d;
826: C(m) = tab[j].c; NEXT(m) = psum[j];
827: psum[j] = m;
828: }
829: }
830: }
831: for ( j = tlen-1, r = 0; j >= 0; j-- )
832: if ( psum[j] ) {
833: _MKDP(n,psum[j],t); _addd_destructive(vl,r,t,&t1); r = t1;
834: }
835: if ( r )
836: r->sugar = p->sugar + m0->dl->td;
837: *pr = r;
838: }
839: }
840:
841: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
842:
1.11 noro 843: void _weyl_mulmm_dup(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.9 noro 844: {
845: P c;
846: DL d,d0,d1,dt;
1.11 noro 847: int i,j,a,b,k,l,n2,s,min,curlen;
1.9 noro 848: struct cdl *p;
1.11 noro 849: static Q *ctab;
1.9 noro 850: static struct cdl *tab;
851: static int tablen;
852: static struct cdl *tmptab;
853: static int tmptablen;
854:
855: if ( !m0 || !m1 ) {
856: rtab[0].c = 0;
857: rtab[0].d = 0;
858: return;
859: }
860: mulp(vl,C(m0),C(m1),&c);
861: d0 = m0->dl; d1 = m1->dl;
862: n2 = n>>1;
863: curlen = 1;
864:
865: _NEWDL(d,n);
866: if ( n & 1 )
867: /* offset of h-degree */
868: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
869: else
870: d->td = 0;
871: rtab[0].c = c;
872: rtab[0].d = d;
873:
874: if ( rtablen > tmptablen ) {
875: if ( tmptab ) GC_free(tmptab);
876: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
877: tmptablen = rtablen;
878: }
879:
880: for ( i = 0; i < n2; i++ ) {
881: a = d0->d[i]; b = d1->d[n2+i];
882: k = d0->d[n2+i]; l = d1->d[i];
1.12 noro 883:
884: /* degree of xi^a*(Di^k*xi^l)*Di^b */
885: a += l;
886: b += k;
887: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
888:
1.9 noro 889: if ( !k || !l ) {
890: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
891: if ( p->c ) {
892: dt = p->d;
893: dt->d[i] = a;
894: dt->d[n2+i] = b;
895: dt->td += s;
896: }
897: }
898: curlen *= k+1;
899: continue;
900: }
901: if ( k+1 > tablen ) {
902: if ( tab ) GC_free(tab);
903: if ( ctab ) GC_free(ctab);
904: tablen = k+1;
905: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
1.11 noro 906: ctab = (Q *)MALLOC(tablen*sizeof(P));
1.9 noro 907: }
908: /* compute xi^a*(Di^k*xi^l)*Di^b */
909: min = MIN(k,l);
910: mkwc(k,l,ctab);
911: bzero(tab,(k+1)*sizeof(struct cdl));
912: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
913: if ( n & 1 )
914: for ( j = 0; j <= min; j++ ) {
915: _NEWDL(d,n);
1.12 noro 916: d->d[i] = a-j; d->d[n2+i] = b-j;
1.9 noro 917: d->td = s;
1.12 noro 918: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.9 noro 919: tab[j].d = d;
1.11 noro 920: tab[j].c = (P)ctab[j];
1.9 noro 921: }
922: else
923: for ( j = 0; j <= min; j++ ) {
924: _NEWDL(d,n);
1.12 noro 925: d->d[i] = a-j; d->d[n2+i] = b-j;
926: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.9 noro 927: tab[j].d = d;
1.11 noro 928: tab[j].c = (P)ctab[j];
1.9 noro 929: }
930: #if 0
931: _comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
932: for ( j = 0; j < curlen; j++ )
933: if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
934: for ( j = 0; j <= min; j++ )
935: if ( tab[j].d ) { _FREEDL(tab[j].d); }
936: curlen *= k+1;
937: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
938: #else
939: _comm_muld_tab_destructive(vl,n,rtab,curlen,tab,k+1);
940: for ( j = 0; j <= min; j++ )
941: if ( tab[j].d ) { _FREEDL(tab[j].d); }
942: curlen *= k+1;
943: #endif
944: }
945: }
946:
947: /* direct product of two cdl tables
948: rt[] = [
949: t[0]*t1[0],...,t[n-1]*t1[0],
950: t[0]*t1[1],...,t[n-1]*t1[1],
951: ...
952: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
953: ]
954: */
955:
1.11 noro 956: void _comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.9 noro 957: {
958: int i,j;
959: struct cdl *p;
960: P c;
961: DL d;
962:
963: bzero(rt,n*n1*sizeof(struct cdl));
964: for ( j = 0, p = rt; j < n1; j++ ) {
965: c = t1[j].c;
966: d = t1[j].d;
967: if ( !c )
968: break;
969: for ( i = 0; i < n; i++, p++ ) {
970: if ( t[i].c ) {
971: mulp(vl,t[i].c,c,&p->c);
972: _adddl_dup(nv,t[i].d,d,&p->d);
973: }
974: }
975: }
976: }
977:
1.11 noro 978: void _comm_muld_tab_destructive(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1)
1.9 noro 979: {
980: int i,j;
981: struct cdl *p;
982: P c;
983: DL d;
984:
985: for ( j = 1, p = t+n; j < n1; j++ ) {
986: c = t1[j].c;
987: d = t1[j].d;
988: if ( !c )
989: break;
990: for ( i = 0; i < n; i++, p++ ) {
991: if ( t[i].c ) {
992: mulp(vl,t[i].c,c,&p->c);
993: _adddl_dup(nv,t[i].d,d,&p->d);
994: }
995: }
996: }
997: c = t1[0].c;
998: d = t1[0].d;
999: for ( i = 0, p = t; i < n; i++, p++ )
1000: if ( t[i].c ) {
1001: mulp(vl,t[i].c,c,&p->c);
1002: /* t[i].d += d */
1003: adddl_destructive(nv,t[i].d,d);
1004: }
1.1 noro 1005: }