Annotation of OpenXM_contrib2/asir2018/engine/dist.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:
! 52: #define ORD_REVGRADLEX 0
! 53: #define ORD_GRADLEX 1
! 54: #define ORD_LEX 2
! 55: #define ORD_BREVGRADLEX 3
! 56: #define ORD_BGRADLEX 4
! 57: #define ORD_BLEX 5
! 58: #define ORD_BREVREV 6
! 59: #define ORD_BGRADREV 7
! 60: #define ORD_BLEXREV 8
! 61: #define ORD_ELIM 9
! 62: #define ORD_WEYL_ELIM 10
! 63: #define ORD_HOMO_WW_DRL 11
! 64: #define ORD_DRL_ZIGZAG 12
! 65: #define ORD_HOMO_WW_DRL_ZIGZAG 13
! 66:
! 67: int cmpdl_drl_zigzag(), cmpdl_homo_ww_drl_zigzag();
! 68: int cmpdl_top_weight();
! 69:
! 70: int (*cmpdl)()=cmpdl_revgradlex;
! 71: int (*cmpdl_tie_breaker)();
! 72: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
! 73:
! 74: Obj current_top_weight;
! 75: int current_top_weight_len;
! 76:
! 77: int do_weyl;
! 78:
! 79: int dp_nelim,dp_fcoeffs;
! 80: struct order_spec *dp_current_spec;
! 81: struct modorder_spec *dp_current_modspec;
! 82: int *dp_dl_work;
! 83:
! 84: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr);
! 85: void comm_quod(VL vl,DP p1,DP p2,DP *pr);
! 86: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr);
! 87: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr);
! 88: int create_order_spec(VL vl,Obj obj,struct order_spec **specp);
! 89: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s);
! 90:
! 91: void order_init()
! 92: {
! 93: struct order_spec *spec;
! 94:
! 95: create_order_spec(0,0,&spec);
! 96: initd(spec);
! 97: create_modorder_spec(0,0,&dp_current_modspec);
! 98: }
! 99:
! 100: int has_sfcoef_p(Obj f);
! 101:
! 102: int has_sfcoef(DP f)
! 103: {
! 104: MP t;
! 105:
! 106: if ( !f )
! 107: return 0;
! 108: for ( t = BDY(f); t; t = NEXT(t) )
! 109: if ( has_sfcoef_p(t->c) )
! 110: break;
! 111: return t ? 1 : 0;
! 112: }
! 113:
! 114: int has_sfcoef_p(Obj f)
! 115: {
! 116: DCP dc;
! 117:
! 118: if ( !f )
! 119: return 0;
! 120: else if ( NUM(f) )
! 121: return (NID((Num)f) == N_GFS) ? 1 : 0;
! 122: else if ( POLY(f) ) {
! 123: for ( dc = DC((P)f); dc; dc = NEXT(dc) )
! 124: if ( has_sfcoef_p((Obj)COEF(dc)) )
! 125: return 1;
! 126: return 0;
! 127: } else
! 128: return 0;
! 129: }
! 130:
! 131: extern Obj nd_top_weight;
! 132:
! 133: void reset_top_weight()
! 134: {
! 135: cmpdl = cmpdl_tie_breaker;
! 136: cmpdl_tie_breaker = 0;
! 137: nd_top_weight = 0;
! 138: current_top_weight = 0;
! 139: current_top_weight_len = 0;
! 140: }
! 141:
! 142: void initd(struct order_spec *spec)
! 143: {
! 144: int len,i,k,row;
! 145: Q **mat;
! 146:
! 147: switch ( spec->id ) {
! 148: case 3:
! 149: cmpdl = cmpdl_composite;
! 150: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
! 151: break;
! 152: case 2:
! 153: cmpdl = cmpdl_matrix;
! 154: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
! 155: break;
! 156: case 1:
! 157: cmpdl = cmpdl_order_pair;
! 158: break;
! 159: default:
! 160: switch ( spec->ord.simple ) {
! 161: case ORD_REVGRADLEX:
! 162: cmpdl = cmpdl_revgradlex; break;
! 163: case ORD_GRADLEX:
! 164: cmpdl = cmpdl_gradlex; break;
! 165: case ORD_BREVGRADLEX:
! 166: cmpdl = cmpdl_brevgradlex; break;
! 167: case ORD_BGRADLEX:
! 168: cmpdl = cmpdl_bgradlex; break;
! 169: case ORD_BLEX:
! 170: cmpdl = cmpdl_blex; break;
! 171: case ORD_BREVREV:
! 172: cmpdl = cmpdl_brevrev; break;
! 173: case ORD_BGRADREV:
! 174: cmpdl = cmpdl_bgradrev; break;
! 175: case ORD_BLEXREV:
! 176: cmpdl = cmpdl_blexrev; break;
! 177: case ORD_ELIM:
! 178: cmpdl = cmpdl_elim; break;
! 179: case ORD_WEYL_ELIM:
! 180: cmpdl = cmpdl_weyl_elim; break;
! 181: case ORD_HOMO_WW_DRL:
! 182: cmpdl = cmpdl_homo_ww_drl; break;
! 183: case ORD_DRL_ZIGZAG:
! 184: cmpdl = cmpdl_drl_zigzag; break;
! 185: case ORD_HOMO_WW_DRL_ZIGZAG:
! 186: cmpdl = cmpdl_homo_ww_drl_zigzag; break;
! 187: case ORD_LEX: default:
! 188: cmpdl = cmpdl_lex; break;
! 189: }
! 190: break;
! 191: }
! 192: if ( current_top_weight ) {
! 193: cmpdl_tie_breaker = cmpdl;
! 194: cmpdl = cmpdl_top_weight;
! 195: if ( OID(current_top_weight) == O_VECT ) {
! 196: mat = (Q **)&BDY((VECT)current_top_weight);
! 197: row = 1;
! 198: } else {
! 199: mat = (Q **)BDY((MAT)current_top_weight);
! 200: row = ((MAT)current_top_weight)->row;
! 201: }
! 202: for ( k = 0, len = 0; k < row; k++ )
! 203: for ( i = 0; i < spec->nv; i++ )
! 204: if ( mat[k][i] )
! 205: len = MAX(mpz_size(BDY((Z)mat[k][i])),len);
! 206: current_top_weight_len = len;
! 207: }
! 208: dp_current_spec = spec;
! 209: }
! 210:
! 211: int dpm_ispot;
! 212:
! 213: /* type=0 => TOP, type=1 => POT */
! 214: void initdpm(struct order_spec *spec,int type)
! 215: {
! 216: int len,i,k,row;
! 217: Q **mat;
! 218:
! 219: initd(spec);
! 220: dpm_ispot = type;
! 221: }
! 222:
! 223: void ptod(VL vl,VL dvl,P p,DP *pr)
! 224: {
! 225: int n,i,j,k;
! 226: VL tvl;
! 227: V v;
! 228: DL d;
! 229: MP m;
! 230: DCP dc;
! 231: DCP *w;
! 232: DP r,s,t,u;
! 233: P x,c;
! 234:
! 235: if ( !p )
! 236: *pr = 0;
! 237: else if ( OID(p) > O_P )
! 238: error("ptod : only polynomials can be converted.");
! 239: else {
! 240: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
! 241: if ( NUM(p) ) {
! 242: NEWDL(d,n);
! 243: NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
! 244: } else {
! 245: for ( i = 0, tvl = dvl, v = VR(p);
! 246: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
! 247: if ( !tvl ) {
! 248: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
! 249: w = (DCP *)ALLOCA(k*sizeof(DCP));
! 250: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
! 251: w[j] = dc;
! 252:
! 253: for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
! 254: ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
! 255: muldc(vl,t,(Obj)c,&r); addd(vl,r,s,&t); s = t;
! 256: }
! 257: *pr = s;
! 258: } else {
! 259: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
! 260: w = (DCP *)ALLOCA(k*sizeof(DCP));
! 261: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
! 262: w[j] = dc;
! 263:
! 264: for ( j = k-1, s = 0; j >= 0; j-- ) {
! 265: ptod(vl,dvl,COEF(w[j]),&t);
! 266: NEWDL(d,n); d->d[i] = QTOS(DEG(w[j]));
! 267: d->td = MUL_WEIGHT(d->d[i],i);
! 268: NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
! 269: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
! 270: }
! 271: *pr = s;
! 272: }
! 273: }
! 274: }
! 275: #if 0
! 276: if ( !dp_fcoeffs && has_sfcoef(*pr) )
! 277: dp_fcoeffs = N_GFS;
! 278: #endif
! 279: }
! 280:
! 281: void dtop(VL vl,VL dvl,DP p,Obj *pr)
! 282: {
! 283: int n,i,j,k;
! 284: DL d;
! 285: MP m;
! 286: MP *a;
! 287: P r;
! 288: Obj t,w,s,u;
! 289: Z q;
! 290: VL tvl;
! 291:
! 292: if ( !p )
! 293: *pr = 0;
! 294: else {
! 295: for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
! 296: a = (MP *)ALLOCA(k*sizeof(MP));
! 297: for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
! 298: a[j] = m;
! 299:
! 300: for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
! 301: m = a[j];
! 302: t = C(m);
! 303: if ( NUM(t) && NID((Num)t) == N_M ) {
! 304: mptop((P)t,(P *)&u); t = u;
! 305: }
! 306: for ( i = 0, d = m->dl, tvl = dvl;
! 307: i < n; tvl = NEXT(tvl), i++ ) {
! 308: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,(P *)&u);
! 309: arf_mul(vl,t,(Obj)u,&w); t = w;
! 310: }
! 311: arf_add(vl,s,t,&u); s = u;
! 312: }
! 313: *pr = s;
! 314: }
! 315: }
! 316:
! 317: void nodetod(NODE node,DP *dp)
! 318: {
! 319: NODE t;
! 320: int len,i,td;
! 321: Q e;
! 322: DL d;
! 323: MP m;
! 324: DP u;
! 325:
! 326: for ( t = node, len = 0; t; t = NEXT(t), len++ );
! 327: NEWDL(d,len);
! 328: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
! 329: e = (Q)BDY(t);
! 330: if ( !e )
! 331: d->d[i] = 0;
! 332: else if ( !NUM(e) || !RATN(e) || !INT(e) )
! 333: error("nodetod : invalid input");
! 334: else {
! 335: d->d[i] = QTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
! 336: }
! 337: }
! 338: d->td = td;
! 339: NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0;
! 340: MKDP(len,m,u); u->sugar = td; *dp = u;
! 341: }
! 342:
! 343: void nodetodpm(NODE node,Obj pos,DPM *dp)
! 344: {
! 345: NODE t;
! 346: int len,i,td;
! 347: Q e;
! 348: DL d;
! 349: DMM m;
! 350: DPM u;
! 351:
! 352: for ( t = node, len = 0; t; t = NEXT(t), len++ );
! 353: NEWDL(d,len);
! 354: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
! 355: e = (Q)BDY(t);
! 356: if ( !e )
! 357: d->d[i] = 0;
! 358: else if ( !NUM(e) || !RATN(e) || !INT(e) )
! 359: error("nodetodpm : invalid input");
! 360: else {
! 361: d->d[i] = QTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
! 362: }
! 363: }
! 364: d->td = td;
! 365: NEWDMM(m); m->dl = d; m->pos = QTOS((Q)pos); C(m) = (Obj)ONE; NEXT(m) = 0;
! 366: MKDPM(len,m,u); u->sugar = td; *dp = u;
! 367: }
! 368:
! 369: void dtodpm(DP d,int pos,DPM *dp)
! 370: {
! 371: DMM mr0,mr;
! 372: MP m;
! 373:
! 374: if ( !d ) *dp = 0;
! 375: else {
! 376: for ( m = BDY(d), mr0 = 0; m; m = NEXT(m) ) {
! 377: NEXTDMM(mr0,mr);
! 378: mr->dl = m->dl;
! 379: mr->pos = pos;
! 380: C(mr) = C(m);
! 381: }
! 382: MKDPM(d->nv,mr0,*dp); (*dp)->sugar = d->sugar;
! 383: }
! 384: }
! 385:
! 386: int sugard(MP m)
! 387: {
! 388: int s;
! 389:
! 390: for ( s = 0; m; m = NEXT(m) )
! 391: s = MAX(s,m->dl->td);
! 392: return s;
! 393: }
! 394:
! 395: void addd(VL vl,DP p1,DP p2,DP *pr)
! 396: {
! 397: int n;
! 398: MP m1,m2,mr=0,mr0;
! 399: Obj t;
! 400: DL d;
! 401:
! 402: if ( !p1 )
! 403: *pr = p2;
! 404: else if ( !p2 )
! 405: *pr = p1;
! 406: else {
! 407: if ( OID(p1) <= O_R ) {
! 408: n = NV(p2); NEWDL(d,n);
! 409: NEWMP(m1); m1->dl = d; C(m1) = (Obj)p1; NEXT(m1) = 0;
! 410: MKDP(n,m1,p1); (p1)->sugar = 0;
! 411: }
! 412: if ( OID(p2) <= O_R ) {
! 413: n = NV(p1); NEWDL(d,n);
! 414: NEWMP(m2); m2->dl = d; C(m2) = (Obj)p2; NEXT(m2) = 0;
! 415: MKDP(n,m2,p2); (p2)->sugar = 0;
! 416: }
! 417: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
! 418: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
! 419: case 0:
! 420: arf_add(vl,C(m1),C(m2),&t);
! 421: if ( t ) {
! 422: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
! 423: }
! 424: m1 = NEXT(m1); m2 = NEXT(m2); break;
! 425: case 1:
! 426: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
! 427: m1 = NEXT(m1); break;
! 428: case -1:
! 429: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
! 430: m2 = NEXT(m2); break;
! 431: }
! 432: if ( !mr0 )
! 433: if ( m1 )
! 434: mr0 = m1;
! 435: else if ( m2 )
! 436: mr0 = m2;
! 437: else {
! 438: *pr = 0;
! 439: return;
! 440: }
! 441: else if ( m1 )
! 442: NEXT(mr) = m1;
! 443: else if ( m2 )
! 444: NEXT(mr) = m2;
! 445: else
! 446: NEXT(mr) = 0;
! 447: MKDP(NV(p1),mr0,*pr);
! 448: if ( *pr )
! 449: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 450: }
! 451: }
! 452:
! 453: /* for F4 symbolic reduction */
! 454:
! 455: void symb_addd(DP p1,DP p2,DP *pr)
! 456: {
! 457: int n;
! 458: MP m1,m2,mr=0,mr0;
! 459:
! 460: if ( !p1 )
! 461: *pr = p2;
! 462: else if ( !p2 )
! 463: *pr = p1;
! 464: else {
! 465: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
! 466: NEXTMP(mr0,mr); C(mr) = (Obj)ONE;
! 467: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
! 468: case 0:
! 469: mr->dl = m1->dl;
! 470: m1 = NEXT(m1); m2 = NEXT(m2); break;
! 471: case 1:
! 472: mr->dl = m1->dl;
! 473: m1 = NEXT(m1); break;
! 474: case -1:
! 475: mr->dl = m2->dl;
! 476: m2 = NEXT(m2); break;
! 477: }
! 478: }
! 479: if ( !mr0 )
! 480: if ( m1 )
! 481: mr0 = m1;
! 482: else if ( m2 )
! 483: mr0 = m2;
! 484: else {
! 485: *pr = 0;
! 486: return;
! 487: }
! 488: else if ( m1 )
! 489: NEXT(mr) = m1;
! 490: else if ( m2 )
! 491: NEXT(mr) = m2;
! 492: else
! 493: NEXT(mr) = 0;
! 494: MKDP(NV(p1),mr0,*pr);
! 495: if ( *pr )
! 496: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 497: }
! 498: }
! 499:
! 500: /*
! 501: * destructive merge of two list
! 502: *
! 503: * p1, p2 : list of DL
! 504: * return : a merged list
! 505: */
! 506:
! 507: NODE symb_merge(NODE m1,NODE m2,int n)
! 508: {
! 509: NODE top=0,prev,cur,m=0,t;
! 510: DL d1,d2;
! 511:
! 512: if ( !m1 )
! 513: return m2;
! 514: else if ( !m2 )
! 515: return m1;
! 516: else {
! 517: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
! 518: case 0:
! 519: top = m1; m = NEXT(m2);
! 520: break;
! 521: case 1:
! 522: top = m1; m = m2;
! 523: break;
! 524: case -1:
! 525: top = m2; m = m1;
! 526: break;
! 527: }
! 528: prev = top; cur = NEXT(top);
! 529: /* BDY(prev) > BDY(m) always holds */
! 530: while ( cur && m ) {
! 531: d1 = (DL)BDY(cur);
! 532: d2 = (DL)BDY(m);
! 533: #if 1
! 534: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
! 535: #else
! 536: /* XXX only valid for DRL */
! 537: if ( d1->td > d2->td )
! 538: c = 1;
! 539: else if ( d1->td < d2->td )
! 540: c = -1;
! 541: else {
! 542: for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
! 543: if ( i < 0 )
! 544: c = 0;
! 545: else if ( d1->d[i] < d2->d[i] )
! 546: c = 1;
! 547: else
! 548: c = -1;
! 549: }
! 550: switch ( c ) {
! 551: #endif
! 552: case 0:
! 553: m = NEXT(m);
! 554: prev = cur; cur = NEXT(cur);
! 555: break;
! 556: case 1:
! 557: t = NEXT(cur); NEXT(cur) = m; m = t;
! 558: prev = cur; cur = NEXT(cur);
! 559: break;
! 560: case -1:
! 561: NEXT(prev) = m; m = cur;
! 562: prev = NEXT(prev); cur = NEXT(prev);
! 563: break;
! 564: }
! 565: }
! 566: if ( !cur )
! 567: NEXT(prev) = m;
! 568: return top;
! 569: }
! 570: }
! 571:
! 572: void _adddl(int n,DL d1,DL d2,DL d3)
! 573: {
! 574: int i;
! 575:
! 576: d3->td = d1->td+d2->td;
! 577: for ( i = 0; i < n; i++ )
! 578: d3->d[i] = d1->d[i]+d2->d[i];
! 579: }
! 580:
! 581: /* m1 <- m1 U dl*f, destructive */
! 582:
! 583: NODE mul_dllist(DL dl,DP f);
! 584:
! 585: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
! 586: {
! 587: NODE top,prev,cur,n1;
! 588: DP g;
! 589: DL t,s;
! 590: MP m;
! 591:
! 592: if ( !m1 )
! 593: return mul_dllist(dl,f);
! 594: else if ( !f )
! 595: return m1;
! 596: else {
! 597: m = BDY(f);
! 598: NEWDL_NOINIT(t,n);
! 599: _adddl(n,m->dl,dl,t);
! 600: top = m1; prev = 0; cur = m1;
! 601: while ( m ) {
! 602: switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
! 603: case 0:
! 604: prev = cur; cur = NEXT(cur);
! 605: if ( !cur ) {
! 606: MKDP(n,m,g);
! 607: NEXT(prev) = mul_dllist(dl,g);
! 608: return top;
! 609: }
! 610: m = NEXT(m);
! 611: if ( m ) _adddl(n,m->dl,dl,t);
! 612: break;
! 613: case 1:
! 614: prev = cur; cur = NEXT(cur);
! 615: if ( !cur ) {
! 616: MKDP(n,m,g);
! 617: NEXT(prev) = mul_dllist(dl,g);
! 618: return top;
! 619: }
! 620: break;
! 621: case -1:
! 622: NEWDL_NOINIT(s,n);
! 623: s->td = t->td;
! 624: bcopy(t->d,s->d,n*sizeof(int));
! 625: NEWNODE(n1);
! 626: n1->body = (pointer)s;
! 627: NEXT(n1) = cur;
! 628: if ( !prev ) {
! 629: top = n1; cur = n1;
! 630: } else {
! 631: NEXT(prev) = n1; prev = n1;
! 632: }
! 633: m = NEXT(m);
! 634: if ( m ) _adddl(n,m->dl,dl,t);
! 635: break;
! 636: }
! 637: }
! 638: return top;
! 639: }
! 640: }
! 641:
! 642: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
! 643: {
! 644: DLBUCKET top,prev,cur,m,t;
! 645:
! 646: if ( !m1 )
! 647: return m2;
! 648: else if ( !m2 )
! 649: return m1;
! 650: else {
! 651: if ( m1->td == m2->td ) {
! 652: top = m1;
! 653: BDY(top) = symb_merge(BDY(top),BDY(m2),n);
! 654: m = NEXT(m2);
! 655: } else if ( m1->td > m2->td ) {
! 656: top = m1; m = m2;
! 657: } else {
! 658: top = m2; m = m1;
! 659: }
! 660: prev = top; cur = NEXT(top);
! 661: /* prev->td > m->td always holds */
! 662: while ( cur && m ) {
! 663: if ( cur->td == m->td ) {
! 664: BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
! 665: m = NEXT(m);
! 666: prev = cur; cur = NEXT(cur);
! 667: } else if ( cur->td > m->td ) {
! 668: t = NEXT(cur); NEXT(cur) = m; m = t;
! 669: prev = cur; cur = NEXT(cur);
! 670: } else {
! 671: NEXT(prev) = m; m = cur;
! 672: prev = NEXT(prev); cur = NEXT(prev);
! 673: }
! 674: }
! 675: if ( !cur )
! 676: NEXT(prev) = m;
! 677: return top;
! 678: }
! 679: }
! 680:
! 681: void subd(VL vl,DP p1,DP p2,DP *pr)
! 682: {
! 683: DP t;
! 684:
! 685: if ( !p2 )
! 686: *pr = p1;
! 687: else {
! 688: chsgnd(p2,&t); addd(vl,p1,t,pr);
! 689: }
! 690: }
! 691:
! 692: void chsgnd(DP p,DP *pr)
! 693: {
! 694: MP m,mr=0,mr0;
! 695: Obj r;
! 696:
! 697: if ( !p )
! 698: *pr = 0;
! 699: else if ( OID(p) <= O_R ) {
! 700: arf_chsgn((Obj)p,&r); *pr = (DP)r;
! 701: } else {
! 702: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 703: NEXTMP(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->dl = m->dl;
! 704: }
! 705: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 706: if ( *pr )
! 707: (*pr)->sugar = p->sugar;
! 708: }
! 709: }
! 710:
! 711: void muld(VL vl,DP p1,DP p2,DP *pr)
! 712: {
! 713: if ( ! do_weyl )
! 714: comm_muld(vl,p1,p2,pr);
! 715: else
! 716: weyl_muld(vl,p1,p2,pr);
! 717: }
! 718:
! 719: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
! 720: {
! 721: MP m;
! 722: DP s,t,u;
! 723: int i,l,l1;
! 724: static MP *w;
! 725: static int wlen;
! 726:
! 727: if ( !p1 || !p2 )
! 728: *pr = 0;
! 729: else if ( OID(p1) != O_DP )
! 730: muldc(vl,p2,(Obj)p1,pr);
! 731: else if ( OID(p2) != O_DP )
! 732: muldc(vl,p1,(Obj)p2,pr);
! 733: else {
! 734: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
! 735: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 736: if ( l1 < l ) {
! 737: t = p1; p1 = p2; p2 = t;
! 738: l = l1;
! 739: }
! 740: if ( l > wlen ) {
! 741: if ( w ) GCFREE(w);
! 742: w = (MP *)MALLOC(l*sizeof(MP));
! 743: wlen = l;
! 744: }
! 745: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 746: w[i] = m;
! 747: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 748: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
! 749: }
! 750: bzero(w,l*sizeof(MP));
! 751: *pr = s;
! 752: }
! 753: }
! 754:
! 755: /* discard terms which is not a multiple of dl */
! 756:
! 757: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
! 758: {
! 759: MP m;
! 760: DP s,t,u;
! 761: int i,l,l1;
! 762: static MP *w;
! 763: static int wlen;
! 764:
! 765: if ( !p1 || !p2 )
! 766: *pr = 0;
! 767: else if ( OID(p1) != O_DP )
! 768: muldc_trunc(vl,p2,(Obj)p1,dl,pr);
! 769: else if ( OID(p2) != O_DP )
! 770: muldc_trunc(vl,p1,(Obj)p2,dl,pr);
! 771: else {
! 772: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
! 773: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
! 774: if ( l1 < l ) {
! 775: t = p1; p1 = p2; p2 = t;
! 776: l = l1;
! 777: }
! 778: if ( l > wlen ) {
! 779: if ( w ) GCFREE(w);
! 780: w = (MP *)MALLOC(l*sizeof(MP));
! 781: wlen = l;
! 782: }
! 783: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
! 784: w[i] = m;
! 785: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 786: muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
! 787: }
! 788: bzero(w,l*sizeof(MP));
! 789: *pr = s;
! 790: }
! 791: }
! 792:
! 793: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
! 794: {
! 795: MP m=0,m0;
! 796: DP s,t;
! 797: int i,n,sugar;
! 798: DL d1,d2,d;
! 799: Q a,b;
! 800:
! 801: if ( !p2 )
! 802: error("comm_quod : invalid input");
! 803: if ( !p1 )
! 804: *pr = 0;
! 805: else {
! 806: n = NV(p1);
! 807: d2 = BDY(p2)->dl;
! 808: m0 = 0;
! 809: sugar = p1->sugar;
! 810: while ( p1 ) {
! 811: d1 = BDY(p1)->dl;
! 812: NEWDL(d,n);
! 813: d->td = d1->td - d2->td;
! 814: for ( i = 0; i < n; i++ )
! 815: d->d[i] = d1->d[i]-d2->d[i];
! 816: NEXTMP(m0,m);
! 817: m->dl = d;
! 818: divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
! 819: C(m) = (Obj)b;
! 820: muldm_trunc(vl,p2,m,d2,&t);
! 821: addd(vl,p1,t,&s); p1 = s;
! 822: C(m) = (Obj)a;
! 823: }
! 824: if ( m0 ) {
! 825: NEXT(m) = 0; MKDP(n,m0,*pr);
! 826: } else
! 827: *pr = 0;
! 828: /* XXX */
! 829: if ( *pr )
! 830: (*pr)->sugar = sugar - d2->td;
! 831: }
! 832: }
! 833:
! 834: void muldm(VL vl,DP p,MP m0,DP *pr)
! 835: {
! 836: MP m,mr=0,mr0;
! 837: Obj c;
! 838: DL d;
! 839: int n;
! 840:
! 841: if ( !p )
! 842: *pr = 0;
! 843: else {
! 844: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
! 845: m; m = NEXT(m) ) {
! 846: NEXTMP(mr0,mr);
! 847: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
! 848: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
! 849: else
! 850: arf_mul(vl,C(m),c,&C(mr));
! 851: adddl(n,m->dl,d,&mr->dl);
! 852: }
! 853: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 854: if ( *pr )
! 855: (*pr)->sugar = p->sugar + m0->dl->td;
! 856: }
! 857: }
! 858:
! 859: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
! 860: {
! 861: MP m,mr=0,mr0;
! 862: Obj c;
! 863: DL d,tdl;
! 864: int n,i;
! 865:
! 866: if ( !p )
! 867: *pr = 0;
! 868: else {
! 869: n = NV(p);
! 870: NEWDL(tdl,n);
! 871: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
! 872: m; m = NEXT(m) ) {
! 873: _adddl(n,m->dl,d,tdl);
! 874: for ( i = 0; i < n; i++ )
! 875: if ( tdl->d[i] < dl->d[i] )
! 876: break;
! 877: if ( i < n )
! 878: continue;
! 879: NEXTMP(mr0,mr);
! 880: mr->dl = tdl;
! 881: NEWDL(tdl,n);
! 882: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
! 883: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
! 884: else
! 885: arf_mul(vl,C(m),(Obj)c,&C(mr));
! 886: }
! 887: if ( mr0 ) {
! 888: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 889: } else
! 890: *pr = 0;
! 891: if ( *pr )
! 892: (*pr)->sugar = p->sugar + m0->dl->td;
! 893: }
! 894: }
! 895:
! 896: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
! 897: {
! 898: MP m;
! 899: DP s,t,u;
! 900: int i,l;
! 901: static MP *w;
! 902: static int wlen;
! 903:
! 904: if ( !p1 || !p2 )
! 905: *pr = 0;
! 906: else if ( OID(p1) != O_DP )
! 907: muldc(vl,p2,(Obj)p1,pr);
! 908: else if ( OID(p2) != O_DP )
! 909: muldc(vl,p1,(Obj)p2,pr);
! 910: else {
! 911: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
! 912: if ( l > wlen ) {
! 913: if ( w ) GCFREE(w);
! 914: w = (MP *)MALLOC(l*sizeof(MP));
! 915: wlen = l;
! 916: }
! 917: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
! 918: w[i] = m;
! 919: for ( s = 0, i = l-1; i >= 0; i-- ) {
! 920: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
! 921: }
! 922: bzero(w,l*sizeof(MP));
! 923: *pr = s;
! 924: }
! 925: }
! 926:
! 927: void actm(VL vl,int nv,MP m1,MP m2,DP *pr)
! 928: {
! 929: DL d1,d2,d;
! 930: int n2,i,j,k;
! 931: Z jq,c,c1;
! 932: MP m;
! 933: Obj t;
! 934:
! 935: d1 = m1->dl;
! 936: d2 = m2->dl;
! 937: for ( i = 0; i < nv; i++ )
! 938: if ( d1->d[i] > d2->d[i] ) {
! 939: *pr = 0; return;
! 940: }
! 941: NEWDL(d,nv);
! 942: c = ONE;
! 943: for ( i = 0; i < nv; i++ ) {
! 944: for ( j = d2->d[i], k = d1->d[i]; k > 0; k--, j-- ) {
! 945: STOQ(j,jq); mulz(c,jq,&c1); c = c1;
! 946: }
! 947: d->d[i] = d2->d[i]-d1->d[i];
! 948: }
! 949: arf_mul(vl,C(m1),C(m2),&t);
! 950: NEWMP(m);
! 951: arf_mul(vl,(Obj)c,t,&C(m));
! 952: m->dl = d;
! 953: MKDP(nv,m,*pr);
! 954: }
! 955:
! 956: void weyl_actd(VL vl,DP p1,DP p2,DP *pr)
! 957: {
! 958: int n;
! 959: MP m1,m2;
! 960: DP d,r,s;
! 961:
! 962: if ( !p1 || !p2 ) *pr = 0;
! 963: else {
! 964: n = NV(p1);
! 965: r = 0;
! 966: for ( m1 = BDY(p1); m1; m1 = NEXT(m1) )
! 967: for ( m2 = BDY(p2); m2; m2 = NEXT(m2) ) {
! 968: actm(vl,n,m1,m2,&d);
! 969: addd(vl,r,d,&s); r = s;
! 970: }
! 971: *pr = r;
! 972: }
! 973: }
! 974:
! 975: /* monomial * polynomial */
! 976:
! 977: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
! 978: {
! 979: DP r,t,t1;
! 980: MP m;
! 981: DL d0;
! 982: int n,n2,l,i,j,tlen;
! 983: static MP *w,*psum;
! 984: static struct cdl *tab;
! 985: static int wlen;
! 986: static int rtlen;
! 987:
! 988: if ( !p )
! 989: *pr = 0;
! 990: else {
! 991: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 992: if ( l > wlen ) {
! 993: if ( w ) GCFREE(w);
! 994: w = (MP *)MALLOC(l*sizeof(MP));
! 995: wlen = l;
! 996: }
! 997: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 998: w[i] = m;
! 999:
! 1000: n = NV(p); n2 = n>>1;
! 1001: d0 = m0->dl;
! 1002: for ( i = 0, tlen = 1; i < n2; i++ )
! 1003: tlen *= d0->d[n2+i]+1;
! 1004: if ( tlen > rtlen ) {
! 1005: if ( tab ) GCFREE(tab);
! 1006: if ( psum ) GCFREE(psum);
! 1007: rtlen = tlen;
! 1008: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
! 1009: psum = (MP *)MALLOC(rtlen*sizeof(MP));
! 1010: }
! 1011: bzero(psum,tlen*sizeof(MP));
! 1012: for ( i = l-1; i >= 0; i-- ) {
! 1013: bzero(tab,tlen*sizeof(struct cdl));
! 1014: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
! 1015: for ( j = 0; j < tlen; j++ ) {
! 1016: if ( tab[j].c ) {
! 1017: NEWMP(m); m->dl = tab[j].d; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
! 1018: psum[j] = m;
! 1019: }
! 1020: }
! 1021: }
! 1022: for ( j = tlen-1, r = 0; j >= 0; j-- )
! 1023: if ( psum[j] ) {
! 1024: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
! 1025: }
! 1026: if ( r )
! 1027: r->sugar = p->sugar + m0->dl->td;
! 1028: *pr = r;
! 1029: }
! 1030: }
! 1031:
! 1032: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
! 1033: /* rtab : array of length (e0+1)*(e1+1)*... */
! 1034:
! 1035: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
! 1036: {
! 1037: Obj c,c0,c1;
! 1038: DL d,d0,d1,dt;
! 1039: int i,j,a,b,k,l,n2,s,min,curlen;
! 1040: struct cdl *p;
! 1041: static Z *ctab;
! 1042: static struct cdl *tab;
! 1043: static int tablen;
! 1044: static struct cdl *tmptab;
! 1045: static int tmptablen;
! 1046:
! 1047:
! 1048: if ( !m0 || !m1 ) {
! 1049: rtab[0].c = 0;
! 1050: rtab[0].d = 0;
! 1051: return;
! 1052: }
! 1053: c0 = C(m0); c1 = C(m1);
! 1054: arf_mul(vl,c0,c1,&c);
! 1055: d0 = m0->dl; d1 = m1->dl;
! 1056: n2 = n>>1;
! 1057: curlen = 1;
! 1058: NEWDL(d,n);
! 1059: if ( n & 1 )
! 1060: /* offset of h-degree */
! 1061: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
! 1062: else
! 1063: d->td = 0;
! 1064: rtab[0].c = c;
! 1065: rtab[0].d = d;
! 1066:
! 1067: if ( rtablen > tmptablen ) {
! 1068: if ( tmptab ) GCFREE(tmptab);
! 1069: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
! 1070: tmptablen = rtablen;
! 1071: }
! 1072: for ( i = 0; i < n2; i++ ) {
! 1073: a = d0->d[i]; b = d1->d[n2+i];
! 1074: k = d0->d[n2+i]; l = d1->d[i];
! 1075:
! 1076: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 1077: a += l;
! 1078: b += k;
! 1079: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
! 1080:
! 1081: if ( !k || !l ) {
! 1082: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
! 1083: if ( p->c ) {
! 1084: dt = p->d;
! 1085: dt->d[i] = a;
! 1086: dt->d[n2+i] = b;
! 1087: dt->td += s;
! 1088: }
! 1089: }
! 1090: curlen *= k+1;
! 1091: continue;
! 1092: }
! 1093: if ( k+1 > tablen ) {
! 1094: if ( tab ) GCFREE(tab);
! 1095: if ( ctab ) GCFREE(ctab);
! 1096: tablen = k+1;
! 1097: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
! 1098: ctab = (Z *)MALLOC(tablen*sizeof(Q));
! 1099: }
! 1100: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 1101: min = MIN(k,l);
! 1102: mkwc(k,l,ctab);
! 1103: bzero(tab,(k+1)*sizeof(struct cdl));
! 1104: if ( n & 1 )
! 1105: for ( j = 0; j <= min; j++ ) {
! 1106: NEWDL(d,n);
! 1107: d->d[i] = a-j; d->d[n2+i] = b-j;
! 1108: d->td = s;
! 1109: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
! 1110: tab[j].d = d;
! 1111: tab[j].c = (Obj)ctab[j];
! 1112: }
! 1113: else
! 1114: for ( j = 0; j <= min; j++ ) {
! 1115: NEWDL(d,n);
! 1116: d->d[i] = a-j; d->d[n2+i] = b-j;
! 1117: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
! 1118: tab[j].d = d;
! 1119: tab[j].c = (Obj)ctab[j];
! 1120: }
! 1121: bzero(ctab,(min+1)*sizeof(Q));
! 1122: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
! 1123: curlen *= k+1;
! 1124: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
! 1125: }
! 1126: }
! 1127:
! 1128: /* direct product of two cdl tables
! 1129: rt[] = [
! 1130: t[0]*t1[0],...,t[n-1]*t1[0],
! 1131: t[0]*t1[1],...,t[n-1]*t1[1],
! 1132: ...
! 1133: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
! 1134: ]
! 1135: */
! 1136:
! 1137: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
! 1138: {
! 1139: int i,j;
! 1140: struct cdl *p;
! 1141: Obj c;
! 1142: DL d;
! 1143:
! 1144: bzero(rt,n*n1*sizeof(struct cdl));
! 1145: for ( j = 0, p = rt; j < n1; j++ ) {
! 1146: c = (Obj)t1[j].c;
! 1147: d = t1[j].d;
! 1148: if ( !c )
! 1149: break;
! 1150: for ( i = 0; i < n; i++, p++ ) {
! 1151: if ( t[i].c ) {
! 1152: arf_mul(vl,(Obj)t[i].c,c,(Obj *)&p->c);
! 1153: adddl(nv,t[i].d,d,&p->d);
! 1154: }
! 1155: }
! 1156: }
! 1157: }
! 1158:
! 1159: void muldc(VL vl,DP p,Obj c,DP *pr)
! 1160: {
! 1161: MP m,mr=0,mr0;
! 1162:
! 1163: if ( !p || !c )
! 1164: *pr = 0;
! 1165: else if ( NUM(c) && UNIQ((Q)c) )
! 1166: *pr = p;
! 1167: else if ( NUM(c) && MUNIQ((Q)c) )
! 1168: chsgnd(p,pr);
! 1169: else {
! 1170: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 1171: NEXTMP(mr0,mr);
! 1172: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
! 1173: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
! 1174: else
! 1175: arf_mul(vl,C(m),c,&C(mr));
! 1176: mr->dl = m->dl;
! 1177: }
! 1178: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 1179: if ( *pr )
! 1180: (*pr)->sugar = p->sugar;
! 1181: }
! 1182: }
! 1183:
! 1184: void divdc(VL vl,DP p,Obj c,DP *pr)
! 1185: {
! 1186: Obj inv;
! 1187:
! 1188: arf_div(vl,(Obj)ONE,c,&inv);
! 1189: muld(vl,p,(DP)inv,pr);
! 1190: }
! 1191:
! 1192: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr)
! 1193: {
! 1194: MP m,mr=0,mr0;
! 1195: DL mdl;
! 1196: int i,n;
! 1197:
! 1198: if ( !p || !c ) {
! 1199: *pr = 0; return;
! 1200: }
! 1201: n = NV(p);
! 1202: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 1203: mdl = m->dl;
! 1204: for ( i = 0; i < n; i++ )
! 1205: if ( mdl->d[i] < dl->d[i] )
! 1206: break;
! 1207: if ( i < n )
! 1208: break;
! 1209: NEXTMP(mr0,mr);
! 1210: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
! 1211: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
! 1212: else
! 1213: arf_mul(vl,C(m),c,&C(mr));
! 1214: mr->dl = m->dl;
! 1215: }
! 1216: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 1217: if ( *pr )
! 1218: (*pr)->sugar = p->sugar;
! 1219: }
! 1220:
! 1221: void divsdc(VL vl,DP p,P c,DP *pr)
! 1222: {
! 1223: MP m,mr=0,mr0;
! 1224:
! 1225: if ( !c )
! 1226: error("disvsdc : division by 0");
! 1227: else if ( !p )
! 1228: *pr = 0;
! 1229: else if ( OID(c) > O_P )
! 1230: error("divsdc : invalid argument");
! 1231: else {
! 1232: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 1233: NEXTMP(mr0,mr); divsp(vl,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl;
! 1234: }
! 1235: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
! 1236: if ( *pr )
! 1237: (*pr)->sugar = p->sugar;
! 1238: }
! 1239: }
! 1240:
! 1241: void adddl(int n,DL d1,DL d2,DL *dr)
! 1242: {
! 1243: DL dt;
! 1244: int i;
! 1245:
! 1246: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
! 1247: dt->td = d1->td + d2->td;
! 1248: for ( i = 0; i < n; i++ )
! 1249: dt->d[i] = d1->d[i]+d2->d[i];
! 1250: }
! 1251:
! 1252: /* d1 += d2 */
! 1253:
! 1254: void adddl_destructive(int n,DL d1,DL d2)
! 1255: {
! 1256: int i;
! 1257:
! 1258: d1->td += d2->td;
! 1259: for ( i = 0; i < n; i++ )
! 1260: d1->d[i] += d2->d[i];
! 1261: }
! 1262:
! 1263: int compd(VL vl,DP p1,DP p2)
! 1264: {
! 1265: int n,t;
! 1266: MP m1,m2;
! 1267:
! 1268: if ( !p1 )
! 1269: return p2 ? -1 : 0;
! 1270: else if ( !p2 )
! 1271: return 1;
! 1272: else if ( NV(p1) != NV(p2) ) {
! 1273: error("compd : size mismatch");
! 1274: return 0; /* XXX */
! 1275: } else {
! 1276: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
! 1277: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
! 1278: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
! 1279: (t = arf_comp(vl,C(m1),C(m2)) ) )
! 1280: return t;
! 1281: if ( m1 )
! 1282: return 1;
! 1283: else if ( m2 )
! 1284: return -1;
! 1285: else
! 1286: return 0;
! 1287: }
! 1288: }
! 1289:
! 1290: int cmpdl_lex(int n,DL d1,DL d2)
! 1291: {
! 1292: int i;
! 1293:
! 1294: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
! 1295: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
! 1296: }
! 1297:
! 1298: int cmpdl_revlex(int n,DL d1,DL d2)
! 1299: {
! 1300: int i;
! 1301:
! 1302: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
! 1303: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
! 1304: }
! 1305:
! 1306: int cmpdl_gradlex(int n,DL d1,DL d2)
! 1307: {
! 1308: if ( d1->td > d2->td )
! 1309: return 1;
! 1310: else if ( d1->td < d2->td )
! 1311: return -1;
! 1312: else
! 1313: return cmpdl_lex(n,d1,d2);
! 1314: }
! 1315:
! 1316: int cmpdl_revgradlex(int n,DL d1,DL d2)
! 1317: {
! 1318: register int i,c;
! 1319: register int *p1,*p2;
! 1320:
! 1321: if ( d1->td > d2->td )
! 1322: return 1;
! 1323: else if ( d1->td < d2->td )
! 1324: return -1;
! 1325: else {
! 1326: i = n-1;
! 1327: p1 = d1->d+n-1;
! 1328: p2 = d2->d+n-1;
! 1329: while ( i >= 7 ) {
! 1330: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1331: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1332: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1333: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1334: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1335: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1336: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1337: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1338: i -= 8;
! 1339: }
! 1340: switch ( i ) {
! 1341: case 6:
! 1342: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1343: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1344: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1345: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1346: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1347: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1348: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1349: return 0;
! 1350: case 5:
! 1351: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1352: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1353: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1354: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1355: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1356: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1357: return 0;
! 1358: case 4:
! 1359: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1360: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1361: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1362: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1363: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1364: return 0;
! 1365: case 3:
! 1366: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1367: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1368: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1369: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1370: return 0;
! 1371: case 2:
! 1372: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1373: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1374: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1375: return 0;
! 1376: case 1:
! 1377: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1378: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1379: return 0;
! 1380: case 0:
! 1381: c = (*p1--) - (*p2--); if ( c ) goto LAST;
! 1382: return 0;
! 1383: default:
! 1384: return 0;
! 1385: }
! 1386: LAST:
! 1387: if ( c > 0 ) return -1;
! 1388: else return 1;
! 1389: }
! 1390: }
! 1391:
! 1392: int cmpdl_blex(int n,DL d1,DL d2)
! 1393: {
! 1394: int c;
! 1395:
! 1396: if ( (c = cmpdl_lex(n-1,d1,d2)) )
! 1397: return c;
! 1398: else {
! 1399: c = d1->d[n-1] - d2->d[n-1];
! 1400: return c > 0 ? 1 : c < 0 ? -1 : 0;
! 1401: }
! 1402: }
! 1403:
! 1404: int cmpdl_bgradlex(int n,DL d1,DL d2)
! 1405: {
! 1406: int e1,e2,c;
! 1407:
! 1408: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
! 1409: if ( e1 > e2 )
! 1410: return 1;
! 1411: else if ( e1 < e2 )
! 1412: return -1;
! 1413: else {
! 1414: c = cmpdl_lex(n-1,d1,d2);
! 1415: if ( c )
! 1416: return c;
! 1417: else
! 1418: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
! 1419: }
! 1420: }
! 1421:
! 1422: int cmpdl_brevgradlex(int n,DL d1,DL d2)
! 1423: {
! 1424: int e1,e2,c;
! 1425:
! 1426: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
! 1427: if ( e1 > e2 )
! 1428: return 1;
! 1429: else if ( e1 < e2 )
! 1430: return -1;
! 1431: else {
! 1432: c = cmpdl_revlex(n-1,d1,d2);
! 1433: if ( c )
! 1434: return c;
! 1435: else
! 1436: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
! 1437: }
! 1438: }
! 1439:
! 1440: int cmpdl_brevrev(int n,DL d1,DL d2)
! 1441: {
! 1442: int e1,e2,f1,f2,c,i;
! 1443:
! 1444: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
! 1445: e1 += d1->d[i]; e2 += d2->d[i];
! 1446: }
! 1447: f1 = d1->td - e1; f2 = d2->td - e2;
! 1448: if ( e1 > e2 )
! 1449: return 1;
! 1450: else if ( e1 < e2 )
! 1451: return -1;
! 1452: else {
! 1453: c = cmpdl_revlex(dp_nelim,d1,d2);
! 1454: if ( c )
! 1455: return c;
! 1456: else if ( f1 > f2 )
! 1457: return 1;
! 1458: else if ( f1 < f2 )
! 1459: return -1;
! 1460: else {
! 1461: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
! 1462: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
! 1463: }
! 1464: }
! 1465: }
! 1466:
! 1467: int cmpdl_bgradrev(int n,DL d1,DL d2)
! 1468: {
! 1469: int e1,e2,f1,f2,c,i;
! 1470:
! 1471: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
! 1472: e1 += d1->d[i]; e2 += d2->d[i];
! 1473: }
! 1474: f1 = d1->td - e1; f2 = d2->td - e2;
! 1475: if ( e1 > e2 )
! 1476: return 1;
! 1477: else if ( e1 < e2 )
! 1478: return -1;
! 1479: else {
! 1480: c = cmpdl_lex(dp_nelim,d1,d2);
! 1481: if ( c )
! 1482: return c;
! 1483: else if ( f1 > f2 )
! 1484: return 1;
! 1485: else if ( f1 < f2 )
! 1486: return -1;
! 1487: else {
! 1488: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
! 1489: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
! 1490: }
! 1491: }
! 1492: }
! 1493:
! 1494: int cmpdl_blexrev(int n,DL d1,DL d2)
! 1495: {
! 1496: int e1,e2,f1,f2,c,i;
! 1497:
! 1498: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
! 1499: e1 += d1->d[i]; e2 += d2->d[i];
! 1500: }
! 1501: f1 = d1->td - e1; f2 = d2->td - e2;
! 1502: c = cmpdl_lex(dp_nelim,d1,d2);
! 1503: if ( c )
! 1504: return c;
! 1505: else if ( f1 > f2 )
! 1506: return 1;
! 1507: else if ( f1 < f2 )
! 1508: return -1;
! 1509: else {
! 1510: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
! 1511: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
! 1512: }
! 1513: }
! 1514:
! 1515: int cmpdl_elim(int n,DL d1,DL d2)
! 1516: {
! 1517: int e1,e2,i;
! 1518:
! 1519: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
! 1520: e1 += d1->d[i]; e2 += d2->d[i];
! 1521: }
! 1522: if ( e1 > e2 )
! 1523: return 1;
! 1524: else if ( e1 < e2 )
! 1525: return -1;
! 1526: else
! 1527: return cmpdl_revgradlex(n,d1,d2);
! 1528: }
! 1529:
! 1530: int cmpdl_weyl_elim(int n,DL d1,DL d2)
! 1531: {
! 1532: int e1,e2,i;
! 1533:
! 1534: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
! 1535: e1 += d1->d[n-i]; e2 += d2->d[n-i];
! 1536: }
! 1537: if ( e1 > e2 )
! 1538: return 1;
! 1539: else if ( e1 < e2 )
! 1540: return -1;
! 1541: else if ( d1->td > d2->td )
! 1542: return 1;
! 1543: else if ( d1->td < d2->td )
! 1544: return -1;
! 1545: else return -cmpdl_revlex(n,d1,d2);
! 1546: }
! 1547:
! 1548: /*
! 1549: a special ordering
! 1550: 1. total order
! 1551: 2. (-w,w) for the first 2*m variables
! 1552: 3. DRL for the first 2*m variables
! 1553: */
! 1554:
! 1555: extern int *current_weyl_weight_vector;
! 1556:
! 1557: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
! 1558: {
! 1559: int e1,e2,m,i;
! 1560: int *p1,*p2;
! 1561:
! 1562: if ( d1->td > d2->td )
! 1563: return 1;
! 1564: else if ( d1->td < d2->td )
! 1565: return -1;
! 1566:
! 1567: m = n>>1;
! 1568: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
! 1569: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
! 1570: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
! 1571: }
! 1572: if ( e1 > e2 )
! 1573: return 1;
! 1574: else if ( e1 < e2 )
! 1575: return -1;
! 1576:
! 1577: e1 = d1->td - d1->d[n-1];
! 1578: e2 = d2->td - d2->d[n-1];
! 1579: if ( e1 > e2 )
! 1580: return 1;
! 1581: else if ( e1 < e2 )
! 1582: return -1;
! 1583:
! 1584: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
! 1585: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
! 1586: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
! 1587: }
! 1588:
! 1589: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
! 1590: {
! 1591: int i,t,m;
! 1592: int *p1,*p2;
! 1593:
! 1594: if ( d1->td > d2->td )
! 1595: return 1;
! 1596: else if ( d1->td < d2->td )
! 1597: return -1;
! 1598: else {
! 1599: m = n>>1;
! 1600: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
! 1601: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
! 1602: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
! 1603: }
! 1604: return 0;
! 1605: }
! 1606: }
! 1607:
! 1608: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
! 1609: {
! 1610: int e1,e2,m,i,t;
! 1611: int *p1,*p2;
! 1612:
! 1613: if ( d1->td > d2->td )
! 1614: return 1;
! 1615: else if ( d1->td < d2->td )
! 1616: return -1;
! 1617:
! 1618: m = n>>1;
! 1619: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
! 1620: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
! 1621: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
! 1622: }
! 1623: if ( e1 > e2 )
! 1624: return 1;
! 1625: else if ( e1 < e2 )
! 1626: return -1;
! 1627:
! 1628: e1 = d1->td - d1->d[n-1];
! 1629: e2 = d2->td - d2->d[n-1];
! 1630: if ( e1 > e2 )
! 1631: return 1;
! 1632: else if ( e1 < e2 )
! 1633: return -1;
! 1634:
! 1635: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
! 1636: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
! 1637: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
! 1638: }
! 1639: return 0;
! 1640: }
! 1641:
! 1642: int cmpdl_order_pair(int n,DL d1,DL d2)
! 1643: {
! 1644: int e1,e2,i,j,l;
! 1645: int *t1,*t2;
! 1646: int len,head;
! 1647: struct order_pair *pair;
! 1648:
! 1649: len = dp_current_spec->ord.block.length;
! 1650: if ( n != dp_current_spec->nv )
! 1651: error("cmpdl_order_pair : incompatible order specification");
! 1652: pair = dp_current_spec->ord.block.order_pair;
! 1653:
! 1654: head = 0;
! 1655: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
! 1656: l = pair[i].length;
! 1657: switch ( pair[i].order ) {
! 1658: case 0:
! 1659: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
! 1660: e1 += MUL_WEIGHT(t1[j],head+j);
! 1661: e2 += MUL_WEIGHT(t2[j],head+j);
! 1662: }
! 1663: if ( e1 > e2 )
! 1664: return 1;
! 1665: else if ( e1 < e2 )
! 1666: return -1;
! 1667: else {
! 1668: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
! 1669: if ( j >= 0 )
! 1670: return t1[j] < t2[j] ? 1 : -1;
! 1671: }
! 1672: break;
! 1673: case 1:
! 1674: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
! 1675: e1 += MUL_WEIGHT(t1[j],head+j);
! 1676: e2 += MUL_WEIGHT(t2[j],head+j);
! 1677: }
! 1678: if ( e1 > e2 )
! 1679: return 1;
! 1680: else if ( e1 < e2 )
! 1681: return -1;
! 1682: else {
! 1683: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
! 1684: if ( j < l )
! 1685: return t1[j] > t2[j] ? 1 : -1;
! 1686: }
! 1687: break;
! 1688: case 2:
! 1689: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
! 1690: if ( j < l )
! 1691: return t1[j] > t2[j] ? 1 : -1;
! 1692: break;
! 1693: default:
! 1694: error("cmpdl_order_pair : invalid order"); break;
! 1695: }
! 1696: t1 += l; t2 += l; head += l;
! 1697: }
! 1698: return 0;
! 1699: }
! 1700:
! 1701: int cmpdl_composite(int nv,DL d1,DL d2)
! 1702: {
! 1703: int n,i,j,k,start,s,len;
! 1704: int *dw;
! 1705: struct sparse_weight *sw;
! 1706: struct weight_or_block *worb;
! 1707: int *w,*t1,*t2;
! 1708:
! 1709: n = dp_current_spec->ord.composite.length;
! 1710: worb = dp_current_spec->ord.composite.w_or_b;
! 1711: w = dp_dl_work;
! 1712: for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
! 1713: w[i] = t1[i]-t2[i];
! 1714: for ( i = 0; i < n; i++, worb++ ) {
! 1715: len = worb->length;
! 1716: switch ( worb->type ) {
! 1717: case IS_DENSE_WEIGHT:
! 1718: dw = worb->body.dense_weight;
! 1719: for ( j = 0, s = 0; j < len; j++ )
! 1720: s += dw[j]*w[j];
! 1721: if ( s > 0 ) return 1;
! 1722: else if ( s < 0 ) return -1;
! 1723: break;
! 1724: case IS_SPARSE_WEIGHT:
! 1725: sw = worb->body.sparse_weight;
! 1726: for ( j = 0, s = 0; j < len; j++ )
! 1727: s += sw[j].value*w[sw[j].pos];
! 1728: if ( s > 0 ) return 1;
! 1729: else if ( s < 0 ) return -1;
! 1730: break;
! 1731: case IS_BLOCK:
! 1732: start = worb->body.block.start;
! 1733: switch ( worb->body.block.order ) {
! 1734: case 0:
! 1735: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
! 1736: s += MUL_WEIGHT(w[k],k);
! 1737: }
! 1738: if ( s > 0 ) return 1;
! 1739: else if ( s < 0 ) return -1;
! 1740: else {
! 1741: for ( j = k-1; j >= start && w[j] == 0; j-- );
! 1742: if ( j >= start )
! 1743: return w[j] < 0 ? 1 : -1;
! 1744: }
! 1745: break;
! 1746: case 1:
! 1747: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
! 1748: s += MUL_WEIGHT(w[k],k);
! 1749: }
! 1750: if ( s > 0 ) return 1;
! 1751: else if ( s < 0 ) return -1;
! 1752: else {
! 1753: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
! 1754: if ( j < len )
! 1755: return w[j] > 0 ? 1 : -1;
! 1756: }
! 1757: break;
! 1758: case 2:
! 1759: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
! 1760: if ( j < len )
! 1761: return w[j] > 0 ? 1 : -1;
! 1762: break;
! 1763: }
! 1764: break;
! 1765: }
! 1766: }
! 1767: return 0;
! 1768: }
! 1769:
! 1770: int cmpdl_matrix(int n,DL d1,DL d2)
! 1771: {
! 1772: int *v,*w,*t1,*t2;
! 1773: int s,i,j,len;
! 1774: int **matrix;
! 1775:
! 1776: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
! 1777: w[i] = t1[i]-t2[i];
! 1778: len = dp_current_spec->ord.matrix.row;
! 1779: matrix = dp_current_spec->ord.matrix.matrix;
! 1780: for ( j = 0; j < len; j++ ) {
! 1781: v = matrix[j];
! 1782: for ( i = 0, s = 0; i < n; i++ )
! 1783: s += v[i]*w[i];
! 1784: if ( s > 0 )
! 1785: return 1;
! 1786: else if ( s < 0 )
! 1787: return -1;
! 1788: }
! 1789: return 0;
! 1790: }
! 1791:
! 1792: int cmpdl_top_weight(int n,DL d1,DL d2)
! 1793: {
! 1794: int *w;
! 1795: Z **mat;
! 1796: Z *a;
! 1797: mpz_t sum;
! 1798: int len,i,sgn,tsgn,row,k;
! 1799: int *t1,*t2;
! 1800:
! 1801: w = (int *)ALLOCA(n*sizeof(int));
! 1802: len = current_top_weight_len+3;
! 1803: t1 = d1->d; t2 = d2->d;
! 1804: for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
! 1805: mpz_init_set_ui(sum,0);
! 1806: if ( OID(current_top_weight) == O_VECT ) {
! 1807: mat = (Z **)&BDY((VECT)current_top_weight);
! 1808: row = 1;
! 1809: } else {
! 1810: mat = (Z **)BDY((MAT)current_top_weight);
! 1811: row = ((MAT)current_top_weight)->row;
! 1812: }
! 1813: for ( k = 0; k < row; k++ ) {
! 1814: a = mat[k];
! 1815: for ( i = 0; i < n; i++ ) {
! 1816: if ( !a[i] || !w[i] ) continue;
! 1817: if ( w[i] > 0 )
! 1818: mpz_addmul_ui(sum,BDY(a[i]),(unsigned int)w[i]);
! 1819: else
! 1820: mpz_submul_ui(sum,BDY(a[i]),(unsigned int)(-w[i]));
! 1821: }
! 1822: sgn = mpz_sgn(sum);
! 1823: if ( sgn > 0 ) return 1;
! 1824: else if ( sgn < 0 ) return -1;
! 1825: }
! 1826: return (*cmpdl_tie_breaker)(n,d1,d2);
! 1827: }
! 1828:
! 1829: GeoBucket create_bucket()
! 1830: {
! 1831: GeoBucket g;
! 1832:
! 1833: g = CALLOC(1,sizeof(struct oGeoBucket));
! 1834: g->m = 32;
! 1835: return g;
! 1836: }
! 1837:
! 1838: int length(NODE d);
! 1839:
! 1840: void add_bucket(GeoBucket g,NODE d,int nv)
! 1841: {
! 1842: int l,k,m;
! 1843:
! 1844: l = length(d);
! 1845: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
! 1846: /* 2^(k-1) < l <= 2^k */
! 1847: d = symb_merge(g->body[k],d,nv);
! 1848: for ( ; length(d) > (1<<(k)); k++ ) {
! 1849: g->body[k] = 0;
! 1850: d = symb_merge(g->body[k+1],d,nv);
! 1851: }
! 1852: g->body[k] = d;
! 1853: g->m = MAX(g->m,k);
! 1854: }
! 1855:
! 1856: DL remove_head_bucket(GeoBucket g,int nv)
! 1857: {
! 1858: int j,i,c,m;
! 1859: DL d;
! 1860:
! 1861: j = -1;
! 1862: m = g->m;
! 1863: for ( i = 0; i <= m; i++ ) {
! 1864: if ( !g->body[i] )
! 1865: continue;
! 1866: if ( j < 0 ) j = i;
! 1867: else {
! 1868: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
! 1869: if ( c > 0 )
! 1870: j = i;
! 1871: else if ( c == 0 )
! 1872: g->body[i] = NEXT(g->body[i]);
! 1873: }
! 1874: }
! 1875: if ( j < 0 )
! 1876: return 0;
! 1877: else {
! 1878: d = g->body[j]->body;
! 1879: g->body[j] = NEXT(g->body[j]);
! 1880: return d;
! 1881: }
! 1882: }
! 1883:
! 1884: /* DPV functions */
! 1885:
! 1886: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
! 1887: {
! 1888: int i,len;
! 1889: DP *e;
! 1890:
! 1891: if ( !p1 || !p2 )
! 1892: error("adddv : invalid argument");
! 1893: else if ( p1->len != p2->len )
! 1894: error("adddv : size mismatch");
! 1895: else {
! 1896: len = p1->len;
! 1897: e = (DP *)MALLOC(p1->len*sizeof(DP));
! 1898: for ( i = 0; i < len; i++ )
! 1899: addd(vl,p1->body[i],p2->body[i],&e[i]);
! 1900: MKDPV(len,e,*pr);
! 1901: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 1902: }
! 1903: }
! 1904:
! 1905: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
! 1906: {
! 1907: int i,len;
! 1908: DP *e;
! 1909:
! 1910: if ( !p1 || !p2 )
! 1911: error("subdv : invalid argument");
! 1912: else if ( p1->len != p2->len )
! 1913: error("subdv : size mismatch");
! 1914: else {
! 1915: len = p1->len;
! 1916: e = (DP *)MALLOC(p1->len*sizeof(DP));
! 1917: for ( i = 0; i < len; i++ )
! 1918: subd(vl,p1->body[i],p2->body[i],&e[i]);
! 1919: MKDPV(len,e,*pr);
! 1920: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 1921: }
! 1922: }
! 1923:
! 1924: void chsgndv(DPV p1,DPV *pr)
! 1925: {
! 1926: int i,len;
! 1927: DP *e;
! 1928:
! 1929: if ( !p1 )
! 1930: error("subdv : invalid argument");
! 1931: else {
! 1932: len = p1->len;
! 1933: e = (DP *)MALLOC(p1->len*sizeof(DP));
! 1934: for ( i = 0; i < len; i++ )
! 1935: chsgnd(p1->body[i],&e[i]);
! 1936: MKDPV(len,e,*pr);
! 1937: (*pr)->sugar = p1->sugar;
! 1938: }
! 1939: }
! 1940:
! 1941: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
! 1942: {
! 1943: int i,len;
! 1944: DP *e;
! 1945:
! 1946: len = p2->len;
! 1947: e = (DP *)MALLOC(p2->len*sizeof(DP));
! 1948: if ( !p1 ) {
! 1949: MKDPV(len,e,*pr);
! 1950: (*pr)->sugar = 0;
! 1951: } else {
! 1952: for ( i = 0; i < len; i++ )
! 1953: muld(vl,p1,p2->body[i],&e[i]);
! 1954: MKDPV(len,e,*pr);
! 1955: (*pr)->sugar = p1->sugar + p2->sugar;
! 1956: }
! 1957: }
! 1958:
! 1959: int compdv(VL vl,DPV p1,DPV p2)
! 1960: {
! 1961: int i,t,len;
! 1962:
! 1963: if ( p1->len != p2->len ) {
! 1964: error("compdv : size mismatch");
! 1965: return 0; /* XXX */
! 1966: } else {
! 1967: len = p1->len;
! 1968: for ( i = 0; i < len; i++ )
! 1969: if ( (t = compd(vl,p1->body[i],p2->body[i])) )
! 1970: return t;
! 1971: return 0;
! 1972: }
! 1973: }
! 1974:
! 1975: int ni_next(int *a,int n)
! 1976: {
! 1977: int i,j,k,kj;
! 1978:
! 1979: /* find the first nonzero a[j] */
! 1980: for ( j = 0; j < n && a[j] == 0; j++ );
! 1981: /* find the first zero a[k] after a[j] */
! 1982: for ( k = j; k < n && a[k] == 1; k++ );
! 1983: if ( k == n ) return 0;
! 1984: /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
! 1985: /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
! 1986: kj = k-j-1;
! 1987: for ( i = 0; i < kj; i++ ) a[i] = 1;
! 1988: for ( ; i < k; i++ ) a[i] = 0;
! 1989: a[k] = 1;
! 1990: return 1;
! 1991: }
! 1992:
! 1993: int comp_nbm(NBM a,NBM b)
! 1994: {
! 1995: int d,i,ai,bi;
! 1996: unsigned int *ab,*bb;
! 1997:
! 1998: if ( a->d > b->d ) return 1;
! 1999: else if ( a->d < b->d ) return -1;
! 2000: else {
! 2001: d = a->d; ab = a->b; bb = b->b;
! 2002: #if 0
! 2003: w = (d+31)/32;
! 2004: for ( i = 0; i < w; i++ )
! 2005: if ( ab[i] > bb[i] ) return 1;
! 2006: else if ( ab[i] < bb[i] ) return -1;
! 2007: #else
! 2008: for ( i = 0; i < d; i++ ) {
! 2009: ai = NBM_GET(ab,i);
! 2010: bi = NBM_GET(bb,i);
! 2011: if ( ai > bi ) return 1;
! 2012: else if ( ai < bi ) return -1;
! 2013: }
! 2014: #endif
! 2015: return 0;
! 2016: }
! 2017: }
! 2018:
! 2019: NBM mul_nbm(NBM a,NBM b)
! 2020: {
! 2021: int ad,bd,d,i,j;
! 2022: unsigned int *ab,*bb,*mb;
! 2023: NBM m;
! 2024:
! 2025: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
! 2026: d = ad + bd;
! 2027: NEWNBM(m); NEWNBMBDY(m,d);
! 2028: m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
! 2029: j = 0;
! 2030: for ( i = 0; i < ad; i++, j++ )
! 2031: if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
! 2032: else NBM_CLR(mb,j);
! 2033: for ( i = 0; i < bd; i++, j++ )
! 2034: if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
! 2035: else NBM_CLR(mb,j);
! 2036: return m;
! 2037: }
! 2038:
! 2039: NBP nbmtonbp(NBM m)
! 2040: {
! 2041: NODE n;
! 2042: NBP u;
! 2043:
! 2044: MKNODE(n,m,0);
! 2045: MKNBP(u,n);
! 2046: return u;
! 2047: }
! 2048:
! 2049: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
! 2050:
! 2051: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
! 2052: {
! 2053: int i,d1;
! 2054: NBM t;
! 2055:
! 2056: if ( !a->d ) error("separate_nbm : invalid argument");
! 2057:
! 2058: if ( a0 ) {
! 2059: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
! 2060: *a0 = nbmtonbp(t);
! 2061: }
! 2062:
! 2063: if ( ah ) {
! 2064: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
! 2065: if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
! 2066: else NBM_CLR(t->b,0);
! 2067: *ah = nbmtonbp(t);
! 2068: }
! 2069:
! 2070: if ( ar ) {
! 2071: d1 = a->d-1;
! 2072: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
! 2073: for ( i = 0; i < d1; i++ ) {
! 2074: if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
! 2075: else NBM_CLR(t->b,i);
! 2076: }
! 2077: *ar = nbmtonbp(t);
! 2078: }
! 2079:
! 2080: return a->c;
! 2081: }
! 2082:
! 2083: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
! 2084:
! 2085: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
! 2086: {
! 2087: int i,d,d1;
! 2088: NBM t;
! 2089:
! 2090: if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
! 2091:
! 2092: if ( a0 ) {
! 2093: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
! 2094: *a0 = nbmtonbp(t);
! 2095: }
! 2096:
! 2097: d1 = a->d-1;
! 2098: if ( at ) {
! 2099: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
! 2100: if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
! 2101: else NBM_CLR(t->b,0);
! 2102: *at = nbmtonbp(t);
! 2103: }
! 2104:
! 2105: if ( ar ) {
! 2106: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
! 2107: for ( i = 0; i < d1; i++ ) {
! 2108: if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
! 2109: else NBM_CLR(t->b,i);
! 2110: }
! 2111: *ar = nbmtonbp(t);
! 2112: }
! 2113:
! 2114: return a->c;
! 2115: }
! 2116:
! 2117: NBP make_xky(int k)
! 2118: {
! 2119: int k1,i;
! 2120: NBM t;
! 2121:
! 2122: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
! 2123: k1 = k-1;
! 2124: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
! 2125: NBM_CLR(t->b,i);
! 2126: return nbmtonbp(t);
! 2127: }
! 2128:
! 2129: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
! 2130:
! 2131: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
! 2132: {
! 2133: int i,d1,k,k1;
! 2134: NBM t;
! 2135:
! 2136: if ( !a->d )
! 2137: error("separate_nbm : invalid argument");
! 2138: for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
! 2139: if ( i == a->d )
! 2140: error("separate_nbm : invalid argument");
! 2141: k1 = i;
! 2142: k = i+1;
! 2143:
! 2144: if ( a0 ) {
! 2145: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
! 2146: *a0 = nbmtonbp(t);
! 2147: }
! 2148:
! 2149: if ( ah ) {
! 2150: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
! 2151: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
! 2152: NBM_CLR(t->b,i);
! 2153: *ah = nbmtonbp(t);
! 2154: }
! 2155:
! 2156: if ( ar ) {
! 2157: d1 = a->d-k;
! 2158: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
! 2159: for ( i = 0; i < d1; i++ ) {
! 2160: if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
! 2161: else NBM_CLR(t->b,i);
! 2162: }
! 2163: *ar = nbmtonbp(t);
! 2164: }
! 2165:
! 2166: return a->c;
! 2167: }
! 2168:
! 2169: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
! 2170: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
! 2171: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
! 2172: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
! 2173:
! 2174: NBP shuffle_mul_nbm(NBM a,NBM b)
! 2175: {
! 2176: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
! 2177: P ac,bc,c;
! 2178:
! 2179: if ( !a->d || !b->d )
! 2180: u = nbmtonbp(mul_nbm(a,b));
! 2181: else {
! 2182: ac = separate_nbm(a,&a0,&ah,&ar);
! 2183: bc = separate_nbm(b,&b0,&bh,&br);
! 2184: mulp(CO,ac,bc,&c);
! 2185: shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
! 2186: shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
! 2187: addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
! 2188: }
! 2189: return u;
! 2190: }
! 2191:
! 2192: NBP harmonic_mul_nbm(NBM a,NBM b)
! 2193: {
! 2194: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
! 2195: P ac,bc,c;
! 2196:
! 2197: if ( !a->d || !b->d )
! 2198: u = nbmtonbp(mul_nbm(a,b));
! 2199: else {
! 2200: mulp(CO,a->c,b->c,&c);
! 2201: ac = separate_xky_nbm(a,&a0,&ah,&ar);
! 2202: bc = separate_xky_nbm(b,&b0,&bh,&br);
! 2203: mulp(CO,ac,bc,&c);
! 2204: harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
! 2205: harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
! 2206: abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
! 2207: harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
! 2208: addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
! 2209: }
! 2210: return u;
! 2211:
! 2212: }
! 2213:
! 2214: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
! 2215: {
! 2216: NODE b1,b2,br=0,br0;
! 2217: NBM m1,m2,m;
! 2218: P c;
! 2219:
! 2220: if ( !p1 )
! 2221: *rp = p2;
! 2222: else if ( !p2 )
! 2223: *rp = p1;
! 2224: else {
! 2225: for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
! 2226: m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
! 2227: switch ( comp_nbm(m1,m2) ) {
! 2228: case 0:
! 2229: addp(CO,m1->c,m2->c,&c);
! 2230: if ( c ) {
! 2231: NEXTNODE(br0,br);
! 2232: NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
! 2233: BDY(br) = (pointer)m;
! 2234: }
! 2235: b1 = NEXT(b1); b2 = NEXT(b2); break;
! 2236: case 1:
! 2237: NEXTNODE(br0,br); BDY(br) = BDY(b1);
! 2238: b1 = NEXT(b1); break;
! 2239: case -1:
! 2240: NEXTNODE(br0,br); BDY(br) = BDY(b2);
! 2241: b2 = NEXT(b2); break;
! 2242: }
! 2243: }
! 2244: if ( !br0 )
! 2245: if ( b1 )
! 2246: br0 = b1;
! 2247: else if ( b2 )
! 2248: br0 = b2;
! 2249: else {
! 2250: *rp = 0;
! 2251: return;
! 2252: }
! 2253: else if ( b1 )
! 2254: NEXT(br) = b1;
! 2255: else if ( b2 )
! 2256: NEXT(br) = b2;
! 2257: else
! 2258: NEXT(br) = 0;
! 2259: MKNBP(*rp,br0);
! 2260: }
! 2261: }
! 2262:
! 2263: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
! 2264: {
! 2265: NBP t;
! 2266:
! 2267: chsgnnbp(p2,&t);
! 2268: addnbp(vl,p1,t,rp);
! 2269: }
! 2270:
! 2271: void chsgnnbp(NBP p,NBP *rp)
! 2272: {
! 2273: NODE r0,r=0,b;
! 2274: NBM m,m1;
! 2275:
! 2276: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
! 2277: NEXTNODE(r0,r);
! 2278: m = (NBM)BDY(b);
! 2279: NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
! 2280: BDY(r) = m1;
! 2281: }
! 2282: if ( r0 ) NEXT(r) = 0;
! 2283: MKNBP(*rp,r0);
! 2284: }
! 2285:
! 2286: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
! 2287: {
! 2288: NODE b,n;
! 2289: NBP r,t,s;
! 2290: NBM m;
! 2291:
! 2292: if ( !p1 || !p2 ) {
! 2293: *rp = 0; return;
! 2294: }
! 2295: if ( OID(p1) != O_NBP ) {
! 2296: if ( !POLY(p1) )
! 2297: error("mulnbp : invalid argument");
! 2298: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
! 2299: MKNODE(n,m,0); MKNBP(p1,n);
! 2300: }
! 2301: if ( OID(p2) != O_NBP ) {
! 2302: if ( !POLY(p2) )
! 2303: error("mulnbp : invalid argument");
! 2304: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
! 2305: MKNODE(n,m,0); MKNBP(p2,n);
! 2306: }
! 2307: if ( length(BDY(p1)) < length(BDY(p2)) ) {
! 2308: for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
! 2309: mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
! 2310: addnbp(vl,r,t,&s); r = s;
! 2311: }
! 2312: *rp = r;
! 2313: } else {
! 2314: for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
! 2315: mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
! 2316: addnbp(vl,r,t,&s); r = s;
! 2317: }
! 2318: *rp = r;
! 2319: }
! 2320: }
! 2321:
! 2322: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
! 2323: {
! 2324: NODE b,r0,r=0;
! 2325:
! 2326: if ( !p ) *rp = 0;
! 2327: else {
! 2328: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
! 2329: NEXTNODE(r0,r);
! 2330: BDY(r) = mul_nbm(m,(NBM)BDY(b));
! 2331: }
! 2332: if ( r0 ) NEXT(r) = 0;
! 2333: MKNBP(*rp,r0);
! 2334: }
! 2335: }
! 2336:
! 2337: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
! 2338: {
! 2339: NODE b,r0,r=0;
! 2340:
! 2341: if ( !p ) *rp = 0;
! 2342: else {
! 2343: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
! 2344: NEXTNODE(r0,r);
! 2345: BDY(r) = mul_nbm((NBM)BDY(b),m);
! 2346: }
! 2347: if ( r0 ) NEXT(r) = 0;
! 2348: MKNBP(*rp,r0);
! 2349: }
! 2350: }
! 2351:
! 2352: void pwrnbp(VL vl,NBP a,Z q,NBP *c)
! 2353: {
! 2354: NBP a1,a2;
! 2355: Z q1,r1,two;
! 2356: NBM m;
! 2357: NODE r;
! 2358:
! 2359: if ( !q ) {
! 2360: NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
! 2361: MKNODE(r,m,0); MKNBP(*c,r);
! 2362: } else if ( !a )
! 2363: *c = 0;
! 2364: else if ( UNIQ(q) )
! 2365: *c = a;
! 2366: else {
! 2367: STOQ(2,two);
! 2368: divqrz(q,two,&q1,&r1);
! 2369: pwrnbp(vl,a,q1,&a1);
! 2370: mulnbp(vl,a1,a1,&a2);
! 2371: if ( r1 )
! 2372: mulnbp(vl,a2,a,c);
! 2373: else
! 2374: *c = a2;
! 2375: }
! 2376: }
! 2377:
! 2378: int compnbp(VL vl,NBP p1,NBP p2)
! 2379: {
! 2380: NODE n1,n2;
! 2381: NBM m1,m2;
! 2382: int t;
! 2383:
! 2384: if ( !p1 )
! 2385: return p2 ? -1 : 0;
! 2386: else if ( !p2 )
! 2387: return 1;
! 2388: else {
! 2389: for ( n1 = BDY(p1), n2 = BDY(p2);
! 2390: n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
! 2391: m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
! 2392: if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
! 2393: return t;
! 2394: }
! 2395: if ( n1 )
! 2396: return 1;
! 2397: else if ( n2 )
! 2398: return -1;
! 2399: else
! 2400: return 0;
! 2401: }
! 2402: }
! 2403:
! 2404: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
! 2405: {
! 2406: NODE b1,b2,n;
! 2407: NBP r,t,s;
! 2408: NBM m;
! 2409:
! 2410: if ( !p1 || !p2 ) {
! 2411: *rp = 0; return;
! 2412: }
! 2413: if ( OID(p1) != O_NBP ) {
! 2414: if ( !POLY(p1) )
! 2415: error("shuffle_mulnbp : invalid argument");
! 2416: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
! 2417: MKNODE(n,m,0); MKNBP(p1,n);
! 2418: }
! 2419: if ( OID(p2) != O_NBP ) {
! 2420: if ( !POLY(p2) )
! 2421: error("shuffle_mulnbp : invalid argument");
! 2422: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
! 2423: MKNODE(n,m,0); MKNBP(p2,n);
! 2424: }
! 2425: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
! 2426: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
! 2427: t = shuffle_mul_nbm(m,(NBM)BDY(b2));
! 2428: addnbp(vl,r,t,&s); r = s;
! 2429: }
! 2430: *rp = r;
! 2431: }
! 2432:
! 2433: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
! 2434: {
! 2435: NODE b1,b2,n;
! 2436: NBP r,t,s;
! 2437: NBM m;
! 2438:
! 2439: if ( !p1 || !p2 ) {
! 2440: *rp = 0; return;
! 2441: }
! 2442: if ( OID(p1) != O_NBP ) {
! 2443: if ( !POLY(p1) )
! 2444: error("harmonic_mulnbp : invalid argument");
! 2445: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
! 2446: MKNODE(n,m,0); MKNBP(p1,n);
! 2447: }
! 2448: if ( OID(p2) != O_NBP ) {
! 2449: if ( !POLY(p2) )
! 2450: error("harmonic_mulnbp : invalid argument");
! 2451: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
! 2452: MKNODE(n,m,0); MKNBP(p2,n);
! 2453: }
! 2454: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
! 2455: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
! 2456: t = harmonic_mul_nbm(m,(NBM)BDY(b2));
! 2457: addnbp(vl,r,t,&s); r = s;
! 2458: }
! 2459: *rp = r;
! 2460: }
! 2461:
! 2462: #if 0
! 2463: NBP shuffle_mul_nbm(NBM a,NBM b)
! 2464: {
! 2465: int ad,bd,d,i,ai,bi,bit,s;
! 2466: int *ab,*bb,*wmb,*w;
! 2467: NBM wm,tm;
! 2468: P c,c1;
! 2469: NODE r,t,t1,p;
! 2470: NBP u;
! 2471:
! 2472: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
! 2473: d = ad + bd;
! 2474: w = (int *)ALLOCA(d*sizeof(int));
! 2475: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2476: for ( i = 0; i < ad; i++ ) w[i] = 1;
! 2477: for ( ; i < d; i++ ) w[i] = 0;
! 2478: mulp(CO,a->c,b->c,&c);
! 2479: r = 0;
! 2480: do {
! 2481: wm->d = d; wm->c = c;
! 2482: ai = 0; bi = 0;
! 2483: for ( i = 0; i < d; i++ ) {
! 2484: if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
! 2485: else { bit = NBM_GET(bb,bi); bi++; }
! 2486: if ( bit ) NBM_SET(wmb,i);
! 2487: else NBM_CLR(wmb,i);
! 2488: }
! 2489: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
! 2490: tm = (NBM)BDY(t);
! 2491: s = comp_nbm(tm,wm);
! 2492: if ( s < 0 ) {
! 2493: /* insert */
! 2494: MKNODE(t1,wm,t);
! 2495: if ( !p ) r = t1;
! 2496: else NEXT(p) = t1;
! 2497: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2498: break;
! 2499: } else if ( s == 0 ) {
! 2500: /* add coefs */
! 2501: addp(CO,tm->c,c,&c1);
! 2502: if ( c1 ) tm->c = c1;
! 2503: else NEXT(p) = NEXT(t);
! 2504: break;
! 2505: }
! 2506: }
! 2507: if ( !t ) {
! 2508: /* append */
! 2509: MKNODE(t1,wm,t);
! 2510: if ( !p ) r = t1;
! 2511: else NEXT(p) = t1;
! 2512: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2513: }
! 2514: } while ( ni_next(w,d) );
! 2515: MKNBP(u,r);
! 2516: return u;
! 2517: }
! 2518:
! 2519: int nbmtoxky(NBM a,int *b)
! 2520: {
! 2521: int d,i,j,k;
! 2522: int *p;
! 2523:
! 2524: d = a->d; p = a->b;
! 2525: for ( i = j = 0, k = 1; i < d; i++ ) {
! 2526: if ( !NBM_GET(p,i) ) {
! 2527: b[j++] = k;
! 2528: k = 1;
! 2529: } else k++;
! 2530: }
! 2531: return j;
! 2532: }
! 2533:
! 2534: NBP harmonic_mul_nbm(NBM a,NBM b)
! 2535: {
! 2536: int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
! 2537: int i,j,k,ia,ib,s;
! 2538: int *wa,*wb,*w,*wab,*wa1,*wmb;
! 2539: P c,c1;
! 2540: NBM wm,tm;
! 2541: NODE r,t1,t,p;
! 2542: NBP u;
! 2543:
! 2544: da = a->d; db = b->d; d = da+db;
! 2545: wa = (int *)ALLOCA(da*sizeof(int));
! 2546: wb = (int *)ALLOCA(db*sizeof(int));
! 2547: la = nbmtoxky(a,wa);
! 2548: lb = nbmtoxky(b,wb);
! 2549: mulp(CO,a->c,b->c,&c);
! 2550: /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
! 2551: /* lmax : total length */
! 2552: lmax = la+lb;
! 2553: lmin = la>lb?la:lb;
! 2554: w = (int *)ALLOCA(lmax*sizeof(int));
! 2555: /* position of a+b */
! 2556: wab = (int *)ALLOCA(lmax*sizeof(int));
! 2557: /* position of a */
! 2558: wa1 = (int *)ALLOCA(lmax*sizeof(int));
! 2559: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2560: for ( l = lmin, r = 0; l <= lmax; l++ ) {
! 2561: lab = lmax - l;
! 2562: la1 = la - lab;
! 2563: lb1 = lb - lab;
! 2564: lab1 = l-lab;
! 2565: /* partion l into three parts: a, b, a+b */
! 2566: /* initialize wab */
! 2567: for ( i = 0; i < lab; i++ ) wab[i] = 1;
! 2568: for ( ; i < l; i++ ) wab[i] = 0;
! 2569: do {
! 2570: /* initialize wa1 */
! 2571: for ( i = 0; i < la1; i++ ) wa1[i] = 1;
! 2572: for ( ; i < lab1; i++ ) wa1[i] = 0;
! 2573: do {
! 2574: ia = 0; ib = 0;
! 2575: for ( i = j = 0; i < l; i++ )
! 2576: if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
! 2577: else if ( wa1[j++] ) w[i] = wa[ia++];
! 2578: else w[i] = wb[ib++];
! 2579: for ( i = j = 0; i < l; i++ ) {
! 2580: for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
! 2581: NBM_CLR(wmb,j); j++;
! 2582: }
! 2583: wm->d = j; wm->c = c;
! 2584: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
! 2585: tm = (NBM)BDY(t);
! 2586: s = comp_nbm(tm,wm);
! 2587: if ( s < 0 ) {
! 2588: /* insert */
! 2589: MKNODE(t1,wm,t);
! 2590: if ( !p ) r = t1;
! 2591: else NEXT(p) = t1;
! 2592: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2593: break;
! 2594: } else if ( s == 0 ) {
! 2595: /* add coefs */
! 2596: addp(CO,tm->c,c,&c1);
! 2597: if ( c1 ) tm->c = c1;
! 2598: else NEXT(p) = NEXT(t);
! 2599: break;
! 2600: }
! 2601: }
! 2602: if ( !t ) {
! 2603: /* append */
! 2604: MKNODE(t1,wm,t);
! 2605: if ( !p ) r = t1;
! 2606: else NEXT(p) = t1;
! 2607: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2608: }
! 2609: } while ( ni_next(wa1,lab1) );
! 2610: } while ( ni_next(wab,l) );
! 2611: }
! 2612: MKNBP(u,r);
! 2613: return u;
! 2614: }
! 2615: #endif
! 2616:
! 2617: /* DPM functions */
! 2618:
! 2619: int compdmm(int n,DMM m1,DMM m2)
! 2620: {
! 2621: int t;
! 2622:
! 2623: if ( dpm_ispot ) {
! 2624: if ( m1->pos < m2->pos ) return 1;
! 2625: else if ( m1->pos > m2->pos ) return -1;
! 2626: else return (*cmpdl)(n,m1->dl,m2->dl);
! 2627: } else {
! 2628: t = (*cmpdl)(n,m1->dl,m2->dl);
! 2629: if ( t ) return t;
! 2630: else if ( m1->pos < m2->pos ) return 1;
! 2631: else if ( m1->pos > m2->pos ) return -1;
! 2632: else return 0;
! 2633: }
! 2634: }
! 2635:
! 2636: void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
! 2637: {
! 2638: int n;
! 2639: DMM m1,m2,mr=0,mr0;
! 2640: Obj t;
! 2641: DL d;
! 2642:
! 2643: if ( !p1 )
! 2644: *pr = p2;
! 2645: else if ( !p2 )
! 2646: *pr = p1;
! 2647: else {
! 2648: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
! 2649: switch ( compdmm(n,m1,m2) ) {
! 2650: case 0:
! 2651: arf_add(vl,C(m1),C(m2),&t);
! 2652: if ( t ) {
! 2653: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = t;
! 2654: }
! 2655: m1 = NEXT(m1); m2 = NEXT(m2); break;
! 2656: case 1:
! 2657: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = C(m1);
! 2658: m1 = NEXT(m1); break;
! 2659: case -1:
! 2660: NEXTDMM(mr0,mr); mr->pos = m2->pos; mr->dl = m2->dl; C(mr) = C(m2);
! 2661: m2 = NEXT(m2); break;
! 2662: }
! 2663: if ( !mr0 )
! 2664: if ( m1 )
! 2665: mr0 = m1;
! 2666: else if ( m2 )
! 2667: mr0 = m2;
! 2668: else {
! 2669: *pr = 0;
! 2670: return;
! 2671: }
! 2672: else if ( m1 )
! 2673: NEXT(mr) = m1;
! 2674: else if ( m2 )
! 2675: NEXT(mr) = m2;
! 2676: else
! 2677: NEXT(mr) = 0;
! 2678: MKDPM(NV(p1),mr0,*pr);
! 2679: if ( *pr )
! 2680: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 2681: }
! 2682: }
! 2683:
! 2684: void subdpm(VL vl,DPM p1,DPM p2,DPM *pr)
! 2685: {
! 2686: DPM t;
! 2687:
! 2688: if ( !p2 )
! 2689: *pr = p1;
! 2690: else {
! 2691: chsgndpm(p2,&t); adddpm(vl,p1,t,pr);
! 2692: }
! 2693: }
! 2694:
! 2695: void chsgndpm(DPM p,DPM *pr)
! 2696: {
! 2697: DMM m,mr=0,mr0;
! 2698: Obj r;
! 2699:
! 2700: if ( !p )
! 2701: *pr = 0;
! 2702: else {
! 2703: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2704: NEXTDMM(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->pos = m->pos; mr->dl = m->dl;
! 2705: }
! 2706: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
! 2707: if ( *pr )
! 2708: (*pr)->sugar = p->sugar;
! 2709: }
! 2710: }
! 2711:
! 2712: void mulcdpm(VL vl,Obj c,DPM p,DPM *pr)
! 2713: {
! 2714: DMM m,mr=0,mr0;
! 2715:
! 2716: if ( !p || !c )
! 2717: *pr = 0;
! 2718: else if ( NUM(c) && UNIQ((Q)c) )
! 2719: *pr = p;
! 2720: else if ( NUM(c) && MUNIQ((Q)c) )
! 2721: chsgndpm(p,pr);
! 2722: else {
! 2723: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2724: NEXTDMM(mr0,mr);
! 2725: arf_mul(vl,C(m),c,&C(mr));
! 2726: mr->pos = m->pos;
! 2727: mr->dl = m->dl;
! 2728: }
! 2729: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
! 2730: if ( *pr )
! 2731: (*pr)->sugar = p->sugar;
! 2732: }
! 2733: }
! 2734:
! 2735: void comm_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
! 2736: {
! 2737: DMM m,mr=0,mr0;
! 2738: DL d;
! 2739: Obj c;
! 2740: int n;
! 2741:
! 2742: if ( !p )
! 2743: *pr = 0;
! 2744: else {
! 2745: n = NV(p);
! 2746: d = m0->dl;
! 2747: c = C(m0);
! 2748: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2749: NEXTDMM(mr0,mr);
! 2750: arf_mul(vl,C(m),c,&C(mr));
! 2751: mr->pos = m->pos;
! 2752: adddl(n,m->dl,d,&mr->dl);
! 2753: }
! 2754: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
! 2755: if ( *pr )
! 2756: (*pr)->sugar = p->sugar;
! 2757: }
! 2758: }
! 2759:
! 2760: void weyl_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
! 2761: {
! 2762: DPM r,t,t1;
! 2763: DMM m;
! 2764: DL d0;
! 2765: int n,n2,l,i,j,tlen;
! 2766: struct oMP mp;
! 2767: static DMM *w,*psum;
! 2768: static struct cdl *tab;
! 2769: static int wlen;
! 2770: static int rtlen;
! 2771:
! 2772: if ( !p )
! 2773: *pr = 0;
! 2774: else {
! 2775: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 2776: if ( l > wlen ) {
! 2777: if ( w ) GCFREE(w);
! 2778: w = (DMM *)MALLOC(l*sizeof(DMM));
! 2779: wlen = l;
! 2780: }
! 2781: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 2782: w[i] = m;
! 2783:
! 2784: n = NV(p); n2 = n>>1;
! 2785: d0 = m0->dl;
! 2786: for ( i = 0, tlen = 1; i < n2; i++ )
! 2787: tlen *= d0->d[n2+i]+1;
! 2788: if ( tlen > rtlen ) {
! 2789: if ( tab ) GCFREE(tab);
! 2790: if ( psum ) GCFREE(psum);
! 2791: rtlen = tlen;
! 2792: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
! 2793: psum = (DMM *)MALLOC(rtlen*sizeof(DMM));
! 2794: }
! 2795: bzero(psum,tlen*sizeof(DMM));
! 2796: for ( i = l-1; i >= 0; i-- ) {
! 2797: bzero(tab,tlen*sizeof(struct cdl));
! 2798: mp.dl = w[i]->dl; mp.c = C(w[i]); mp.next = 0;
! 2799: weyl_mulmm(vl,m0,&mp,n,tab,tlen);
! 2800: for ( j = 0; j < tlen; j++ ) {
! 2801: if ( tab[j].c ) {
! 2802: NEWDMM(m); m->dl = tab[j].d; m->pos = w[i]->pos; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
! 2803: psum[j] = m;
! 2804: }
! 2805: }
! 2806: }
! 2807: for ( j = tlen-1, r = 0; j >= 0; j-- )
! 2808: if ( psum[j] ) {
! 2809: MKDPM(n,psum[j],t); adddpm(vl,r,t,&t1); r = t1;
! 2810: }
! 2811: if ( r )
! 2812: r->sugar = p->sugar + m0->dl->td;
! 2813: *pr = r;
! 2814: }
! 2815: }
! 2816:
! 2817: void mulobjdpm(VL vl,Obj p1,DPM p2,DPM *pr)
! 2818: {
! 2819: MP m;
! 2820: DPM s,t,u;
! 2821:
! 2822: if ( !p1 || !p2 )
! 2823: *pr = 0;
! 2824: else if ( OID(p1) != O_DP )
! 2825: mulcdpm(vl,p1,p2,pr);
! 2826: else {
! 2827: s = 0;
! 2828: for ( m = BDY((DP)p1); m; m = NEXT(m) ) {
! 2829: if ( do_weyl )
! 2830: weyl_mulmpdpm(vl,m,p2,&t);
! 2831: else
! 2832: comm_mulmpdpm(vl,m,p2,&t);
! 2833: adddpm(vl,s,t,&u); s = u;
! 2834: }
! 2835: *pr = s;
! 2836: }
! 2837: }
! 2838:
! 2839: int compdpm(VL vl,DPM p1,DPM p2)
! 2840: {
! 2841: int n,t;
! 2842: DMM m1,m2;
! 2843:
! 2844: if ( !p1 )
! 2845: return p2 ? -1 : 0;
! 2846: else if ( !p2 )
! 2847: return 1;
! 2848: else if ( NV(p1) != NV(p2) ) {
! 2849: error("compdpm : size mismatch");
! 2850: return 0; /* XXX */
! 2851: } else {
! 2852: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
! 2853: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
! 2854: if ( (t = compdmm(n,m1,m2)) ||
! 2855: (t = arf_comp(vl,C(m1),C(m2)) ) )
! 2856: return t;
! 2857: if ( m1 )
! 2858: return 1;
! 2859: else if ( m2 )
! 2860: return -1;
! 2861: else
! 2862: return 0;
! 2863: }
! 2864: }
! 2865:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dist.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine