Annotation of OpenXM_contrib2/asir2000/engine/distm.c, Revision 1.21
1.5 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.6 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.5 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.21 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.20 2017/08/31 02:36:21 noro Exp $
1.5 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "inline.h"
52:
53: extern int (*cmpdl)();
1.2 noro 54: extern int do_weyl;
55:
1.9 noro 56: void ptomd(VL vl,int mod,VL dvl,P p,DP *pr)
1.1 noro 57: {
1.21 ! noro 58: P t;
1.1 noro 59:
1.21 ! noro 60: ptomp(mod,p,&t);
! 61: mptomd(vl,mod,dvl,t,pr);
1.1 noro 62: }
63:
1.9 noro 64: void mptomd(VL vl,int mod,VL dvl,P p,DP *pr)
1.1 noro 65: {
1.21 ! noro 66: int n,i;
! 67: VL tvl;
! 68: V v;
! 69: DL d;
! 70: MP m;
! 71: DCP dc;
! 72: DP r,s,t,u;
! 73: P x,c;
! 74:
! 75: if ( !p )
! 76: *pr = 0;
! 77: else {
! 78: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
! 79: if ( NUM(p) ) {
! 80: NEWDL(d,n);
! 81: NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr);
! 82: } else {
! 83: for ( i = 0, tvl = dvl, v = VR(p);
! 84: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
! 85: if ( !tvl ) {
! 86: for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
! 87: mptomd(vl,mod,dvl,COEF(dc),&t); pwrmp(vl,mod,x,DEG(dc),&c);
! 88: mulmdc(vl,mod,t,c,&r); addmd(vl,mod,r,s,&t); s = t;
! 89: }
! 90: *pr = s;
! 91: } else {
! 92: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
! 93: mptomd(vl,mod,dvl,COEF(dc),&t);
! 94: NEWDL(d,n); d->d[i] = QTOS(DEG(dc));
! 95: d->td = MUL_WEIGHT(d->d[i],i);
! 96: NEWMP(m); m->dl = d; C(m) = (Obj)ONEM; NEXT(m) = 0; MKDP(n,m,u);
! 97: comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t;
! 98: }
! 99: *pr = s;
! 100: }
! 101: }
! 102: }
1.18 noro 103: }
104:
105: void mdtodp(DP p,DP *pr)
106: {
1.21 ! noro 107: MP m,mr0,mr;
1.18 noro 108:
1.21 ! noro 109: if ( !p )
! 110: *pr = 0;
! 111: else {
! 112: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
! 113: NEXTMP(mr0,mr); mr->dl = m->dl;
! 114: mptop((P)C(m),(P *)&C(mr));
! 115: }
! 116: NEXT(mr) = 0;
! 117: MKDP(NV(p),mr0,*pr);
! 118: (*pr)->sugar = p->sugar;
! 119: }
1.18 noro 120: }
121:
122: void _mdtodp(DP p,DP *pr)
123: {
1.21 ! noro 124: MP m,mr0,mr;
! 125: int i;
! 126: Q q;
! 127:
! 128: if ( !p )
! 129: *pr = 0;
! 130: else {
! 131: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
! 132: NEXTMP(mr0,mr); mr->dl = m->dl;
! 133: i = ITOS(m->c); STOQ(i,q); C(mr) = (Obj)q;
! 134: }
! 135: NEXT(mr) = 0;
! 136: MKDP(NV(p),mr0,*pr);
! 137: (*pr)->sugar = p->sugar;
! 138: }
1.1 noro 139: }
140:
1.9 noro 141: void mdtop(VL vl,int mod,VL dvl,DP p,P *pr)
1.1 noro 142: {
1.21 ! noro 143: int n,i;
! 144: DL d;
! 145: MP m;
! 146: P r,s,t,u,w;
! 147: Q q;
! 148: VL tvl;
! 149:
! 150: if ( !p )
! 151: *pr = 0;
! 152: else {
! 153: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
! 154: for ( i = 0, t = (P)C(m), d = m->dl, tvl = dvl;
! 155: i < n; tvl = NEXT(tvl), i++ ) {
! 156: MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
! 157: mulmp(vl,mod,t,u,&w); t = w;
! 158: }
! 159: addmp(vl,mod,s,t,&u); s = u;
! 160: }
! 161: mptop(s,pr);
! 162: }
1.1 noro 163: }
164:
1.9 noro 165: void addmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.1 noro 166: {
1.21 ! noro 167: int n;
! 168: MP m1,m2,mr,mr0;
! 169: P t;
! 170: int tmp;
! 171: MQ q;
! 172:
! 173: if ( !p1 )
! 174: *pr = p2;
! 175: else if ( !p2 )
! 176: *pr = p1;
! 177: else {
! 178: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
! 179: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
! 180: case 0:
! 181: if ( NUM(C(m1)) && NUM(C(m2)) ) {
! 182: tmp = (CONT((MQ)C(m1))+CONT((MQ)C(m2))) - mod;
! 183: if ( tmp < 0 )
! 184: tmp += mod;
! 185: if ( tmp ) {
! 186: STOMQ(tmp,q); t = (P)q;
! 187: } else
! 188: t = 0;
! 189: } else
! 190: addmp(vl,mod,(P)C(m1),(P)C(m2),&t);
! 191: if ( t ) {
! 192: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = (Obj)t;
! 193: }
! 194: m1 = NEXT(m1); m2 = NEXT(m2); break;
! 195: case 1:
! 196: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
! 197: m1 = NEXT(m1); break;
! 198: case -1:
! 199: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
! 200: m2 = NEXT(m2); break;
! 201: }
! 202: if ( !mr0 )
! 203: if ( m1 )
! 204: mr0 = m1;
! 205: else if ( m2 )
! 206: mr0 = m2;
! 207: else {
! 208: *pr = 0;
! 209: return;
! 210: }
! 211: else if ( m1 )
! 212: NEXT(mr) = m1;
! 213: else if ( m2 )
! 214: NEXT(mr) = m2;
! 215: else
! 216: NEXT(mr) = 0;
! 217: MKDP(NV(p1),mr0,*pr);
! 218: if ( *pr )
! 219: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 220: }
1.1 noro 221: }
222:
1.9 noro 223: void submd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.1 noro 224: {
1.21 ! noro 225: DP t;
1.1 noro 226:
1.21 ! noro 227: if ( !p2 )
! 228: *pr = p1;
! 229: else {
! 230: chsgnmd(mod,p2,&t); addmd(vl,mod,p1,t,pr);
! 231: }
1.1 noro 232: }
233:
1.9 noro 234: void chsgnmd(int mod,DP p,DP *pr)
1.1 noro 235: {
1.21 ! noro 236: MP m,mr,mr0;
1.1 noro 237:
1.21 ! noro 238: if ( !p )
! 239: *pr = 0;
! 240: else {
! 241: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 242: NEXTMP(mr0,mr); chsgnmp(mod,(P)C(m),(P *)&C(mr)); mr->dl = m->dl;
! 243: }
! 244: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 245: if ( *pr )
! 246: (*pr)->sugar = p->sugar;
! 247: }
1.1 noro 248: }
249:
1.9 noro 250: void mulmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.1 noro 251: {
1.21 ! noro 252: if ( !do_weyl )
! 253: comm_mulmd(vl,mod,p1,p2,pr);
! 254: else
! 255: weyl_mulmd(vl,mod,p1,p2,pr);
1.2 noro 256: }
257:
1.9 noro 258: void comm_mulmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.2 noro 259: {
1.21 ! noro 260: MP m;
! 261: DP s,t,u;
! 262: int i,l,l1;
! 263: static MP *w;
! 264: static int wlen;
! 265:
! 266: if ( !p1 || !p2 )
! 267: *pr = 0;
! 268: else if ( OID(p1) <= O_P )
! 269: mulmdc(vl,mod,p2,(P)p1,pr);
! 270: else if ( OID(p2) <= O_P )
! 271: mulmdc(vl,mod,p1,(P)p2,pr);
! 272: else {
! 273: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
! 274: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 275: if ( l1 < l ) {
! 276: t = p1; p1 = p2; p2 = t;
! 277: l = l1;
! 278: }
! 279: if ( l > wlen ) {
! 280: if ( w ) GCFREE(w);
! 281: w = (MP *)MALLOC(l*sizeof(MP));
! 282: wlen = l;
! 283: }
! 284: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 285: w[i] = m;
! 286: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 287: mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
! 288: }
! 289: bzero(w,l*sizeof(MP));
! 290: *pr = s;
! 291: }
1.2 noro 292: }
293:
1.9 noro 294: void weyl_mulmd(VL vl,int mod,DP p1,DP p2,DP *pr)
1.2 noro 295: {
1.21 ! noro 296: MP m;
! 297: DP s,t,u;
! 298: int i,l;
! 299: static MP *w;
! 300: static int wlen;
! 301:
! 302: if ( !p1 || !p2 )
! 303: *pr = 0;
! 304: else if ( OID(p1) <= O_P )
! 305: mulmdc(vl,mod,p2,(P)p1,pr);
! 306: else if ( OID(p2) <= O_P )
! 307: mulmdc(vl,mod,p1,(P)p2,pr);
! 308: else {
! 309: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 310: if ( l > wlen ) {
! 311: if ( w ) GCFREE(w);
! 312: w = (MP *)MALLOC(l*sizeof(MP));
! 313: wlen = l;
! 314: }
! 315: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 316: w[i] = m;
! 317: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 318: weyl_mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
! 319: }
! 320: bzero(w,l*sizeof(MP));
! 321: *pr = s;
! 322: }
1.1 noro 323: }
324:
1.9 noro 325: void mulmdm(VL vl,int mod,DP p,MP m0,DP *pr)
1.1 noro 326: {
1.21 ! noro 327: MP m,mr,mr0;
! 328: P c;
! 329: DL d;
! 330: int n,t;
! 331: MQ q;
! 332:
! 333: if ( !p )
! 334: *pr = 0;
! 335: else {
! 336: for ( mr0 = 0, m = BDY(p), c = (P)C(m0), d = m0->dl, n = NV(p);
! 337: m; m = NEXT(m) ) {
! 338: NEXTMP(mr0,mr);
! 339: if ( NUM(C(m)) && NUM(c) ) {
! 340: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
! 341: STOMQ(t,q); C(mr) = (Obj)q;
! 342: } else
! 343: mulmp(vl,mod,(P)C(m),c,(P *)&C(mr));
! 344: adddl(n,m->dl,d,&mr->dl);
! 345: }
! 346: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 347: if ( *pr )
! 348: (*pr)->sugar = p->sugar + m0->dl->td;
! 349: }
1.1 noro 350: }
351:
1.9 noro 352: void weyl_mulmdm(VL vl,int mod,DP p,MP m0,DP *pr)
1.2 noro 353: {
1.21 ! noro 354: DP r,t,t1;
! 355: MP m;
! 356: int n,l,i;
! 357: static MP *w;
! 358: static int wlen;
! 359:
! 360: if ( !p )
! 361: *pr = 0;
! 362: else {
! 363: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 364: if ( l > wlen ) {
! 365: if ( w ) GCFREE(w);
! 366: w = (MP *)MALLOC(l*sizeof(MP));
! 367: wlen = l;
! 368: }
! 369: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 370: w[i] = m;
! 371: for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
! 372: weyl_mulmmm(vl,mod,w[i],m0,n,&t);
! 373: addmd(vl,mod,r,t,&t1); r = t1;
! 374: }
! 375: bzero(w,l*sizeof(MP));
! 376: if ( r )
! 377: r->sugar = p->sugar + m0->dl->td;
! 378: *pr = r;
! 379: }
1.2 noro 380: }
381:
382: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
383:
1.9 noro 384: void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *pr)
1.2 noro 385: {
1.21 ! noro 386: MP mr,mr0;
! 387: MQ mq;
! 388: DP r,t,t1;
! 389: P c,c0,c1;
! 390: DL d,d0,d1;
! 391: int i,j,a,b,k,l,n2,s,min;
! 392: static int *tab;
! 393: static int tablen;
! 394:
! 395: if ( !m0 || !m1 )
! 396: *pr = 0;
! 397: else {
! 398: c0 = (P)C(m0); c1 = (P)C(m1);
! 399: mulmp(vl,mod,c0,c1,&c);
! 400: d0 = m0->dl; d1 = m1->dl;
! 401: n2 = n>>1;
! 402: if ( n & 1 ) {
! 403: /* homogenized computation; dx-xd=h^2 */
! 404: /* offset of h-degree */
! 405: NEWDL(d,n);
! 406: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
! 407: NEWMP(mr); mr->c = (Obj)ONEM; mr->dl = d; NEXT(mr) = 0;
! 408: MKDP(n,mr,r); r->sugar = d->td;
! 409: } else
! 410: r = (DP)ONEM;
! 411: for ( i = 0; i < n2; i++ ) {
! 412: a = d0->d[i]; b = d1->d[n2+i];
! 413: k = d0->d[n2+i]; l = d1->d[i];
! 414:
! 415: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 416: a += l;
! 417: b += k;
! 418: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
! 419:
! 420: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 421: min = MIN(k,l);
! 422:
! 423: if ( min+1 > tablen ) {
! 424: if ( tab ) GCFREE(tab);
! 425: tab = (int *)MALLOC((min+1)*sizeof(int));
! 426: tablen = min+1;
! 427: }
! 428: mkwcm(k,l,mod,tab);
! 429: if ( n & 1 )
! 430: for ( mr0 = 0, j = 0; j <= min; j++ ) {
! 431: NEXTMP(mr0,mr); NEWDL(d,n);
! 432: d->d[i] = a-j; d->d[n2+i] = b-j;
! 433: d->td = s;
! 434: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
! 435: STOMQ(tab[j],mq); mr->c = (Obj)mq; mr->dl = d;
! 436: }
! 437: else
! 438: for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
! 439: NEXTMP(mr0,mr); NEWDL(d,n);
! 440: d->d[i] = a-j; d->d[n2+i] = b-j;
! 441: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
! 442: s = MAX(s,d->td); /* XXX */
! 443: STOMQ(tab[j],mq); mr->c = (Obj)mq; mr->dl = d;
! 444: }
! 445: bzero(tab,(min+1)*sizeof(int));
! 446: if ( mr0 )
! 447: NEXT(mr) = 0;
! 448: MKDP(n,mr0,t);
! 449: if ( t )
! 450: t->sugar = s;
! 451: comm_mulmd(vl,mod,r,t,&t1); r = t1;
! 452: }
! 453: mulmdc(vl,mod,r,c,pr);
! 454: }
1.2 noro 455: }
456:
1.9 noro 457: void mulmdc(VL vl,int mod,DP p,P c,DP *pr)
1.1 noro 458: {
1.21 ! noro 459: MP m,mr,mr0;
! 460: int t;
! 461: MQ q;
! 462:
! 463: if ( !p || !c )
! 464: *pr = 0;
! 465: else {
! 466: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 467: NEXTMP(mr0,mr);
! 468: if ( NUM(C(m)) && NUM(c) ) {
! 469: t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
! 470: STOMQ(t,q); C(mr) = (Obj)q;
! 471: } else
! 472: mulmp(vl,mod,(P)C(m),c,(P *)&C(mr));
! 473: mr->dl = m->dl;
! 474: }
! 475: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 476: if ( *pr )
! 477: (*pr)->sugar = p->sugar;
! 478: }
1.1 noro 479: }
480:
1.9 noro 481: void divsmdc(VL vl,int mod,DP p,P c,DP *pr)
1.1 noro 482: {
1.21 ! noro 483: MP m,mr,mr0;
1.1 noro 484:
1.21 ! noro 485: if ( !c )
! 486: error("disvsdc : division by 0");
! 487: else if ( !p )
! 488: *pr = 0;
! 489: else {
! 490: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 491: NEXTMP(mr0,mr); divsmp(vl,mod,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl;
! 492: }
! 493: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 494: if ( *pr )
! 495: (*pr)->sugar = p->sugar;
! 496: }
1.1 noro 497: }
498:
1.9 noro 499: void _dtop_mod(VL vl,VL dvl,DP p,P *pr)
1.1 noro 500: {
1.21 ! noro 501: int n,i;
! 502: DL d;
! 503: MP m;
! 504: P r,s,t,u,w;
! 505: Q q;
! 506: VL tvl;
! 507:
! 508: if ( !p )
! 509: *pr = 0;
! 510: else {
! 511: for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
! 512: i = ITOS(m->c); STOQ(i,q); t = (P)q;
! 513: for ( i = 0, d = m->dl, tvl = dvl;
! 514: i < n; tvl = NEXT(tvl), i++ ) {
! 515: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
! 516: mulp(vl,t,u,&w); t = w;
! 517: }
! 518: addp(vl,s,t,&u); s = u;
! 519: }
! 520: *pr = s;
! 521: }
1.1 noro 522: }
523:
1.9 noro 524: void _dp_mod(DP p,int mod,NODE subst,DP *rp)
1.1 noro 525: {
1.21 ! noro 526: MP m,mr,mr0;
! 527: P t,s;
! 528: NODE tn;
! 529:
! 530: if ( !p )
! 531: *rp = 0;
! 532: else {
! 533: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 534: for ( tn = subst, s = (P)m->c; tn; tn = NEXT(NEXT(tn)) ) {
! 535: substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t);
! 536: s = t;
! 537: }
! 538: ptomp(mod,s,&t);
! 539: if ( t ) {
! 540: NEXTMP(mr0,mr); mr->c = (Obj)STOI(CONT((MQ)t)); mr->dl = m->dl;
! 541: }
! 542: }
! 543: if ( mr0 ) {
! 544: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 545: } else
! 546: *rp = 0;
! 547: }
1.1 noro 548: }
549:
1.9 noro 550: void _dp_monic(DP p,int mod,DP *rp)
1.2 noro 551: {
1.21 ! noro 552: MP m,mr,mr0;
! 553: int c,c1;
1.2 noro 554:
1.21 ! noro 555: if ( !p )
! 556: *rp = 0;
! 557: else {
! 558: c = invm(ITOS(BDY(p)->c),mod);
! 559: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 560: c1 = dmar(ITOS(m->c),c,0,mod);
! 561: NEXTMP(mr0,mr); mr->c = (Obj)STOI(c1); mr->dl = m->dl;
! 562: }
! 563: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 564: }
1.1 noro 565: }
566:
1.9 noro 567: void _printdp(DP d)
1.1 noro 568: {
1.21 ! noro 569: int n,i;
! 570: MP m;
! 571: DL dl;
! 572:
! 573: if ( !d ) {
! 574: printf("0"); return;
! 575: }
! 576: for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) {
! 577: printf("%d*<<",ITOS(m->c));
! 578: for ( i = 0, dl = m->dl; i < n-1; i++ )
! 579: printf("%d,",dl->d[i]);
! 580: printf("%d",dl->d[i]);
! 581: printf(">>");
! 582: if ( NEXT(m) )
! 583: printf("+");
! 584: }
1.8 noro 585: }
586:
587: /* merge p1 and p2 into pr */
588:
1.9 noro 589: void addmd_destructive(int mod,DP p1,DP p2,DP *pr)
1.8 noro 590: {
1.21 ! noro 591: int n;
! 592: MP m1,m2,mr,mr0,s;
! 593: int t;
! 594:
! 595: if ( !p1 )
! 596: *pr = p2;
! 597: else if ( !p2 )
! 598: *pr = p1;
! 599: else {
! 600: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
! 601: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
! 602: case 0:
! 603: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
! 604: if ( t < 0 )
! 605: t += mod;
! 606: s = m1; m1 = NEXT(m1);
! 607: if ( t ) {
! 608: NEXTMP2(mr0,mr,s); C(mr) = (Obj)STOI(t);
! 609: }
! 610: s = m2; m2 = NEXT(m2);
! 611: break;
! 612: case 1:
! 613: s = m1; m1 = NEXT(m1); NEXTMP2(mr0,mr,s);
! 614: break;
! 615: case -1:
! 616: s = m2; m2 = NEXT(m2); NEXTMP2(mr0,mr,s);
! 617: break;
! 618: }
! 619: if ( !mr0 )
! 620: if ( m1 )
! 621: mr0 = m1;
! 622: else if ( m2 )
! 623: mr0 = m2;
! 624: else {
! 625: *pr = 0;
! 626: return;
! 627: }
! 628: else if ( m1 )
! 629: NEXT(mr) = m1;
! 630: else if ( m2 )
! 631: NEXT(mr) = m2;
! 632: else
! 633: NEXT(mr) = 0;
! 634: MKDP(NV(p1),mr0,*pr);
! 635: if ( *pr )
! 636: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 637: }
1.8 noro 638: }
639:
1.9 noro 640: void mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.8 noro 641: {
1.21 ! noro 642: if ( !do_weyl )
! 643: comm_mulmd_dup(mod,p1,p2,pr);
! 644: else
! 645: weyl_mulmd_dup(mod,p1,p2,pr);
1.8 noro 646: }
647:
1.9 noro 648: void comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.8 noro 649: {
1.21 ! noro 650: MP m;
! 651: DP s,t,u;
! 652: int i,l,l1;
! 653: static MP *w;
! 654: static int wlen;
! 655:
! 656: if ( !p1 || !p2 )
! 657: *pr = 0;
! 658: else {
! 659: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
! 660: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 661: if ( l1 < l ) {
! 662: t = p1; p1 = p2; p2 = t;
! 663: l = l1;
! 664: }
! 665: if ( l > wlen ) {
! 666: if ( w ) GCFREE(w);
! 667: w = (MP *)MALLOC(l*sizeof(MP));
! 668: wlen = l;
! 669: }
! 670: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 671: w[i] = m;
! 672: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 673: mulmdm_dup(mod,p1,w[i],&t); addmd_destructive(mod,s,t,&u); s = u;
! 674: }
! 675: bzero(w,l*sizeof(MP));
! 676: *pr = s;
! 677: }
1.8 noro 678: }
679:
1.9 noro 680: void weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.8 noro 681: {
1.21 ! noro 682: MP m;
! 683: DP s,t,u;
! 684: int i,l;
! 685: static MP *w;
! 686: static int wlen;
! 687:
! 688: if ( !p1 || !p2 )
! 689: *pr = 0;
! 690: else {
! 691: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
! 692: if ( l > wlen ) {
! 693: if ( w ) GCFREE(w);
! 694: w = (MP *)MALLOC(l*sizeof(MP));
! 695: wlen = l;
! 696: }
! 697: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
! 698: w[i] = m;
! 699: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 700: weyl_mulmdm_dup(mod,w[i],p2,&t); addmd_destructive(mod,s,t,&u); s = u;
! 701: }
! 702: bzero(w,l*sizeof(MP));
! 703: *pr = s;
! 704: }
1.8 noro 705: }
706:
1.9 noro 707: void mulmdm_dup(int mod,DP p,MP m0,DP *pr)
1.8 noro 708: {
1.21 ! noro 709: MP m,mr,mr0;
! 710: DL d,dt,dm;
! 711: int c,n,i;
! 712: int *pt,*p1,*p2;
! 713:
! 714: if ( !p )
! 715: *pr = 0;
! 716: else {
! 717: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
! 718: m; m = NEXT(m) ) {
! 719: NEXTMP(mr0,mr);
! 720: C(mr) = (Obj)STOI(dmar(ITOS(C(m)),c,0,mod));
! 721: NEWDL_NOINIT(dt,n); mr->dl = dt;
! 722: dm = m->dl;
! 723: dt->td = d->td + dm->td;
! 724: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
! 725: *pt++ = *p1++ + *p2++;
! 726: }
! 727: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 728: if ( *pr )
! 729: (*pr)->sugar = p->sugar + m0->dl->td;
! 730: }
1.8 noro 731: }
732:
1.9 noro 733: void weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr)
1.8 noro 734: {
1.21 ! noro 735: DP r,t,t1;
! 736: MP m;
! 737: DL d0;
! 738: int n,n2,l,i,j,tlen;
! 739: static MP *w,*psum;
! 740: static struct cdlm *tab;
! 741: static int wlen;
! 742: static int rtlen;
! 743:
! 744: if ( !p )
! 745: *pr = 0;
! 746: else {
! 747: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 748: if ( l > wlen ) {
! 749: if ( w ) GCFREE(w);
! 750: w = (MP *)MALLOC(l*sizeof(MP));
! 751: wlen = l;
! 752: }
! 753: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 754: w[i] = m;
! 755: n = NV(p); n2 = n>>1;
! 756: d0 = m0->dl;
! 757:
! 758: for ( i = 0, tlen = 1; i < n2; i++ )
! 759: tlen *= d0->d[n2+i]+1;
! 760: if ( tlen > rtlen ) {
! 761: if ( tab ) GCFREE(tab);
! 762: if ( psum ) GCFREE(psum);
! 763: rtlen = tlen;
! 764: tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
! 765: psum = (MP *)MALLOC(rtlen*sizeof(MP));
! 766: }
! 767: bzero(psum,tlen*sizeof(MP));
! 768: for ( i = l-1; i >= 0; i-- ) {
! 769: bzero(tab,tlen*sizeof(struct cdlm));
! 770: weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
! 771: for ( j = 0; j < tlen; j++ ) {
! 772: if ( tab[j].c ) {
! 773: NEWMP(m); m->dl = tab[j].d;
! 774: C(m) = (Obj)STOI(tab[j].c); NEXT(m) = psum[j];
! 775: psum[j] = m;
! 776: }
! 777: }
! 778: }
! 779: for ( j = tlen-1, r = 0; j >= 0; j-- )
! 780: if ( psum[j] ) {
! 781: MKDP(n,psum[j],t); addmd_destructive(mod,r,t,&t1); r = t1;
! 782: }
! 783: if ( r )
! 784: r->sugar = p->sugar + m0->dl->td;
! 785: *pr = r;
! 786: }
1.8 noro 787: }
788:
789: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
790:
1.9 noro 791: void weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct cdlm *rtab,int rtablen)
1.8 noro 792: {
1.21 ! noro 793: int c,c0,c1;
! 794: DL d,d0,d1,dt;
! 795: int i,j,a,b,k,l,n2,s,min,curlen;
! 796: struct cdlm *p;
! 797: static int *ctab;
! 798: static struct cdlm *tab;
! 799: static int tablen;
! 800: static struct cdlm *tmptab;
! 801: static int tmptablen;
! 802:
! 803: if ( !m0 || !m1 ) {
! 804: rtab[0].c = 0;
! 805: rtab[0].d = 0;
! 806: return;
! 807: }
! 808: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
! 809: c = dmar(c0,c1,0,mod);
! 810: d0 = m0->dl; d1 = m1->dl;
! 811: n2 = n>>1;
! 812: curlen = 1;
! 813:
! 814: NEWDL(d,n);
! 815: if ( n & 1 )
! 816: /* offset of h-degree */
! 817: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
! 818: else
! 819: d->td = 0;
! 820: rtab[0].c = c;
! 821: rtab[0].d = d;
! 822:
! 823: if ( rtablen > tmptablen ) {
! 824: if ( tmptab ) GCFREE(tmptab);
! 825: tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
! 826: tmptablen = rtablen;
! 827: }
! 828:
! 829: for ( i = 0; i < n2; i++ ) {
! 830: a = d0->d[i]; b = d1->d[n2+i];
! 831: k = d0->d[n2+i]; l = d1->d[i];
! 832:
! 833: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 834: a += l;
! 835: b += k;
! 836: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
! 837:
! 838: if ( !k || !l ) {
! 839: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
! 840: if ( p->c ) {
! 841: dt = p->d;
! 842: dt->d[i] = a;
! 843: dt->d[n2+i] = b;
! 844: dt->td += s;
! 845: }
! 846: }
! 847: curlen *= k+1;
! 848: continue;
! 849: }
! 850: if ( k+1 > tablen ) {
! 851: if ( tab ) GCFREE(tab);
! 852: if ( ctab ) GCFREE(ctab);
! 853: tablen = k+1;
! 854: tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
! 855: ctab = (int *)MALLOC(tablen*sizeof(int));
! 856: }
! 857: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 858: min = MIN(k,l);
! 859: mkwcm(k,l,mod,ctab);
! 860: bzero(tab,(k+1)*sizeof(struct cdlm));
! 861: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
! 862: if ( n & 1 )
! 863: for ( j = 0; j <= min; j++ ) {
! 864: NEWDL(d,n);
! 865: d->d[i] = a-j; d->d[n2+i] = b-j;
! 866: d->td = s;
! 867: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
! 868: tab[j].d = d;
! 869: tab[j].c = ctab[j];
! 870: }
! 871: else
! 872: for ( j = 0; j <= min; j++ ) {
! 873: NEWDL(d,n);
! 874: d->d[i] = a-j; d->d[n2+i] = b-j;
! 875: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
! 876: tab[j].d = d;
! 877: tab[j].c = ctab[j];
! 878: }
! 879: comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
! 880: curlen *= k+1;
! 881: }
1.8 noro 882: }
883:
1.9 noro 884: void comm_mulmd_tab_destructive(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1)
1.8 noro 885: {
1.21 ! noro 886: int i,j;
! 887: struct cdlm *p;
! 888: int c;
! 889: DL d;
! 890:
! 891: for ( j = 1, p = t+n; j < n1; j++ ) {
! 892: c = t1[j].c;
! 893: d = t1[j].d;
! 894: if ( !c )
! 895: break;
! 896: for ( i = 0; i < n; i++, p++ ) {
! 897: if ( t[i].c ) {
! 898: p->c = dmar(t[i].c,c,0,mod);
! 899: adddl_dup(nv,t[i].d,d,&p->d);
! 900: }
! 901: }
! 902: }
! 903: c = t1[0].c;
! 904: d = t1[0].d;
! 905: for ( i = 0, p = t; i < n; i++, p++ )
! 906: if ( t[i].c ) {
! 907: p->c = dmar(t[i].c,c,0,mod);
! 908: /* t[i].d += d */
! 909: adddl_destructive(nv,t[i].d,d);
! 910: }
1.8 noro 911: }
912:
1.9 noro 913: void adddl_dup(int n,DL d1,DL d2,DL *dr)
1.8 noro 914: {
1.21 ! noro 915: DL dt;
! 916: int i;
1.8 noro 917:
1.21 ! noro 918: NEWDL(dt,n);
! 919: *dr = dt;
! 920: dt->td = d1->td + d2->td;
! 921: for ( i = 0; i < n; i++ )
! 922: dt->d[i] = d1->d[i]+d2->d[i];
1.1 noro 923: }
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