Annotation of OpenXM_contrib2/asir2018/engine/_distm.c, Revision 1.1
1.1 ! 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
! 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
! 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: *
! 48: * $OpenXM$
! 49: */
! 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;
! 60:
! 61: void GC_gcollect();
! 62:
! 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: }
! 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);
! 79: #if SIZEOF_LONG == 8
! 80: if ( dl_len & 1 )
! 81: dl_len += 1;
! 82: #endif
! 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: }
! 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++ ) {
! 98: p = (MP)MALLOC(sizeof(struct oMP));
! 99: p->next = _mp_free_list; _mp_free_list = p;
! 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++ ) {
! 111: p = (DP)MALLOC(sizeof(struct oDP));
! 112: p->body = (MP)_dp_free_list; _dp_free_list = p;
! 113: }
! 114: }
! 115:
! 116: /* merge p1 and p2 into pr */
! 117:
! 118: void _addmd_destructive(int mod,DP p1,DP p2,DP *pr)
! 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) = (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: }
! 170: }
! 171:
! 172: void _mulmd_dup(int mod,DP p1,DP p2,DP *pr)
! 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:
! 180: void _comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
! 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 ) 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: }
! 210: }
! 211:
! 212: void _weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
! 213: {
! 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: }
! 237: }
! 238:
! 239: void _mulmdm_dup(int mod,DP p,MP m0,DP *pr)
! 240: {
! 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: }
! 265: }
! 266:
! 267: void _weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr)
! 268: {
! 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: }
! 321: }
! 322:
! 323: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
! 324:
! 325: void _weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct cdlm *rtab,int rtablen)
! 326: {
! 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: }
! 413: #if 0
! 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));
! 421: #else
! 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;
! 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:
! 439: void _comm_mulmd_tab(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1,struct cdlm *rt)
! 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);
! 455: _adddl_dup(nv,t[i].d,d,&p->d);
! 456: }
! 457: }
! 458: }
! 459: }
! 460:
! 461: void _comm_mulmd_tab_destructive(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1)
! 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);
! 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: }
! 488: }
! 489:
! 490: void dlto_dl(DL d,DL *dr)
! 491: {
! 492: int i,n;
! 493: DL t;
! 494:
! 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];
! 500: }
! 501:
! 502: void _dltodl(DL d,DL *dr)
! 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:
! 514: void _adddl_dup(int n,DL d1,DL d2,DL *dr)
! 515: {
! 516: DL dt;
! 517: int i;
! 518:
! 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];
! 524: }
! 525:
! 526: void _free_dlarray(DL *a,int n)
! 527: {
! 528: int i;
! 529:
! 530: for ( i = 0; i < n; i++ ) { _FREEDL(a[i]); }
! 531: }
! 532:
! 533: void _free_dp(DP f)
! 534: {
! 535: MP m,s;
! 536:
! 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);
! 544: }
! 545:
! 546: void dpto_dp(DP p,DP *r)
! 547: {
! 548: MP m,mr0,mr;
! 549: DL t;
! 550:
! 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: }
! 566: }
! 567:
! 568: void _dptodp(DP p,DP *r)
! 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:
! 586: /*
! 587: * destructive merge of two list
! 588: *
! 589: * p1, p2 : list of DL
! 590: * return : a merged list
! 591: */
! 592:
! 593: NODE _symb_merge(NODE m1,NODE m2,int n)
! 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;
! 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: }
! 635: }
! 636:
! 637: /* merge p1 and p2 into pr */
! 638:
! 639: void _addd_destructive(VL vl,DP p1,DP p2,DP *pr)
! 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,(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: }
! 689: }
! 690:
! 691: void _muld_dup(VL vl,DP p1,DP p2,DP *pr)
! 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:
! 699: void _comm_muld_dup(VL vl,DP p1,DP p2,DP *pr)
! 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 ) 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: }
! 729: }
! 730:
! 731: void _weyl_muld_dup(VL vl,DP p1,DP p2,DP *pr)
! 732: {
! 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: }
! 756: }
! 757:
! 758: void _muldm_dup(VL vl,DP p,MP m0,DP *pr)
! 759: {
! 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: }
! 783: }
! 784:
! 785: void _weyl_muldm_dup(VL vl,MP m0,DP p,DP *pr)
! 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 ) 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: }
! 839: }
! 840:
! 841: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
! 842:
! 843: void _weyl_mulmm_dup(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
! 844: {
! 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 Z *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 = (Z *)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: }
! 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:
! 956: void _comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
! 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 = (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: }
! 976: }
! 977:
! 978: void _comm_muld_tab_destructive(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1)
! 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 = (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: }
! 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 / asir2018 / engine