Annotation of OpenXM_contrib2/asir2000/engine/_distm.c, Revision 1.17
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.17 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.16 2017/08/31 02:36:21 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: {
1.17 ! noro 65: _mp_free_list = 0;
! 66: _dp_free_list = 0;
! 67: _dl_free_list = 0;
! 68: GC_gcollect();
1.6 noro 69: }
1.1 noro 70:
71: void _DL_alloc()
72: {
1.17 ! noro 73: int *p;
! 74: int i,dl_len;
! 75: static int DL_alloc_count;
1.1 noro 76:
1.17 ! noro 77: /* fprintf(stderr,"DL_alloc : %d \n",++DL_alloc_count); */
! 78: dl_len = (current_dl_length+1);
1.14 ohara 79: #if SIZEOF_LONG == 8
1.17 ! noro 80: if ( dl_len & 1 )
! 81: dl_len += 1;
1.8 noro 82: #endif
1.17 ! noro 83: for ( i = 0; i < 128; i++, p += dl_len ) {
! 84: p = (int *)MALLOC(dl_len*sizeof(int));
! 85: *(DL *)p = _dl_free_list;
! 86: _dl_free_list = (DL)p;
! 87: }
1.1 noro 88: }
89:
90: void _MP_alloc()
91: {
1.17 ! noro 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++ ) {
! 98: p = (MP)MALLOC(sizeof(struct oMP));
! 99: p->next = _mp_free_list; _mp_free_list = p;
! 100: }
1.1 noro 101: }
1.17 ! noro 102:
1.1 noro 103: void _DP_alloc()
104: {
1.17 ! noro 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++ ) {
! 111: p = (DP)MALLOC(sizeof(struct oDP));
! 112: p->body = (MP)_dp_free_list; _dp_free_list = p;
! 113: }
1.1 noro 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: {
1.17 ! noro 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) = (Obj)STOI(t);
! 138: } else {
! 139: _FREEDL(s->dl); _FREEMP(s);
! 140: }
! 141: s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
! 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);
! 168: _FREEDP(p1); _FREEDP(p2);
! 169: }
1.1 noro 170: }
171:
1.11 noro 172: void _mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1 noro 173: {
1.17 ! noro 174: if ( !do_weyl )
! 175: _comm_mulmd_dup(mod,p1,p2,pr);
! 176: else
! 177: _weyl_mulmd_dup(mod,p1,p2,pr);
1.1 noro 178: }
179:
1.11 noro 180: void _comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1 noro 181: {
1.17 ! noro 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 ) GCFREE(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: }
1.1 noro 210: }
211:
1.11 noro 212: void _weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1 noro 213: {
1.17 ! noro 214: MP m;
! 215: DP s,t,u;
! 216: int i,l;
! 217: static MP *w;
! 218: static int wlen;
! 219:
! 220: if ( !p1 || !p2 )
! 221: *pr = 0;
! 222: else {
! 223: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
! 224: if ( l > wlen ) {
! 225: if ( w ) GCFREE(w);
! 226: w = (MP *)MALLOC(l*sizeof(MP));
! 227: wlen = l;
! 228: }
! 229: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
! 230: w[i] = m;
! 231: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 232: _weyl_mulmdm_dup(mod,w[i],p2,&t); _addmd_destructive(mod,s,t,&u); s = u;
! 233: }
! 234: bzero(w,l*sizeof(MP));
! 235: *pr = s;
! 236: }
1.1 noro 237: }
238:
1.11 noro 239: void _mulmdm_dup(int mod,DP p,MP m0,DP *pr)
1.1 noro 240: {
1.17 ! noro 241: MP m,mr,mr0;
! 242: DL d,dt,dm;
! 243: int c,n,i,c1,c2;
! 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);
! 252: c1 = ITOS(C(m));
! 253: DMAR(c1,c,0,mod,c2);
! 254: C(mr) = (Obj)STOI(c2);
! 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: }
1.1 noro 265: }
266:
1.11 noro 267: void _weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr)
1.1 noro 268: {
1.17 ! noro 269: DP r,t,t1;
! 270: MP m;
! 271: DL d0;
! 272: int n,n2,l,i,j,tlen;
! 273: static MP *w,*psum;
! 274: static struct cdlm *tab;
! 275: static int wlen;
! 276: static int rtlen;
! 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 ) GCFREE(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;
! 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 ) GCFREE(tab);
! 296: if ( psum ) GCFREE(psum);
! 297: rtlen = tlen;
! 298: tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
! 299: psum = (MP *)MALLOC(rtlen*sizeof(MP));
! 300: }
! 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) = (Obj)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: }
! 317: if ( r )
! 318: r->sugar = p->sugar + m0->dl->td;
! 319: *pr = r;
! 320: }
1.1 noro 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.17 ! noro 327: int c,c0,c1;
! 328: DL d,d0,d1,dt;
! 329: int i,j,a,b,k,l,n2,s,min,curlen;
! 330: struct cdlm *p;
! 331: static int *ctab;
! 332: static struct cdlm *tab;
! 333: static int tablen;
! 334: static struct cdlm *tmptab;
! 335: static int tmptablen;
! 336:
! 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 ) GCFREE(tmptab);
! 359: tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
! 360: tmptablen = rtablen;
! 361: }
! 362:
! 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:
! 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:
! 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 ) GCFREE(tab);
! 386: if ( ctab ) GCFREE(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 */
! 396: if ( n & 1 )
! 397: for ( j = 0; j <= min; j++ ) {
! 398: _NEWDL(d,n);
! 399: d->d[i] = a-j; d->d[n2+i] = b-j;
! 400: d->td = s;
! 401: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
! 402: tab[j].d = d;
! 403: tab[j].c = ctab[j];
! 404: }
! 405: else
! 406: for ( j = 0; j <= min; j++ ) {
! 407: _NEWDL(d,n);
! 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 */
! 410: tab[j].d = d;
! 411: tab[j].c = ctab[j];
! 412: }
1.4 noro 413: #if 0
1.17 ! noro 414: _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
! 415: for ( j = 0; j < curlen; j++ )
! 416: if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
! 417: for ( j = 0; j <= min; j++ )
! 418: if ( tab[j].d ) { _FREEDL(tab[j].d); }
! 419: curlen *= k+1;
! 420: bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
1.4 noro 421: #else
1.17 ! noro 422: _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
! 423: for ( j = 0; j <= min; j++ )
! 424: if ( tab[j].d ) { _FREEDL(tab[j].d); }
! 425: curlen *= k+1;
1.4 noro 426: #endif
1.17 ! noro 427: }
1.4 noro 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: {
1.17 ! noro 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);
! 455: _adddl_dup(nv,t[i].d,d,&p->d);
! 456: }
! 457: }
! 458: }
1.4 noro 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: {
1.17 ! noro 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);
! 476: _adddl_dup(nv,t[i].d,d,&p->d);
! 477: }
! 478: }
! 479: }
! 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: {
1.17 ! noro 492: int i,n;
! 493: DL t;
1.1 noro 494:
1.17 ! noro 495: n = current_dl_length;
! 496: _NEWDL(t,n); *dr = t;
! 497: t->td = d->td;
! 498: for ( i = 0; i < n; i++ )
! 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: {
1.17 ! noro 504: int i,n;
! 505: DL t;
1.1 noro 506:
1.17 ! noro 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];
1.1 noro 512: }
513:
1.11 noro 514: void _adddl_dup(int n,DL d1,DL d2,DL *dr)
1.1 noro 515: {
1.17 ! noro 516: DL dt;
! 517: int i;
1.1 noro 518:
1.17 ! noro 519: _NEWDL(dt,n);
! 520: *dr = dt;
! 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: {
1.17 ! noro 528: int i;
1.4 noro 529:
1.17 ! 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: {
1.17 ! noro 535: MP m,s;
1.1 noro 536:
1.17 ! noro 537: if ( !f )
! 538: return;
! 539: m = f->body;
! 540: while ( m ) {
! 541: s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
! 542: }
! 543: _FREEDP(f);
1.5 noro 544: }
545:
1.11 noro 546: void dpto_dp(DP p,DP *r)
1.5 noro 547: {
1.17 ! noro 548: MP m,mr0,mr;
! 549: DL t;
1.5 noro 550:
1.17 ! noro 551: if ( !p )
! 552: *r = 0;
! 553: else {
! 554: /* XXX : dummy call to set current_dl_length */
! 555: _NEWDL_NOINIT(t,NV(p));
! 556:
! 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;
! 565: }
1.1 noro 566: }
567:
1.11 noro 568: void _dptodp(DP p,DP *r)
1.1 noro 569: {
1.17 ! noro 570: MP m,mr0,mr;
1.1 noro 571:
1.17 ! noro 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: }
1.1 noro 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: {
1.17 ! noro 595: NODE top,prev,cur,m,t;
1.5 noro 596:
1.17 ! noro 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;
! 611: break;
! 612: }
! 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;
! 629: }
! 630: }
! 631: if ( !cur )
! 632: NEXT(prev) = m;
! 633: return top;
! 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: {
1.17 ! noro 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,(P)C(m1),(P)C(m2),&t);
! 654: s = m1; m1 = NEXT(m1);
! 655: if ( t ) {
! 656: _NEXTMP2(mr0,mr,s); C(mr) = (Obj)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: }
1.9 noro 689: }
690:
1.11 noro 691: void _muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9 noro 692: {
1.17 ! noro 693: if ( !do_weyl )
! 694: _comm_muld_dup(vl,p1,p2,pr);
! 695: else
! 696: _weyl_muld_dup(vl,p1,p2,pr);
1.9 noro 697: }
698:
1.11 noro 699: void _comm_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9 noro 700: {
1.17 ! noro 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 ) GCFREE(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: }
1.9 noro 729: }
730:
1.11 noro 731: void _weyl_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9 noro 732: {
1.17 ! noro 733: MP m;
! 734: DP s,t,u;
! 735: int i,l;
! 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 ) GCFREE(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: }
1.9 noro 756: }
757:
1.11 noro 758: void _muldm_dup(VL vl,DP p,MP m0,DP *pr)
1.9 noro 759: {
1.17 ! noro 760: MP m,mr,mr0;
! 761: DL d,dt,dm;
! 762: P c;
! 763: int n,i;
! 764: int *pt,*p1,*p2;
! 765:
! 766: if ( !p )
! 767: *pr = 0;
! 768: else {
! 769: for ( mr0 = 0, m = BDY(p), c = (P)C(m0), d = m0->dl, n = NV(p);
! 770: m; m = NEXT(m) ) {
! 771: _NEXTMP(mr0,mr);
! 772: mulp(vl,(P)C(m),c,(P *)&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: }
1.9 noro 783: }
784:
1.11 noro 785: void _weyl_muldm_dup(VL vl,MP m0,DP p,DP *pr)
1.9 noro 786: {
1.17 ! noro 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 ) GCFREE(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 ) GCFREE(tab);
! 814: if ( psum ) GCFREE(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: }
1.9 noro 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: {
1.17 ! noro 845: P c;
! 846: DL d,d0,d1,dt;
! 847: int i,j,a,b,k,l,n2,s,min,curlen;
! 848: struct cdl *p;
! 849: static Q *ctab;
! 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,(P)C(m0),(P)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 = (Obj)c;
! 872: rtab[0].d = d;
! 873:
! 874: if ( rtablen > tmptablen ) {
! 875: if ( tmptab ) GCFREE(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];
! 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:
! 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 ) GCFREE(tab);
! 903: if ( ctab ) GCFREE(ctab);
! 904: tablen = k+1;
! 905: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
! 906: ctab = (Q *)MALLOC(tablen*sizeof(P));
! 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);
! 916: d->d[i] = a-j; d->d[n2+i] = b-j;
! 917: d->td = s;
! 918: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
! 919: tab[j].d = d;
! 920: tab[j].c = (Obj)ctab[j];
! 921: }
! 922: else
! 923: for ( j = 0; j <= min; j++ ) {
! 924: _NEWDL(d,n);
! 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 */
! 927: tab[j].d = d;
! 928: tab[j].c = (Obj)ctab[j];
! 929: }
1.9 noro 930: #if 0
1.17 ! noro 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));
1.9 noro 938: #else
1.17 ! noro 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;
1.9 noro 943: #endif
1.17 ! noro 944: }
1.9 noro 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: {
1.17 ! noro 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 = (P)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,(P)t[i].c,c,(P *)&p->c);
! 972: _adddl_dup(nv,t[i].d,d,&p->d);
! 973: }
! 974: }
! 975: }
1.9 noro 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: {
1.17 ! noro 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 = (P)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,(P)t[i].c,c,(P *)&p->c);
! 993: _adddl_dup(nv,t[i].d,d,&p->d);
! 994: }
! 995: }
! 996: }
! 997: c = (P)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,(P)t[i].c,c,(P *)&p->c);
! 1002: /* t[i].d += d */
! 1003: adddl_destructive(nv,t[i].d,d);
! 1004: }
1.1 noro 1005: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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