Annotation of OpenXM_contrib2/asir2000/engine/dist.c, Revision 1.53
1.8 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.9 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.8 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.53 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.52 2017/08/31 02:36:21 noro Exp $
1.8 noro 49: */
1.1 noro 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
1.12 noro 62: #define ORD_WEYL_ELIM 10
1.13 noro 63: #define ORD_HOMO_WW_DRL 11
1.21 noro 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();
1.43 noro 68: int cmpdl_top_weight();
1.1 noro 69:
70: int (*cmpdl)()=cmpdl_revgradlex;
1.43 noro 71: int (*cmpdl_tie_breaker)();
1.1 noro 72: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
73:
1.50 noro 74: Obj current_top_weight;
75: int current_top_weight_len;
76:
1.2 noro 77: int do_weyl;
78:
1.1 noro 79: int dp_nelim,dp_fcoeffs;
1.27 noro 80: struct order_spec *dp_current_spec;
1.31 noro 81: struct modorder_spec *dp_current_modspec;
1.1 noro 82: int *dp_dl_work;
83:
1.24 noro 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);
1.52 noro 87: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr);
1.47 noro 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);
1.29 noro 90:
91: void order_init()
92: {
1.53 ! noro 93: struct order_spec *spec;
1.29 noro 94:
1.53 ! noro 95: create_order_spec(0,0,&spec);
! 96: initd(spec);
! 97: create_modorder_spec(0,0,&dp_current_modspec);
1.29 noro 98: }
1.24 noro 99:
1.52 noro 100: int has_sfcoef_p(Obj f);
1.47 noro 101:
1.22 noro 102: int has_sfcoef(DP f)
1.1 noro 103: {
1.53 ! noro 104: MP t;
1.1 noro 105:
1.53 ! noro 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;
1.1 noro 112: }
113:
1.52 noro 114: int has_sfcoef_p(Obj f)
1.1 noro 115: {
1.53 ! noro 116: DCP dc;
1.1 noro 117:
1.53 ! noro 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
1.52 noro 128: return 0;
1.1 noro 129: }
130:
1.50 noro 131: extern Obj nd_top_weight;
132:
133: void reset_top_weight()
134: {
1.53 ! noro 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;
1.50 noro 140: }
1.43 noro 141:
1.19 noro 142: void initd(struct order_spec *spec)
1.1 noro 143: {
1.53 ! noro 144: int len,i,k,row;
1.48 noro 145: Q **mat;
1.43 noro 146:
1.53 ! noro 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: }
1.48 noro 202: for ( k = 0, len = 0; k < row; k++ )
1.53 ! noro 203: for ( i = 0; i < spec->nv; i++ )
! 204: if ( mat[k][i] )
! 205: len = MAX(PL(NM(mat[k][i])),len);
! 206: current_top_weight_len = len;
! 207: }
! 208: dp_current_spec = spec;
1.1 noro 209: }
210:
1.52 noro 211: int dpm_ispot;
212:
213: /* type=0 => TOP, type=1 => POT */
214: void initdpm(struct order_spec *spec,int type)
215: {
1.53 ! noro 216: int len,i,k,row;
1.52 noro 217: Q **mat;
218:
219: initd(spec);
220: dpm_ispot = type;
221: }
222:
1.19 noro 223: void ptod(VL vl,VL dvl,P p,DP *pr)
1.1 noro 224: {
1.53 ! noro 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;
1.1 noro 234:
1.53 ! noro 235: if ( !p )
! 236: *pr = 0;
1.52 noro 237: else if ( OID(p) > O_P )
238: error("ptod : only polynomials can be converted.");
1.53 ! noro 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: }
1.17 noro 275: #if 0
1.53 ! noro 276: if ( !dp_fcoeffs && has_sfcoef(*pr) )
! 277: dp_fcoeffs = N_GFS;
1.17 noro 278: #endif
1.1 noro 279: }
280:
1.52 noro 281: void dtop(VL vl,VL dvl,DP p,Obj *pr)
1.1 noro 282: {
1.53 ! noro 283: int n,i,j,k;
! 284: DL d;
! 285: MP m;
! 286: MP *a;
! 287: P r;
1.52 noro 288: Obj t,w,s,u;
1.53 ! noro 289: Q q;
! 290: VL tvl;
1.1 noro 291:
1.53 ! noro 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: }
1.1 noro 315: }
316:
1.19 noro 317: void nodetod(NODE node,DP *dp)
1.1 noro 318: {
1.53 ! noro 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;
1.1 noro 341: }
342:
1.52 noro 343: void nodetodpm(NODE node,Obj pos,DPM *dp)
344: {
1.53 ! noro 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;
1.52 noro 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:
1.19 noro 386: int sugard(MP m)
1.1 noro 387: {
1.53 ! noro 388: int s;
1.1 noro 389:
1.53 ! noro 390: for ( s = 0; m; m = NEXT(m) )
! 391: s = MAX(s,m->dl->td);
! 392: return s;
1.1 noro 393: }
394:
1.19 noro 395: void addd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 396: {
1.53 ! noro 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: }
1.1 noro 451: }
452:
453: /* for F4 symbolic reduction */
454:
1.19 noro 455: void symb_addd(DP p1,DP p2,DP *pr)
1.1 noro 456: {
1.53 ! noro 457: int n;
! 458: MP m1,m2,mr=0,mr0;
1.1 noro 459:
1.53 ! noro 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: }
1.3 noro 498: }
499:
500: /*
501: * destructive merge of two list
502: *
503: * p1, p2 : list of DL
504: * return : a merged list
505: */
506:
1.19 noro 507: NODE symb_merge(NODE m1,NODE m2,int n)
1.3 noro 508: {
1.53 ! noro 509: NODE top=0,prev,cur,m=0,t;
! 510: DL d1,d2;
1.3 noro 511:
1.53 ! noro 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);
1.26 noro 533: #if 1
1.53 ! noro 534: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
1.26 noro 535: #else
1.53 ! noro 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 ) {
1.25 noro 551: #endif
1.53 ! noro 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: }
1.23 noro 570: }
571:
572: void _adddl(int n,DL d1,DL d2,DL d3)
573: {
1.53 ! noro 574: int i;
1.23 noro 575:
1.53 ! noro 576: d3->td = d1->td+d2->td;
! 577: for ( i = 0; i < n; i++ )
! 578: d3->d[i] = d1->d[i]+d2->d[i];
1.23 noro 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: {
1.53 ! noro 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: }
1.18 noro 640: }
641:
1.19 noro 642: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
1.18 noro 643: {
1.53 ! noro 644: DLBUCKET top,prev,cur,m,t;
1.18 noro 645:
1.53 ! noro 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: }
1.1 noro 679: }
680:
1.19 noro 681: void subd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 682: {
1.53 ! noro 683: DP t;
1.1 noro 684:
1.53 ! noro 685: if ( !p2 )
! 686: *pr = p1;
! 687: else {
! 688: chsgnd(p2,&t); addd(vl,p1,t,pr);
! 689: }
1.1 noro 690: }
691:
1.19 noro 692: void chsgnd(DP p,DP *pr)
1.1 noro 693: {
1.53 ! noro 694: MP m,mr=0,mr0;
! 695: Obj r;
1.1 noro 696:
1.53 ! noro 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: }
1.1 noro 709: }
710:
1.19 noro 711: void muld(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 712: {
1.53 ! noro 713: if ( ! do_weyl )
! 714: comm_muld(vl,p1,p2,pr);
! 715: else
! 716: weyl_muld(vl,p1,p2,pr);
1.2 noro 717: }
718:
1.19 noro 719: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 720: {
1.53 ! noro 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: }
1.1 noro 753: }
754:
1.24 noro 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: {
1.53 ! noro 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: }
1.24 noro 791: }
792:
793: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
794: {
1.53 ! noro 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: }
1.24 noro 832: }
833:
1.19 noro 834: void muldm(VL vl,DP p,MP m0,DP *pr)
1.1 noro 835: {
1.53 ! noro 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: }
1.2 noro 857: }
858:
1.24 noro 859: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
860: {
1.53 ! noro 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: }
1.24 noro 894: }
895:
1.19 noro 896: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 897: {
1.53 ! noro 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: }
1.2 noro 925: }
926:
1.51 noro 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: Q jq,c,c1;
932: MP m;
1.52 noro 933: Obj t;
1.51 noro 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); mulq(c,jq,&c1); c = c1;
946: }
947: d->d[i] = d2->d[i]-d1->d[i];
948: }
1.52 noro 949: arf_mul(vl,C(m1),C(m2),&t);
1.51 noro 950: NEWMP(m);
1.52 noro 951: arf_mul(vl,(Obj)c,t,&C(m));
1.51 noro 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:
1.10 noro 975: /* monomial * polynomial */
976:
1.19 noro 977: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1.2 noro 978: {
1.53 ! noro 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: }
1.2 noro 1030: }
1031:
1.10 noro 1032: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
1033: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2 noro 1034:
1.19 noro 1035: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.2 noro 1036: {
1.52 noro 1037: Obj c,c0,c1;
1.53 ! noro 1038: DL d,d0,d1,dt;
! 1039: int i,j,a,b,k,l,n2,s,min,curlen;
! 1040: struct cdl *p;
! 1041: static Q *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 = (Q *)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: }
1.10 noro 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:
1.19 noro 1137: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.10 noro 1138: {
1.53 ! noro 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: }
1.1 noro 1157: }
1158:
1.52 noro 1159: void muldc(VL vl,DP p,Obj c,DP *pr)
1.1 noro 1160: {
1.53 ! noro 1161: MP m,mr=0,mr0;
1.1 noro 1162:
1.53 ! noro 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: }
1.24 noro 1182: }
1183:
1.52 noro 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)
1.24 noro 1193: {
1.53 ! noro 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;
1.1 noro 1219: }
1220:
1.19 noro 1221: void divsdc(VL vl,DP p,P c,DP *pr)
1.1 noro 1222: {
1.53 ! noro 1223: MP m,mr=0,mr0;
1.1 noro 1224:
1.53 ! noro 1225: if ( !c )
! 1226: error("disvsdc : division by 0");
! 1227: else if ( !p )
! 1228: *pr = 0;
1.52 noro 1229: else if ( OID(p) > O_P )
1.53 ! noro 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: }
1.1 noro 1239: }
1240:
1.19 noro 1241: void adddl(int n,DL d1,DL d2,DL *dr)
1.1 noro 1242: {
1.53 ! noro 1243: DL dt;
! 1244: int i;
1.1 noro 1245:
1.53 ! noro 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];
1.11 noro 1250: }
1251:
1252: /* d1 += d2 */
1253:
1.19 noro 1254: void adddl_destructive(int n,DL d1,DL d2)
1.11 noro 1255: {
1.53 ! noro 1256: int i;
1.11 noro 1257:
1.53 ! noro 1258: d1->td += d2->td;
! 1259: for ( i = 0; i < n; i++ )
! 1260: d1->d[i] += d2->d[i];
1.1 noro 1261: }
1262:
1.19 noro 1263: int compd(VL vl,DP p1,DP p2)
1.1 noro 1264: {
1.53 ! noro 1265: int n,t;
! 1266: MP m1,m2;
1.1 noro 1267:
1.53 ! noro 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: }
1.1 noro 1288: }
1289:
1.19 noro 1290: int cmpdl_lex(int n,DL d1,DL d2)
1.1 noro 1291: {
1.53 ! noro 1292: int i;
1.1 noro 1293:
1.53 ! noro 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);
1.1 noro 1296: }
1297:
1.19 noro 1298: int cmpdl_revlex(int n,DL d1,DL d2)
1.1 noro 1299: {
1.53 ! noro 1300: int i;
1.1 noro 1301:
1.53 ! noro 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);
1.1 noro 1304: }
1305:
1.19 noro 1306: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1 noro 1307: {
1.53 ! noro 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);
1.1 noro 1314: }
1315:
1.19 noro 1316: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1 noro 1317: {
1.53 ! noro 1318: register int i,c;
! 1319: register int *p1,*p2;
1.7 noro 1320:
1.53 ! noro 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: }
1.25 noro 1386: LAST:
1.53 ! noro 1387: if ( c > 0 ) return -1;
! 1388: else return 1;
! 1389: }
1.1 noro 1390: }
1391:
1.19 noro 1392: int cmpdl_blex(int n,DL d1,DL d2)
1.1 noro 1393: {
1.53 ! noro 1394: int c;
1.1 noro 1395:
1.53 ! noro 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: }
1.1 noro 1402: }
1403:
1.19 noro 1404: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1 noro 1405: {
1.53 ! noro 1406: int e1,e2,c;
1.1 noro 1407:
1.53 ! noro 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: }
1.1 noro 1420: }
1421:
1.19 noro 1422: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1 noro 1423: {
1.53 ! noro 1424: int e1,e2,c;
1.1 noro 1425:
1.53 ! noro 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: }
1.1 noro 1438: }
1439:
1.19 noro 1440: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1 noro 1441: {
1.53 ! noro 1442: int e1,e2,f1,f2,c,i;
1.1 noro 1443:
1.53 ! noro 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: }
1.1 noro 1465: }
1466:
1.19 noro 1467: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1 noro 1468: {
1.53 ! noro 1469: int e1,e2,f1,f2,c,i;
1.1 noro 1470:
1.53 ! noro 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: }
1.1 noro 1492: }
1493:
1.19 noro 1494: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1 noro 1495: {
1.53 ! noro 1496: int e1,e2,f1,f2,c,i;
1.1 noro 1497:
1.53 ! noro 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: }
1.1 noro 1513: }
1514:
1.19 noro 1515: int cmpdl_elim(int n,DL d1,DL d2)
1.1 noro 1516: {
1.53 ! noro 1517: int e1,e2,i;
1.1 noro 1518:
1.53 ! noro 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);
1.12 noro 1528: }
1529:
1.19 noro 1530: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12 noro 1531: {
1.53 ! noro 1532: int e1,e2,i;
1.12 noro 1533:
1.53 ! noro 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);
1.13 noro 1546: }
1547:
1548: /*
1.53 ! noro 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
1.13 noro 1553: */
1554:
1.20 noro 1555: extern int *current_weyl_weight_vector;
1.13 noro 1556:
1.19 noro 1557: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13 noro 1558: {
1.53 ! noro 1559: int e1,e2,m,i;
! 1560: int *p1,*p2;
1.13 noro 1561:
1.53 ! noro 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);
1.21 noro 1587: }
1588:
1589: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1590: {
1.53 ! noro 1591: int i,t,m;
! 1592: int *p1,*p2;
1.21 noro 1593:
1.53 ! noro 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: }
1.21 noro 1606: }
1607:
1608: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1609: {
1.53 ! noro 1610: int e1,e2,m,i,t;
! 1611: int *p1,*p2;
1.21 noro 1612:
1.53 ! noro 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;
1.1 noro 1640: }
1641:
1.19 noro 1642: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1 noro 1643: {
1.53 ! noro 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;
1.28 noro 1699: }
1700:
1701: int cmpdl_composite(int nv,DL d1,DL d2)
1702: {
1.53 ! noro 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;
1.1 noro 1768: }
1769:
1.19 noro 1770: int cmpdl_matrix(int n,DL d1,DL d2)
1.1 noro 1771: {
1.53 ! noro 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;
1.25 noro 1790: }
1791:
1.43 noro 1792: int cmpdl_top_weight(int n,DL d1,DL d2)
1793: {
1.53 ! noro 1794: int *w;
! 1795: N sum,wm,wma,t;
! 1796: Q **mat;
! 1797: Q *a;
! 1798: struct oN tn;
! 1799: int len,i,sgn,tsgn,row,k;
! 1800: int *t1,*t2;
! 1801:
! 1802: w = (int *)ALLOCA(n*sizeof(int));
! 1803: len = current_top_weight_len+3;
! 1804: t1 = d1->d; t2 = d2->d;
! 1805: for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
! 1806: sum = (N)W_ALLOC(len); sgn = 0;
! 1807: wm = (N)W_ALLOC(len);
! 1808: wma = (N)W_ALLOC(len);
! 1809: if ( OID(current_top_weight) == O_VECT ) {
! 1810: mat = (Q **)&BDY((VECT)current_top_weight);
! 1811: row = 1;
! 1812: } else {
! 1813: mat = (Q **)BDY((MAT)current_top_weight);
! 1814: row = ((MAT)current_top_weight)->row;
! 1815: }
! 1816: for ( k = 0; k < row; k++ ) {
! 1817: a = mat[k];
! 1818: for ( i = 0; i < n; i++ ) {
! 1819: if ( !a[i] || !w[i] ) continue;
! 1820: tn.p = 1;
! 1821: if ( w[i] > 0 ) {
! 1822: tn.b[0] = w[i]; tsgn = 1;
! 1823: } else {
! 1824: tn.b[0] = -w[i]; tsgn = -1;
! 1825: }
! 1826: _muln(NM(a[i]),&tn,wm);
! 1827: if ( !sgn ) {
! 1828: sgn = tsgn;
! 1829: t = wm; wm = sum; sum = t;
! 1830: } else if ( sgn == tsgn ) {
! 1831: _addn(sum,wm,wma);
! 1832: if ( !PL(wma) )
! 1833: sgn = 0;
! 1834: t = wma; wma = sum; sum = t;
! 1835: } else {
! 1836: sgn *= _subn(sum,wm,wma);
! 1837: t = wma; wma = sum; sum = t;
! 1838: }
! 1839: }
! 1840: if ( sgn > 0 ) return 1;
! 1841: else if ( sgn < 0 ) return -1;
! 1842: }
1.48 noro 1843: return (*cmpdl_tie_breaker)(n,d1,d2);
1.43 noro 1844: }
1845:
1.25 noro 1846: GeoBucket create_bucket()
1847: {
1.53 ! noro 1848: GeoBucket g;
1.25 noro 1849:
1.53 ! noro 1850: g = CALLOC(1,sizeof(struct oGeoBucket));
! 1851: g->m = 32;
! 1852: return g;
1.25 noro 1853: }
1854:
1.47 noro 1855: int length(NODE d);
1856:
1.25 noro 1857: void add_bucket(GeoBucket g,NODE d,int nv)
1858: {
1.53 ! noro 1859: int l,k,m;
1.25 noro 1860:
1.53 ! noro 1861: l = length(d);
! 1862: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
! 1863: /* 2^(k-1) < l <= 2^k */
! 1864: d = symb_merge(g->body[k],d,nv);
! 1865: for ( ; length(d) > (1<<(k)); k++ ) {
! 1866: g->body[k] = 0;
! 1867: d = symb_merge(g->body[k+1],d,nv);
! 1868: }
! 1869: g->body[k] = d;
! 1870: g->m = MAX(g->m,k);
1.25 noro 1871: }
1872:
1873: DL remove_head_bucket(GeoBucket g,int nv)
1874: {
1.53 ! noro 1875: int j,i,c,m;
! 1876: DL d;
1.25 noro 1877:
1.53 ! noro 1878: j = -1;
! 1879: m = g->m;
! 1880: for ( i = 0; i <= m; i++ ) {
! 1881: if ( !g->body[i] )
! 1882: continue;
! 1883: if ( j < 0 ) j = i;
! 1884: else {
! 1885: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
! 1886: if ( c > 0 )
! 1887: j = i;
! 1888: else if ( c == 0 )
! 1889: g->body[i] = NEXT(g->body[i]);
! 1890: }
! 1891: }
! 1892: if ( j < 0 )
! 1893: return 0;
! 1894: else {
! 1895: d = g->body[j]->body;
! 1896: g->body[j] = NEXT(g->body[j]);
! 1897: return d;
! 1898: }
1.31 noro 1899: }
1900:
1901: /* DPV functions */
1902:
1903: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
1904: {
1.53 ! noro 1905: int i,len;
! 1906: DP *e;
1.31 noro 1907:
1.53 ! noro 1908: if ( !p1 || !p2 )
! 1909: error("adddv : invalid argument");
! 1910: else if ( p1->len != p2->len )
! 1911: error("adddv : size mismatch");
! 1912: else {
! 1913: len = p1->len;
! 1914: e = (DP *)MALLOC(p1->len*sizeof(DP));
! 1915: for ( i = 0; i < len; i++ )
! 1916: addd(vl,p1->body[i],p2->body[i],&e[i]);
! 1917: MKDPV(len,e,*pr);
! 1918: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 1919: }
1.31 noro 1920: }
1921:
1922: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
1923: {
1.53 ! noro 1924: int i,len;
! 1925: DP *e;
1.31 noro 1926:
1.53 ! noro 1927: if ( !p1 || !p2 )
! 1928: error("subdv : invalid argument");
! 1929: else if ( p1->len != p2->len )
! 1930: error("subdv : size mismatch");
! 1931: else {
! 1932: len = p1->len;
! 1933: e = (DP *)MALLOC(p1->len*sizeof(DP));
! 1934: for ( i = 0; i < len; i++ )
! 1935: subd(vl,p1->body[i],p2->body[i],&e[i]);
! 1936: MKDPV(len,e,*pr);
! 1937: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 1938: }
1.31 noro 1939: }
1940:
1941: void chsgndv(DPV p1,DPV *pr)
1942: {
1.53 ! noro 1943: int i,len;
! 1944: DP *e;
1.31 noro 1945:
1.53 ! noro 1946: if ( !p1 )
! 1947: error("subdv : invalid argument");
! 1948: else {
! 1949: len = p1->len;
! 1950: e = (DP *)MALLOC(p1->len*sizeof(DP));
! 1951: for ( i = 0; i < len; i++ )
! 1952: chsgnd(p1->body[i],&e[i]);
! 1953: MKDPV(len,e,*pr);
! 1954: (*pr)->sugar = p1->sugar;
! 1955: }
1.31 noro 1956: }
1957:
1958: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
1959: {
1.53 ! noro 1960: int i,len;
! 1961: DP *e;
1.31 noro 1962:
1.53 ! noro 1963: len = p2->len;
! 1964: e = (DP *)MALLOC(p2->len*sizeof(DP));
! 1965: if ( !p1 ) {
! 1966: MKDPV(len,e,*pr);
! 1967: (*pr)->sugar = 0;
! 1968: } else {
! 1969: for ( i = 0; i < len; i++ )
! 1970: muld(vl,p1,p2->body[i],&e[i]);
! 1971: MKDPV(len,e,*pr);
! 1972: (*pr)->sugar = p1->sugar + p2->sugar;
! 1973: }
1.31 noro 1974: }
1975:
1976: int compdv(VL vl,DPV p1,DPV p2)
1977: {
1.53 ! noro 1978: int i,t,len;
1.31 noro 1979:
1.53 ! noro 1980: if ( p1->len != p2->len ) {
! 1981: error("compdv : size mismatch");
! 1982: return 0; /* XXX */
! 1983: } else {
! 1984: len = p1->len;
! 1985: for ( i = 0; i < len; i++ )
! 1986: if ( (t = compd(vl,p1->body[i],p2->body[i])) )
! 1987: return t;
! 1988: return 0;
! 1989: }
1.33 noro 1990: }
1991:
1992: int ni_next(int *a,int n)
1993: {
1.53 ! noro 1994: int i,j,k,kj;
1.33 noro 1995:
1.53 ! noro 1996: /* find the first nonzero a[j] */
! 1997: for ( j = 0; j < n && a[j] == 0; j++ );
! 1998: /* find the first zero a[k] after a[j] */
! 1999: for ( k = j; k < n && a[k] == 1; k++ );
! 2000: if ( k == n ) return 0;
! 2001: /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
! 2002: /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
! 2003: kj = k-j-1;
! 2004: for ( i = 0; i < kj; i++ ) a[i] = 1;
! 2005: for ( ; i < k; i++ ) a[i] = 0;
! 2006: a[k] = 1;
! 2007: return 1;
1.33 noro 2008: }
2009:
2010: int comp_nbm(NBM a,NBM b)
2011: {
1.53 ! noro 2012: int d,i,ai,bi;
! 2013: int *ab,*bb;
1.33 noro 2014:
1.53 ! noro 2015: if ( a->d > b->d ) return 1;
! 2016: else if ( a->d < b->d ) return -1;
! 2017: else {
! 2018: d = a->d; ab = a->b; bb = b->b;
1.41 noro 2019: #if 0
1.53 ! noro 2020: w = (d+31)/32;
! 2021: for ( i = 0; i < w; i++ )
! 2022: if ( ab[i] > bb[i] ) return 1;
! 2023: else if ( ab[i] < bb[i] ) return -1;
1.41 noro 2024: #else
1.53 ! noro 2025: for ( i = 0; i < d; i++ ) {
! 2026: ai = NBM_GET(ab,i);
! 2027: bi = NBM_GET(bb,i);
! 2028: if ( ai > bi ) return 1;
! 2029: else if ( ai < bi ) return -1;
! 2030: }
1.41 noro 2031: #endif
1.53 ! noro 2032: return 0;
! 2033: }
1.33 noro 2034: }
2035:
2036: NBM mul_nbm(NBM a,NBM b)
2037: {
1.53 ! noro 2038: int ad,bd,d,i,j;
! 2039: int *ab,*bb,*mb;
! 2040: NBM m;
! 2041:
! 2042: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
! 2043: d = ad + bd;
! 2044: NEWNBM(m); NEWNBMBDY(m,d);
! 2045: m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
! 2046: j = 0;
! 2047: for ( i = 0; i < ad; i++, j++ )
! 2048: if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
! 2049: else NBM_CLR(mb,j);
! 2050: for ( i = 0; i < bd; i++, j++ )
! 2051: if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
! 2052: else NBM_CLR(mb,j);
! 2053: return m;
1.33 noro 2054: }
2055:
1.37 noro 2056: NBP nbmtonbp(NBM m)
2057: {
1.53 ! noro 2058: NODE n;
! 2059: NBP u;
1.37 noro 2060:
1.53 ! noro 2061: MKNODE(n,m,0);
! 2062: MKNBP(u,n);
! 2063: return u;
1.37 noro 2064: }
2065:
2066: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
2067:
1.40 noro 2068: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 2069: {
1.53 ! noro 2070: int i,d1;
! 2071: NBM t;
1.37 noro 2072:
1.53 ! noro 2073: if ( !a->d ) error("separate_nbm : invalid argument");
1.37 noro 2074:
1.53 ! noro 2075: if ( a0 ) {
! 2076: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
! 2077: *a0 = nbmtonbp(t);
! 2078: }
! 2079:
! 2080: if ( ah ) {
! 2081: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
! 2082: if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
! 2083: else NBM_CLR(t->b,0);
! 2084: *ah = nbmtonbp(t);
! 2085: }
1.42 noro 2086:
1.53 ! noro 2087: if ( ar ) {
! 2088: d1 = a->d-1;
! 2089: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
! 2090: for ( i = 0; i < d1; i++ ) {
! 2091: if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
! 2092: else NBM_CLR(t->b,i);
! 2093: }
! 2094: *ar = nbmtonbp(t);
! 2095: }
! 2096:
! 2097: return a->c;
1.42 noro 2098: }
2099:
2100: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
2101:
2102: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
2103: {
1.53 ! noro 2104: int i,d,d1;
! 2105: NBM t;
! 2106:
! 2107: if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
1.42 noro 2108:
1.53 ! noro 2109: if ( a0 ) {
! 2110: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
! 2111: *a0 = nbmtonbp(t);
! 2112: }
! 2113:
! 2114: d1 = a->d-1;
! 2115: if ( at ) {
! 2116: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
! 2117: if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
! 2118: else NBM_CLR(t->b,0);
! 2119: *at = nbmtonbp(t);
! 2120: }
1.42 noro 2121:
1.53 ! noro 2122: if ( ar ) {
! 2123: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
! 2124: for ( i = 0; i < d1; i++ ) {
! 2125: if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
! 2126: else NBM_CLR(t->b,i);
! 2127: }
! 2128: *ar = nbmtonbp(t);
! 2129: }
1.37 noro 2130:
1.53 ! noro 2131: return a->c;
1.37 noro 2132: }
2133:
2134: NBP make_xky(int k)
2135: {
1.53 ! noro 2136: int k1,i;
! 2137: NBM t;
1.37 noro 2138:
1.53 ! noro 2139: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
! 2140: k1 = k-1;
! 2141: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
! 2142: NBM_CLR(t->b,i);
! 2143: return nbmtonbp(t);
1.37 noro 2144: }
2145:
2146: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
2147:
1.40 noro 2148: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 2149: {
1.53 ! noro 2150: int i,d1,k,k1;
! 2151: NBM t;
1.37 noro 2152:
1.53 ! noro 2153: if ( !a->d )
! 2154: error("separate_nbm : invalid argument");
! 2155: for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
! 2156: if ( i == a->d )
! 2157: error("separate_nbm : invalid argument");
! 2158: k1 = i;
! 2159: k = i+1;
! 2160:
! 2161: if ( a0 ) {
! 2162: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
! 2163: *a0 = nbmtonbp(t);
! 2164: }
! 2165:
! 2166: if ( ah ) {
! 2167: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
! 2168: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
! 2169: NBM_CLR(t->b,i);
! 2170: *ah = nbmtonbp(t);
! 2171: }
1.37 noro 2172:
1.53 ! noro 2173: if ( ar ) {
! 2174: d1 = a->d-k;
! 2175: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
! 2176: for ( i = 0; i < d1; i++ ) {
! 2177: if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
! 2178: else NBM_CLR(t->b,i);
! 2179: }
! 2180: *ar = nbmtonbp(t);
! 2181: }
! 2182:
! 2183: return a->c;
1.37 noro 2184: }
1.33 noro 2185:
1.37 noro 2186: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2187: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
1.38 noro 2188: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
2189: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
1.37 noro 2190:
2191: NBP shuffle_mul_nbm(NBM a,NBM b)
2192: {
1.53 ! noro 2193: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
! 2194: P ac,bc,c;
1.37 noro 2195:
1.53 ! noro 2196: if ( !a->d || !b->d )
! 2197: u = nbmtonbp(mul_nbm(a,b));
! 2198: else {
! 2199: ac = separate_nbm(a,&a0,&ah,&ar);
! 2200: bc = separate_nbm(b,&b0,&bh,&br);
! 2201: mulp(CO,ac,bc,&c);
! 2202: shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
! 2203: shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
! 2204: addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
! 2205: }
! 2206: return u;
1.37 noro 2207: }
1.33 noro 2208:
1.37 noro 2209: NBP harmonic_mul_nbm(NBM a,NBM b)
2210: {
1.53 ! noro 2211: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
! 2212: P ac,bc,c;
1.37 noro 2213:
1.53 ! noro 2214: if ( !a->d || !b->d )
! 2215: u = nbmtonbp(mul_nbm(a,b));
! 2216: else {
! 2217: mulp(CO,a->c,b->c,&c);
! 2218: ac = separate_xky_nbm(a,&a0,&ah,&ar);
! 2219: bc = separate_xky_nbm(b,&b0,&bh,&br);
! 2220: mulp(CO,ac,bc,&c);
! 2221: harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
! 2222: harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
! 2223: abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
! 2224: harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
! 2225: addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
! 2226: }
! 2227: return u;
1.37 noro 2228:
2229: }
1.34 noro 2230:
1.33 noro 2231: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2232: {
1.53 ! noro 2233: NODE b1,b2,br=0,br0;
! 2234: NBM m1,m2,m;
! 2235: P c;
! 2236:
! 2237: if ( !p1 )
! 2238: *rp = p2;
! 2239: else if ( !p2 )
! 2240: *rp = p1;
! 2241: else {
! 2242: for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
! 2243: m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
! 2244: switch ( comp_nbm(m1,m2) ) {
! 2245: case 0:
! 2246: addp(CO,m1->c,m2->c,&c);
! 2247: if ( c ) {
! 2248: NEXTNODE(br0,br);
! 2249: NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
! 2250: BDY(br) = (pointer)m;
! 2251: }
! 2252: b1 = NEXT(b1); b2 = NEXT(b2); break;
! 2253: case 1:
! 2254: NEXTNODE(br0,br); BDY(br) = BDY(b1);
! 2255: b1 = NEXT(b1); break;
! 2256: case -1:
! 2257: NEXTNODE(br0,br); BDY(br) = BDY(b2);
! 2258: b2 = NEXT(b2); break;
! 2259: }
! 2260: }
! 2261: if ( !br0 )
! 2262: if ( b1 )
! 2263: br0 = b1;
! 2264: else if ( b2 )
! 2265: br0 = b2;
! 2266: else {
! 2267: *rp = 0;
! 2268: return;
! 2269: }
! 2270: else if ( b1 )
! 2271: NEXT(br) = b1;
! 2272: else if ( b2 )
! 2273: NEXT(br) = b2;
! 2274: else
! 2275: NEXT(br) = 0;
! 2276: MKNBP(*rp,br0);
! 2277: }
1.33 noro 2278: }
2279:
2280: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2281: {
1.53 ! noro 2282: NBP t;
1.33 noro 2283:
1.53 ! noro 2284: chsgnnbp(p2,&t);
! 2285: addnbp(vl,p1,t,rp);
1.33 noro 2286: }
2287:
2288: void chsgnnbp(NBP p,NBP *rp)
2289: {
1.53 ! noro 2290: NODE r0,r=0,b;
! 2291: NBM m,m1;
1.33 noro 2292:
1.53 ! noro 2293: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
! 2294: NEXTNODE(r0,r);
! 2295: m = (NBM)BDY(b);
! 2296: NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
! 2297: BDY(r) = m1;
! 2298: }
! 2299: if ( r0 ) NEXT(r) = 0;
! 2300: MKNBP(*rp,r0);
1.33 noro 2301: }
2302:
2303: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2304: {
1.53 ! noro 2305: NODE b,n;
! 2306: NBP r,t,s;
! 2307: NBM m;
! 2308:
! 2309: if ( !p1 || !p2 ) {
! 2310: *rp = 0; return;
! 2311: }
! 2312: if ( OID(p1) != O_NBP ) {
! 2313: if ( !POLY(p1) )
! 2314: error("mulnbp : invalid argument");
! 2315: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
! 2316: MKNODE(n,m,0); MKNBP(p1,n);
! 2317: }
! 2318: if ( OID(p2) != O_NBP ) {
! 2319: if ( !POLY(p2) )
! 2320: error("mulnbp : invalid argument");
! 2321: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
! 2322: MKNODE(n,m,0); MKNBP(p2,n);
! 2323: }
! 2324: if ( length(BDY(p1)) < length(BDY(p2)) ) {
! 2325: for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
! 2326: mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
! 2327: addnbp(vl,r,t,&s); r = s;
! 2328: }
! 2329: *rp = r;
! 2330: } else {
! 2331: for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
! 2332: mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
! 2333: addnbp(vl,r,t,&s); r = s;
! 2334: }
! 2335: *rp = r;
! 2336: }
1.33 noro 2337: }
2338:
2339: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
2340: {
1.53 ! noro 2341: NODE b,r0,r=0;
1.33 noro 2342:
1.53 ! noro 2343: if ( !p ) *rp = 0;
! 2344: else {
! 2345: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
! 2346: NEXTNODE(r0,r);
! 2347: BDY(r) = mul_nbm(m,(NBM)BDY(b));
! 2348: }
! 2349: if ( r0 ) NEXT(r) = 0;
! 2350: MKNBP(*rp,r0);
! 2351: }
1.33 noro 2352: }
2353:
2354: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
2355: {
1.53 ! noro 2356: NODE b,r0,r=0;
1.33 noro 2357:
1.53 ! noro 2358: if ( !p ) *rp = 0;
! 2359: else {
! 2360: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
! 2361: NEXTNODE(r0,r);
! 2362: BDY(r) = mul_nbm((NBM)BDY(b),m);
! 2363: }
! 2364: if ( r0 ) NEXT(r) = 0;
! 2365: MKNBP(*rp,r0);
! 2366: }
1.33 noro 2367: }
2368:
2369: void pwrnbp(VL vl,NBP a,Q q,NBP *c)
2370: {
1.53 ! noro 2371: int t;
! 2372: NBP a1,a2;
! 2373: N n1;
! 2374: Q q1;
! 2375: NBM m;
! 2376: NODE r;
! 2377:
! 2378: if ( !q ) {
! 2379: NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
! 2380: MKNODE(r,m,0); MKNBP(*c,r);
! 2381: } else if ( !a )
! 2382: *c = 0;
! 2383: else if ( UNIQ(q) )
! 2384: *c = a;
! 2385: else {
! 2386: t = divin(NM(q),2,&n1); NTOQ(n1,1,q1);
! 2387: pwrnbp(vl,a,q1,&a1);
! 2388: mulnbp(vl,a1,a1,&a2);
! 2389: if ( t )
! 2390: mulnbp(vl,a2,a,c);
! 2391: else
! 2392: *c = a2;
! 2393: }
1.33 noro 2394: }
2395:
1.38 noro 2396: int compnbp(VL vl,NBP p1,NBP p2)
2397: {
1.53 ! noro 2398: NODE n1,n2;
! 2399: NBM m1,m2;
! 2400: int t;
! 2401:
! 2402: if ( !p1 )
! 2403: return p2 ? -1 : 0;
! 2404: else if ( !p2 )
! 2405: return 1;
! 2406: else {
! 2407: for ( n1 = BDY(p1), n2 = BDY(p2);
! 2408: n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
! 2409: m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
! 2410: if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
! 2411: return t;
! 2412: }
! 2413: if ( n1 )
! 2414: return 1;
! 2415: else if ( n2 )
! 2416: return -1;
! 2417: else
! 2418: return 0;
! 2419: }
1.38 noro 2420: }
2421:
1.33 noro 2422: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2423: {
1.53 ! noro 2424: NODE b1,b2,n;
! 2425: NBP r,t,s;
! 2426: NBM m;
! 2427:
! 2428: if ( !p1 || !p2 ) {
! 2429: *rp = 0; return;
! 2430: }
! 2431: if ( OID(p1) != O_NBP ) {
! 2432: if ( !POLY(p1) )
! 2433: error("shuffle_mulnbp : invalid argument");
! 2434: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
! 2435: MKNODE(n,m,0); MKNBP(p1,n);
! 2436: }
! 2437: if ( OID(p2) != O_NBP ) {
! 2438: if ( !POLY(p2) )
! 2439: error("shuffle_mulnbp : invalid argument");
! 2440: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
! 2441: MKNODE(n,m,0); MKNBP(p2,n);
! 2442: }
! 2443: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
! 2444: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
! 2445: t = shuffle_mul_nbm(m,(NBM)BDY(b2));
! 2446: addnbp(vl,r,t,&s); r = s;
! 2447: }
! 2448: *rp = r;
1.33 noro 2449: }
2450:
1.34 noro 2451: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
1.33 noro 2452: {
1.53 ! noro 2453: NODE b1,b2,n;
! 2454: NBP r,t,s;
! 2455: NBM m;
! 2456:
! 2457: if ( !p1 || !p2 ) {
! 2458: *rp = 0; return;
! 2459: }
! 2460: if ( OID(p1) != O_NBP ) {
! 2461: if ( !POLY(p1) )
! 2462: error("harmonic_mulnbp : invalid argument");
! 2463: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
! 2464: MKNODE(n,m,0); MKNBP(p1,n);
! 2465: }
! 2466: if ( OID(p2) != O_NBP ) {
! 2467: if ( !POLY(p2) )
! 2468: error("harmonic_mulnbp : invalid argument");
! 2469: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
! 2470: MKNODE(n,m,0); MKNBP(p2,n);
! 2471: }
! 2472: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
! 2473: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
! 2474: t = harmonic_mul_nbm(m,(NBM)BDY(b2));
! 2475: addnbp(vl,r,t,&s); r = s;
! 2476: }
! 2477: *rp = r;
1.1 noro 2478: }
1.38 noro 2479:
2480: #if 0
2481: NBP shuffle_mul_nbm(NBM a,NBM b)
2482: {
1.53 ! noro 2483: int ad,bd,d,i,ai,bi,bit,s;
! 2484: int *ab,*bb,*wmb,*w;
! 2485: NBM wm,tm;
! 2486: P c,c1;
! 2487: NODE r,t,t1,p;
! 2488: NBP u;
! 2489:
! 2490: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
! 2491: d = ad + bd;
! 2492: w = (int *)ALLOCA(d*sizeof(int));
! 2493: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2494: for ( i = 0; i < ad; i++ ) w[i] = 1;
! 2495: for ( ; i < d; i++ ) w[i] = 0;
! 2496: mulp(CO,a->c,b->c,&c);
! 2497: r = 0;
! 2498: do {
! 2499: wm->d = d; wm->c = c;
! 2500: ai = 0; bi = 0;
! 2501: for ( i = 0; i < d; i++ ) {
! 2502: if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
! 2503: else { bit = NBM_GET(bb,bi); bi++; }
! 2504: if ( bit ) NBM_SET(wmb,i);
! 2505: else NBM_CLR(wmb,i);
! 2506: }
! 2507: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
! 2508: tm = (NBM)BDY(t);
! 2509: s = comp_nbm(tm,wm);
! 2510: if ( s < 0 ) {
! 2511: /* insert */
! 2512: MKNODE(t1,wm,t);
! 2513: if ( !p ) r = t1;
! 2514: else NEXT(p) = t1;
! 2515: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2516: break;
! 2517: } else if ( s == 0 ) {
! 2518: /* add coefs */
! 2519: addp(CO,tm->c,c,&c1);
! 2520: if ( c1 ) tm->c = c1;
! 2521: else NEXT(p) = NEXT(t);
! 2522: break;
! 2523: }
! 2524: }
! 2525: if ( !t ) {
! 2526: /* append */
! 2527: MKNODE(t1,wm,t);
! 2528: if ( !p ) r = t1;
! 2529: else NEXT(p) = t1;
! 2530: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2531: }
! 2532: } while ( ni_next(w,d) );
! 2533: MKNBP(u,r);
! 2534: return u;
1.38 noro 2535: }
2536:
2537: int nbmtoxky(NBM a,int *b)
2538: {
1.53 ! noro 2539: int d,i,j,k;
! 2540: int *p;
1.38 noro 2541:
1.53 ! noro 2542: d = a->d; p = a->b;
! 2543: for ( i = j = 0, k = 1; i < d; i++ ) {
! 2544: if ( !NBM_GET(p,i) ) {
! 2545: b[j++] = k;
! 2546: k = 1;
! 2547: } else k++;
! 2548: }
! 2549: return j;
1.38 noro 2550: }
2551:
2552: NBP harmonic_mul_nbm(NBM a,NBM b)
2553: {
1.53 ! noro 2554: int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
! 2555: int i,j,k,ia,ib,s;
! 2556: int *wa,*wb,*w,*wab,*wa1,*wmb;
! 2557: P c,c1;
! 2558: NBM wm,tm;
! 2559: NODE r,t1,t,p;
! 2560: NBP u;
! 2561:
! 2562: da = a->d; db = b->d; d = da+db;
! 2563: wa = (int *)ALLOCA(da*sizeof(int));
! 2564: wb = (int *)ALLOCA(db*sizeof(int));
! 2565: la = nbmtoxky(a,wa);
! 2566: lb = nbmtoxky(b,wb);
! 2567: mulp(CO,a->c,b->c,&c);
! 2568: /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
! 2569: /* lmax : total length */
! 2570: lmax = la+lb;
! 2571: lmin = la>lb?la:lb;
! 2572: w = (int *)ALLOCA(lmax*sizeof(int));
! 2573: /* position of a+b */
! 2574: wab = (int *)ALLOCA(lmax*sizeof(int));
! 2575: /* position of a */
! 2576: wa1 = (int *)ALLOCA(lmax*sizeof(int));
! 2577: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2578: for ( l = lmin, r = 0; l <= lmax; l++ ) {
! 2579: lab = lmax - l;
! 2580: la1 = la - lab;
! 2581: lb1 = lb - lab;
! 2582: lab1 = l-lab;
! 2583: /* partion l into three parts: a, b, a+b */
! 2584: /* initialize wab */
! 2585: for ( i = 0; i < lab; i++ ) wab[i] = 1;
! 2586: for ( ; i < l; i++ ) wab[i] = 0;
! 2587: do {
! 2588: /* initialize wa1 */
! 2589: for ( i = 0; i < la1; i++ ) wa1[i] = 1;
! 2590: for ( ; i < lab1; i++ ) wa1[i] = 0;
! 2591: do {
! 2592: ia = 0; ib = 0;
! 2593: for ( i = j = 0; i < l; i++ )
! 2594: if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
! 2595: else if ( wa1[j++] ) w[i] = wa[ia++];
! 2596: else w[i] = wb[ib++];
! 2597: for ( i = j = 0; i < l; i++ ) {
! 2598: for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
! 2599: NBM_CLR(wmb,j); j++;
! 2600: }
! 2601: wm->d = j; wm->c = c;
! 2602: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
! 2603: tm = (NBM)BDY(t);
! 2604: s = comp_nbm(tm,wm);
! 2605: if ( s < 0 ) {
! 2606: /* insert */
! 2607: MKNODE(t1,wm,t);
! 2608: if ( !p ) r = t1;
! 2609: else NEXT(p) = t1;
! 2610: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2611: break;
! 2612: } else if ( s == 0 ) {
! 2613: /* add coefs */
! 2614: addp(CO,tm->c,c,&c1);
! 2615: if ( c1 ) tm->c = c1;
! 2616: else NEXT(p) = NEXT(t);
! 2617: break;
! 2618: }
! 2619: }
! 2620: if ( !t ) {
! 2621: /* append */
! 2622: MKNODE(t1,wm,t);
! 2623: if ( !p ) r = t1;
! 2624: else NEXT(p) = t1;
! 2625: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
! 2626: }
! 2627: } while ( ni_next(wa1,lab1) );
! 2628: } while ( ni_next(wab,l) );
! 2629: }
! 2630: MKNBP(u,r);
! 2631: return u;
1.38 noro 2632: }
2633: #endif
1.52 noro 2634:
2635: /* DPM functions */
2636:
2637: int compdmm(int n,DMM m1,DMM m2)
2638: {
2639: int t;
2640:
2641: if ( dpm_ispot ) {
2642: if ( m1->pos < m2->pos ) return 1;
2643: else if ( m1->pos > m2->pos ) return -1;
2644: else return (*cmpdl)(n,m1->dl,m2->dl);
2645: } else {
2646: t = (*cmpdl)(n,m1->dl,m2->dl);
2647: if ( t ) return t;
2648: else if ( m1->pos < m2->pos ) return 1;
2649: else if ( m1->pos > m2->pos ) return -1;
2650: else return 0;
2651: }
2652: }
2653:
2654: void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
2655: {
1.53 ! noro 2656: int n;
! 2657: DMM m1,m2,mr=0,mr0;
! 2658: Obj t;
! 2659: DL d;
! 2660:
! 2661: if ( !p1 )
! 2662: *pr = p2;
! 2663: else if ( !p2 )
! 2664: *pr = p1;
! 2665: else {
! 2666: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
! 2667: switch ( compdmm(n,m1,m2) ) {
! 2668: case 0:
! 2669: arf_add(vl,C(m1),C(m2),&t);
! 2670: if ( t ) {
! 2671: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = t;
! 2672: }
! 2673: m1 = NEXT(m1); m2 = NEXT(m2); break;
! 2674: case 1:
! 2675: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = C(m1);
! 2676: m1 = NEXT(m1); break;
! 2677: case -1:
! 2678: NEXTDMM(mr0,mr); mr->pos = m2->pos; mr->dl = m2->dl; C(mr) = C(m2);
! 2679: m2 = NEXT(m2); break;
! 2680: }
! 2681: if ( !mr0 )
! 2682: if ( m1 )
! 2683: mr0 = m1;
! 2684: else if ( m2 )
! 2685: mr0 = m2;
! 2686: else {
! 2687: *pr = 0;
! 2688: return;
! 2689: }
! 2690: else if ( m1 )
! 2691: NEXT(mr) = m1;
! 2692: else if ( m2 )
! 2693: NEXT(mr) = m2;
! 2694: else
! 2695: NEXT(mr) = 0;
! 2696: MKDPM(NV(p1),mr0,*pr);
! 2697: if ( *pr )
! 2698: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
! 2699: }
1.52 noro 2700: }
2701:
2702: void subdpm(VL vl,DPM p1,DPM p2,DPM *pr)
2703: {
1.53 ! noro 2704: DPM t;
1.52 noro 2705:
1.53 ! noro 2706: if ( !p2 )
! 2707: *pr = p1;
! 2708: else {
! 2709: chsgndpm(p2,&t); adddpm(vl,p1,t,pr);
! 2710: }
1.52 noro 2711: }
2712:
2713: void chsgndpm(DPM p,DPM *pr)
2714: {
1.53 ! noro 2715: DMM m,mr=0,mr0;
! 2716: Obj r;
1.52 noro 2717:
1.53 ! noro 2718: if ( !p )
! 2719: *pr = 0;
1.52 noro 2720: else {
1.53 ! noro 2721: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2722: NEXTDMM(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->pos = m->pos; mr->dl = m->dl;
! 2723: }
! 2724: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
! 2725: if ( *pr )
! 2726: (*pr)->sugar = p->sugar;
! 2727: }
1.52 noro 2728: }
2729:
2730: void mulcdpm(VL vl,Obj c,DPM p,DPM *pr)
2731: {
1.53 ! noro 2732: DMM m,mr=0,mr0;
1.52 noro 2733:
1.53 ! noro 2734: if ( !p || !c )
! 2735: *pr = 0;
! 2736: else if ( NUM(c) && UNIQ((Q)c) )
! 2737: *pr = p;
! 2738: else if ( NUM(c) && MUNIQ((Q)c) )
! 2739: chsgndpm(p,pr);
! 2740: else {
! 2741: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2742: NEXTDMM(mr0,mr);
! 2743: arf_mul(vl,C(m),c,&C(mr));
1.52 noro 2744: mr->pos = m->pos;
1.53 ! noro 2745: mr->dl = m->dl;
! 2746: }
! 2747: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
! 2748: if ( *pr )
! 2749: (*pr)->sugar = p->sugar;
! 2750: }
1.52 noro 2751: }
2752:
2753: void comm_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
2754: {
1.53 ! noro 2755: DMM m,mr=0,mr0;
1.52 noro 2756: DL d;
2757: Obj c;
2758: int n;
2759:
1.53 ! noro 2760: if ( !p )
! 2761: *pr = 0;
! 2762: else {
1.52 noro 2763: n = NV(p);
2764: d = m0->dl;
2765: c = C(m0);
1.53 ! noro 2766: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2767: NEXTDMM(mr0,mr);
! 2768: arf_mul(vl,C(m),c,&C(mr));
1.52 noro 2769: mr->pos = m->pos;
1.53 ! noro 2770: adddl(n,m->dl,d,&mr->dl);
! 2771: }
! 2772: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
! 2773: if ( *pr )
! 2774: (*pr)->sugar = p->sugar;
! 2775: }
1.52 noro 2776: }
2777:
2778: void weyl_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
2779: {
1.53 ! noro 2780: DPM r,t,t1;
! 2781: DMM m;
! 2782: DL d0;
! 2783: int n,n2,l,i,j,tlen;
1.52 noro 2784: struct oMP mp;
1.53 ! noro 2785: static DMM *w,*psum;
! 2786: static struct cdl *tab;
! 2787: static int wlen;
! 2788: static int rtlen;
! 2789:
! 2790: if ( !p )
! 2791: *pr = 0;
! 2792: else {
! 2793: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
! 2794: if ( l > wlen ) {
! 2795: if ( w ) GCFREE(w);
! 2796: w = (DMM *)MALLOC(l*sizeof(DMM));
! 2797: wlen = l;
! 2798: }
! 2799: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
! 2800: w[i] = m;
! 2801:
! 2802: n = NV(p); n2 = n>>1;
! 2803: d0 = m0->dl;
! 2804: for ( i = 0, tlen = 1; i < n2; i++ )
! 2805: tlen *= d0->d[n2+i]+1;
! 2806: if ( tlen > rtlen ) {
! 2807: if ( tab ) GCFREE(tab);
! 2808: if ( psum ) GCFREE(psum);
! 2809: rtlen = tlen;
! 2810: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
! 2811: psum = (DMM *)MALLOC(rtlen*sizeof(DMM));
! 2812: }
! 2813: bzero(psum,tlen*sizeof(DMM));
! 2814: for ( i = l-1; i >= 0; i-- ) {
! 2815: bzero(tab,tlen*sizeof(struct cdl));
1.52 noro 2816: mp.dl = w[i]->dl; mp.c = C(w[i]); mp.next = 0;
1.53 ! noro 2817: weyl_mulmm(vl,m0,&mp,n,tab,tlen);
! 2818: for ( j = 0; j < tlen; j++ ) {
! 2819: if ( tab[j].c ) {
! 2820: NEWDMM(m); m->dl = tab[j].d; m->pos = w[i]->pos; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
! 2821: psum[j] = m;
! 2822: }
! 2823: }
! 2824: }
! 2825: for ( j = tlen-1, r = 0; j >= 0; j-- )
! 2826: if ( psum[j] ) {
! 2827: MKDPM(n,psum[j],t); adddpm(vl,r,t,&t1); r = t1;
! 2828: }
! 2829: if ( r )
! 2830: r->sugar = p->sugar + m0->dl->td;
! 2831: *pr = r;
! 2832: }
1.52 noro 2833: }
2834:
2835: void mulobjdpm(VL vl,Obj p1,DPM p2,DPM *pr)
2836: {
1.53 ! noro 2837: MP m;
! 2838: DPM s,t,u;
1.52 noro 2839:
1.53 ! noro 2840: if ( !p1 || !p2 )
! 2841: *pr = 0;
! 2842: else if ( OID(p1) != O_DP )
! 2843: mulcdpm(vl,p1,p2,pr);
! 2844: else {
1.52 noro 2845: s = 0;
2846: for ( m = BDY((DP)p1); m; m = NEXT(m) ) {
2847: if ( do_weyl )
1.53 ! noro 2848: weyl_mulmpdpm(vl,m,p2,&t);
1.52 noro 2849: else
1.53 ! noro 2850: comm_mulmpdpm(vl,m,p2,&t);
1.52 noro 2851: adddpm(vl,s,t,&u); s = u;
1.53 ! noro 2852: }
! 2853: *pr = s;
! 2854: }
1.52 noro 2855: }
2856:
2857: int compdpm(VL vl,DPM p1,DPM p2)
2858: {
1.53 ! noro 2859: int n,t;
! 2860: DMM m1,m2;
1.52 noro 2861:
1.53 ! noro 2862: if ( !p1 )
! 2863: return p2 ? -1 : 0;
! 2864: else if ( !p2 )
! 2865: return 1;
! 2866: else if ( NV(p1) != NV(p2) ) {
! 2867: error("compdpm : size mismatch");
! 2868: return 0; /* XXX */
! 2869: } else {
! 2870: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
! 2871: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
! 2872: if ( (t = compdmm(n,m1,m2)) ||
! 2873: (t = arf_comp(vl,C(m1),C(m2)) ) )
! 2874: return t;
! 2875: if ( m1 )
! 2876: return 1;
! 2877: else if ( m2 )
! 2878: return -1;
! 2879: else
! 2880: return 0;
! 2881: }
1.52 noro 2882: }
2883:
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 / asir2000 / engine