Annotation of OpenXM_contrib2/asir2000/builtin/dp-supp.c, Revision 1.67
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.67 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.66 2017/08/31 02:36:20 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
1.16 noro 52: #include "inline.h"
1.1 noro 53: #include "parse.h"
54: #include "ox.h"
55:
1.66 noro 56: #define HMAG(p) (p_mag((P)BDY(p)->c))
1.5 noro 57:
1.1 noro 58: extern int (*cmpdl)();
1.5 noro 59: extern double pz_t_e,pz_t_d,pz_t_d1,pz_t_c;
60: extern int dp_nelim,dp_fcoeffs;
1.7 noro 61: extern int NoGCD;
62: extern int GenTrace;
63: extern NODE TraceList;
64:
1.37 noro 65: int show_orderspec;
66:
67: void print_composite_order_spec(struct order_spec *spec);
1.66 noro 68: void dpm_rest(DPM,DPM *);
1.37 noro 69:
1.7 noro 70: /*
71: * content reduction
72: *
73: */
74:
1.46 noro 75: static NODE RatDenomList;
76:
77: void init_denomlist()
78: {
1.67 ! noro 79: RatDenomList = 0;
1.46 noro 80: }
81:
82: void add_denomlist(P f)
83: {
1.67 ! noro 84: NODE n;
1.46 noro 85:
1.67 ! noro 86: if ( OID(f)==O_P ) {
! 87: MKNODE(n,f,RatDenomList); RatDenomList = n;
! 88: }
1.46 noro 89: }
90:
91: LIST get_denomlist()
92: {
1.67 ! noro 93: LIST l;
1.46 noro 94:
1.67 ! noro 95: MKLIST(l,RatDenomList); RatDenomList = 0;
! 96: return l;
1.46 noro 97: }
98:
1.20 noro 99: void dp_ptozp(DP p,DP *rp)
1.7 noro 100: {
1.67 ! noro 101: MP m,mr,mr0;
! 102: int i,n;
! 103: Q *w;
! 104: Q dvr;
! 105: P t;
! 106:
! 107: if ( !p )
! 108: *rp = 0;
! 109: else {
! 110: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
! 111: w = (Q *)ALLOCA(n*sizeof(Q));
! 112: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
! 113: if ( NUM(m->c) )
! 114: w[i] = (Q)m->c;
! 115: else
! 116: ptozp((P)m->c,1,&w[i],&t);
! 117: sortbynm(w,n);
! 118: qltozl(w,n,&dvr);
! 119: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 120: NEXTMP(mr0,mr); divsp(CO,(P)m->c,(P)dvr,(P *)&mr->c); mr->dl = m->dl;
! 121: }
! 122: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 123: }
1.7 noro 124: }
125:
1.20 noro 126: void dp_ptozp2(DP p0,DP p1,DP *hp,DP *rp)
1.7 noro 127: {
1.67 ! noro 128: DP t,s,h,r;
! 129: MP m,mr,mr0,m0;
1.7 noro 130:
1.67 ! noro 131: addd(CO,p0,p1,&t); dp_ptozp(t,&s);
! 132: if ( !p0 ) {
! 133: h = 0; r = s;
! 134: } else if ( !p1 ) {
! 135: h = s; r = 0;
! 136: } else {
! 137: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
! 138: m = NEXT(m), m0 = NEXT(m0) ) {
! 139: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
! 140: }
! 141: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
! 142: }
! 143: if ( h )
! 144: h->sugar = p0->sugar;
! 145: if ( r )
! 146: r->sugar = p1->sugar;
! 147: *hp = h; *rp = r;
1.39 ohara 148: }
149:
1.66 noro 150: void dpm_ptozp(DPM p,DPM *rp)
151: {
1.67 ! noro 152: DMM m,mr,mr0;
! 153: int i,n;
! 154: Q *w;
! 155: Q dvr;
! 156: P t;
! 157:
! 158: if ( !p )
! 159: *rp = 0;
! 160: else {
! 161: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
! 162: w = (Q *)ALLOCA(n*sizeof(Q));
! 163: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
! 164: if ( NUM(m->c) )
! 165: w[i] = (Q)m->c;
! 166: else
! 167: ptozp((P)m->c,1,&w[i],&t);
! 168: sortbynm(w,n);
! 169: qltozl(w,n,&dvr);
! 170: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 171: NEXTDMM(mr0,mr); divsp(CO,(P)m->c,(P)dvr,(P *)&mr->c); mr->dl = m->dl; mr->pos = m->pos;
! 172: }
! 173: NEXT(mr) = 0; MKDPM(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 174: }
1.66 noro 175: }
176:
177: void dpm_ptozp2(DPM p0,DPM p1,DPM *hp,DPM *rp)
178: {
1.67 ! noro 179: DPM t,s,h,r;
! 180: DMM m,mr,mr0,m0;
1.66 noro 181:
1.67 ! noro 182: adddpm(CO,p0,p1,&t); dpm_ptozp(t,&s);
! 183: if ( !p0 ) {
! 184: h = 0; r = s;
! 185: } else if ( !p1 ) {
! 186: h = s; r = 0;
! 187: } else {
! 188: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
! 189: m = NEXT(m), m0 = NEXT(m0) ) {
! 190: NEXTDMM(mr0,mr); mr->c = m->c; mr->dl = m->dl; mr->pos = m->pos;
! 191: }
! 192: NEXT(mr) = 0; MKDPM(p0->nv,mr0,h); MKDPM(p0->nv,m,r);
! 193: }
! 194: if ( h )
! 195: h->sugar = p0->sugar;
! 196: if ( r )
! 197: r->sugar = p1->sugar;
! 198: *hp = h; *rp = r;
1.66 noro 199: }
200:
201:
1.39 ohara 202: void dp_ptozp3(DP p,Q *dvr,DP *rp)
203: {
1.67 ! noro 204: MP m,mr,mr0;
! 205: int i,n;
! 206: Q *w;
! 207: P t;
! 208:
! 209: if ( !p ) {
! 210: *rp = 0; *dvr = 0;
! 211: }else {
! 212: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
! 213: w = (Q *)ALLOCA(n*sizeof(Q));
! 214: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
! 215: if ( NUM(m->c) )
! 216: w[i] = (Q)m->c;
! 217: else
! 218: ptozp((P)m->c,1,&w[i],&t);
! 219: sortbynm(w,n);
! 220: qltozl(w,n,dvr);
! 221: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 222: NEXTMP(mr0,mr); divsp(CO,(P)m->c,(P)(*dvr),(P *)&mr->c); mr->dl = m->dl;
! 223: }
! 224: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 225: }
1.7 noro 226: }
1.1 noro 227:
1.20 noro 228: void dp_idiv(DP p,Q c,DP *rp)
1.1 noro 229: {
1.67 ! noro 230: Q t;
! 231: N nm,q;
! 232: int sgn,s;
! 233: MP mr0,m,mr;
! 234:
! 235: if ( !p )
! 236: *rp = 0;
! 237: else if ( MUNIQ((Q)c) )
! 238: *rp = p;
! 239: else if ( MUNIQ((Q)c) )
! 240: chsgnd(p,rp);
! 241: else {
! 242: nm = NM(c); sgn = SGN(c);
! 243: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 244: NEXTMP(mr0,mr);
! 245:
! 246: divsn(NM((Q)(m->c)),nm,&q);
! 247: s = sgn*SGN((Q)(m->c));
! 248: NTOQ(q,s,t);
! 249: mr->c = (Obj)t;
! 250: mr->dl = m->dl;
! 251: }
! 252: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
! 253: if ( *rp )
! 254: (*rp)->sugar = p->sugar;
! 255: }
1.1 noro 256: }
257:
1.20 noro 258: void dp_mbase(NODE hlist,NODE *mbase)
1.1 noro 259: {
1.67 ! noro 260: DL *dl;
! 261: DL d;
1.65 noro 262: int *t;
1.67 ! noro 263: int i,j,k,n,nvar,td;
1.1 noro 264:
1.67 ! noro 265: n = length(hlist); nvar = ((DP)BDY(hlist))->nv;
! 266: dl = (DL *)MALLOC(n*sizeof(DL));
! 267: NEWDL(d,nvar); *mbase = 0;
! 268: for ( i = 0; i < n; i++, hlist = NEXT(hlist) ) {
! 269: dl[i] = BDY((DP)BDY(hlist))->dl;
1.65 noro 270: /* trivial ideal check */
1.67 ! noro 271: if ( (*cmpdl)(nvar,d,dl[i]) == 0 ) {
1.65 noro 272: return;
273: }
274: }
275: /* zero-dim. ideal check */
276: for ( i = 0; i < nvar; i++ ) {
277: for ( j = 0; j < n; j++ ) {
278: for ( k = 0, t = dl[j]->d; k < nvar; k++ )
279: if ( k != i && t[k] != 0 ) break;
280: if ( k == nvar ) break;
281: }
282: if ( j == n )
283: error("dp_mbase : input ideal is not zero-dimensional");
284: }
1.67 ! noro 285: while ( 1 ) {
! 286: insert_to_node(d,mbase,nvar);
! 287: for ( i = nvar-1; i >= 0; ) {
! 288: d->d[i]++;
! 289: d->td += MUL_WEIGHT(1,i);
! 290: for ( j = 0; j < n; j++ ) {
! 291: if ( _dl_redble(dl[j],d,nvar) )
! 292: break;
! 293: }
! 294: if ( j < n ) {
! 295: for ( j = nvar-1; j >= i; j-- )
! 296: d->d[j] = 0;
! 297: for ( j = 0, td = 0; j < i; j++ )
! 298: td += MUL_WEIGHT(d->d[j],j);
! 299: d->td = td;
! 300: i--;
! 301: } else
! 302: break;
! 303: }
! 304: if ( i < 0 )
! 305: break;
! 306: }
1.1 noro 307: }
308:
1.20 noro 309: int _dl_redble(DL d1,DL d2,int nvar)
1.1 noro 310: {
1.67 ! noro 311: int i;
1.1 noro 312:
1.67 ! noro 313: if ( d1->td > d2->td )
! 314: return 0;
! 315: for ( i = 0; i < nvar; i++ )
! 316: if ( d1->d[i] > d2->d[i] )
! 317: break;
! 318: if ( i < nvar )
! 319: return 0;
! 320: else
! 321: return 1;
1.1 noro 322: }
323:
1.20 noro 324: void insert_to_node(DL d,NODE *n,int nvar)
1.1 noro 325: {
1.67 ! noro 326: DL d1;
! 327: MP m;
! 328: DP dp;
! 329: NODE n0,n1,n2;
! 330:
! 331: NEWDL(d1,nvar); d1->td = d->td;
! 332: bcopy((char *)d->d,(char *)d1->d,nvar*sizeof(int));
! 333: NEWMP(m); m->dl = d1; m->c = (Obj)ONE; NEXT(m) = 0;
! 334: MKDP(nvar,m,dp); dp->sugar = d->td;
! 335: if ( !(*n) ) {
! 336: MKNODE(n1,dp,0); *n = n1;
! 337: } else {
! 338: for ( n1 = *n, n0 = 0; n1; n0 = n1, n1 = NEXT(n1) )
! 339: if ( (*cmpdl)(nvar,d,BDY((DP)BDY(n1))->dl) > 0 ) {
! 340: MKNODE(n2,dp,n1);
! 341: if ( !n0 )
! 342: *n = n2;
! 343: else
! 344: NEXT(n0) = n2;
! 345: break;
! 346: }
! 347: if ( !n1 ) {
! 348: MKNODE(n2,dp,0); NEXT(n0) = n2;
! 349: }
! 350: }
1.1 noro 351: }
352:
1.20 noro 353: void dp_vtod(Q *c,DP p,DP *rp)
1.1 noro 354: {
1.67 ! noro 355: MP mr0,m,mr;
! 356: int i;
1.1 noro 357:
1.67 ! noro 358: if ( !p )
! 359: *rp = 0;
! 360: else {
! 361: for ( mr0 = 0, m = BDY(p), i = 0; m; m = NEXT(m), i++ ) {
! 362: NEXTMP(mr0,mr); mr->c = (Obj)c[i]; mr->dl = m->dl;
! 363: }
! 364: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
! 365: (*rp)->sugar = p->sugar;
! 366: }
1.1 noro 367: }
368:
1.8 noro 369: extern int mpi_mag;
370: extern int PCoeffs;
371:
1.20 noro 372: void dp_ptozp_d(DP p,DP *rp)
1.1 noro 373: {
1.67 ! noro 374: int i,j,k,l,n,nsep;
! 375: MP m;
! 376: NODE tn,n0,n1,n2,n3;
! 377: struct oVECT v;
! 378: VECT c,cs;
! 379: VECT qi,ri;
! 380: LIST *qr;
! 381: Obj dmy;
! 382: Q d0,d1,gcd,a,u,u1;
! 383: Q *q,*r;
! 384: STRING iqr_v;
! 385: pointer *b;
! 386: N qn,gn;
! 387: double get_rtime();
! 388: int blen;
! 389: NODE dist;
! 390: int ndist;
! 391: double t0;
! 392: double t_e,t_d,t_d1,t_c;
! 393: extern int DP_NFStat;
! 394: extern LIST Dist;
! 395: void Pox_rpc();
! 396: void Pox_pop_local();
! 397:
! 398: if ( !p )
! 399: *rp = 0;
! 400: else {
! 401: if ( PCoeffs ) {
! 402: dp_ptozp(p,rp); return;
! 403: }
! 404: if ( !Dist || p_mag((P)BDY(p)->c) <= mpi_mag ) {
! 405: dist = 0; ndist = 0;
! 406: if ( DP_NFStat ) fprintf(asir_out,"L");
! 407: } else {
! 408: dist = BDY(Dist); ndist = length(dist);
! 409: if ( DP_NFStat ) fprintf(asir_out,"D");
! 410: }
! 411: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
! 412: nsep = ndist + 1;
! 413: if ( n <= nsep ) {
! 414: dp_ptozp(p,rp); return;
! 415: }
! 416: t0 = get_rtime();
! 417: dp_dtov(p,&c);
! 418: igcdv_estimate(c,&d0);
! 419: t_e = get_rtime()-t0;
! 420: t0 = get_rtime();
! 421: dp_dtov(p,&c);
! 422: sepvect(c,nsep,&cs);
! 423: MKSTR(iqr_v,"iqr");
! 424: qr = (LIST *)CALLOC(nsep,sizeof(LIST));
! 425: q = (Q *)CALLOC(n,sizeof(Q));
! 426: r = (Q *)CALLOC(n,sizeof(Q));
! 427: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
! 428: MKNODE(n3,d0,0); MKNODE(n2,b[i],n3);
! 429: MKNODE(n1,iqr_v,n2); MKNODE(n0,BDY(tn),n1);
! 430: Pox_rpc(n0,&dmy);
! 431: }
! 432: iqrv(b[i],d0,&qr[i]);
! 433: dp_dtov(p,&c);
! 434: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
! 435: Pox_pop_local(tn,&qr[i]);
! 436: if ( OID(qr[i]) == O_ERR ) {
! 437: printexpr(CO,(Obj)qr[i]);
! 438: error("dp_ptozp_d : aborted");
! 439: }
! 440: }
! 441: t_d = get_rtime()-t0;
! 442: t_d1 = t_d/n;
! 443: t0 = get_rtime();
! 444: for ( i = j = 0; i < nsep; i++ ) {
! 445: tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
! 446: for ( k = 0, l = qi->len; k < l; k++, j++ ) {
! 447: q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
! 448: }
! 449: }
! 450: v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
! 451: if ( d1 ) {
! 452: gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
! 453: divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
! 454: for ( i = 0; i < n; i++ ) {
! 455: mulq(a,q[i],&u);
! 456: if ( r[i] ) {
! 457: divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
! 458: addq(u,u1,&q[i]);
! 459: } else
! 460: q[i] = u;
! 461: }
! 462: } else
! 463: gcd = d0;
! 464: dp_vtod(q,p,rp);
! 465: t_c = get_rtime()-t0;
! 466: blen=p_mag((P)gcd);
! 467: pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
! 468: if ( 0 )
! 469: fprintf(stderr,"(%d,%d)",p_mag((P)d0)-blen,blen);
! 470: }
1.1 noro 471: }
472:
1.20 noro 473: void dp_ptozp2_d(DP p0,DP p1,DP *hp,DP *rp)
1.1 noro 474: {
1.67 ! noro 475: DP t,s,h,r;
! 476: MP m,mr,mr0,m0;
1.1 noro 477:
1.67 ! noro 478: addd(CO,p0,p1,&t); dp_ptozp_d(t,&s);
! 479: if ( !p0 ) {
! 480: h = 0; r = s;
! 481: } else if ( !p1 ) {
! 482: h = s; r = 0;
! 483: } else {
! 484: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
! 485: m = NEXT(m), m0 = NEXT(m0) ) {
! 486: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
! 487: }
! 488: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
! 489: }
! 490: if ( h )
! 491: h->sugar = p0->sugar;
! 492: if ( r )
! 493: r->sugar = p1->sugar;
! 494: *hp = h; *rp = r;
1.5 noro 495: }
496:
1.22 noro 497: int have_sf_coef(P p)
498: {
1.67 ! noro 499: DCP dc;
1.22 noro 500:
1.67 ! noro 501: if ( !p )
! 502: return 0;
! 503: else if ( NUM(p) )
! 504: return NID((Num)p) == N_GFS ? 1 : 0;
! 505: else {
! 506: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 507: if ( have_sf_coef(COEF(dc)) )
! 508: return 1;
! 509: return 0;
! 510: }
1.22 noro 511: }
512:
1.25 noro 513: void head_coef(P p,Num *c)
514: {
1.67 ! noro 515: if ( !p )
! 516: *c = 0;
! 517: else if ( NUM(p) )
! 518: *c = (Num)p;
! 519: else
! 520: head_coef(COEF(DC(p)),c);
1.25 noro 521: }
522:
523: void dp_monic_sf(DP p,DP *rp)
524: {
1.67 ! noro 525: Num c;
1.25 noro 526:
1.67 ! noro 527: if ( !p )
! 528: *rp = 0;
! 529: else {
! 530: head_coef((P)BDY(p)->c,&c);
! 531: divsdc(CO,p,(P)c,rp);
! 532: }
1.25 noro 533: }
534:
1.20 noro 535: void dp_prim(DP p,DP *rp)
1.5 noro 536: {
1.67 ! noro 537: P t,g;
! 538: DP p1;
! 539: MP m,mr,mr0;
! 540: int i,n;
! 541: P *w;
! 542: Q *c;
! 543: Q dvr;
! 544: NODE tn;
! 545:
! 546: if ( !p )
! 547: *rp = 0;
! 548: else if ( dp_fcoeffs == N_GFS ) {
! 549: for ( m = BDY(p); m; m = NEXT(m) )
! 550: if ( OID(m->c) == O_N ) {
! 551: /* GCD of coeffs = 1 */
! 552: dp_monic_sf(p,rp);
! 553: return;
! 554: } else break;
! 555: /* compute GCD over the finite fieid */
! 556: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
! 557: w = (P *)ALLOCA(n*sizeof(P));
! 558: for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ )
! 559: w[i] = (P)m->c;
! 560: gcdsf(CO,w,n,&g);
! 561: if ( NUM(g) )
! 562: dp_monic_sf(p,rp);
! 563: else {
! 564: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 565: NEXTMP(mr0,mr); divsp(CO,(P)m->c,g,(P *)&mr->c); mr->dl = m->dl;
! 566: }
! 567: NEXT(mr) = 0; MKDP(p->nv,mr0,p1); p1->sugar = p->sugar;
! 568: dp_monic_sf(p1,rp);
! 569: }
! 570: return;
! 571: } else if ( dp_fcoeffs )
! 572: *rp = p;
! 573: else if ( NoGCD )
! 574: dp_ptozp(p,rp);
! 575: else {
! 576: dp_ptozp(p,&p1); p = p1;
! 577: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
! 578: if ( n == 1 ) {
! 579: m = BDY(p);
! 580: NEWMP(mr); mr->dl = m->dl; mr->c = (Obj)ONE; NEXT(mr) = 0;
! 581: MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
! 582: return;
! 583: }
! 584: w = (P *)ALLOCA(n*sizeof(P));
! 585: c = (Q *)ALLOCA(n*sizeof(Q));
! 586: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
! 587: if ( NUM(m->c) ) {
! 588: c[i] = (Q)m->c; w[i] = (P)ONE;
! 589: } else
! 590: ptozp((P)m->c,1,&c[i],&w[i]);
! 591: qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
! 592: if ( NUM(g) )
! 593: *rp = p;
! 594: else {
! 595: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 596: NEXTMP(mr0,mr); divsp(CO,(P)m->c,g,(P *)&mr->c); mr->dl = m->dl;
! 597: }
! 598: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 599: add_denomlist(g);
! 600: }
! 601: }
1.5 noro 602: }
603:
1.20 noro 604: void heu_nezgcdnpz(VL vl,P *pl,int m,P *pr)
1.5 noro 605: {
1.67 ! noro 606: int i,r;
! 607: P gcd,t,s1,s2,u;
! 608: Q rq;
! 609: DCP dc;
! 610: extern int DP_Print;
! 611:
! 612: while ( 1 ) {
! 613: for ( i = 0, s1 = 0; i < m; i++ ) {
! 614: r = random(); UTOQ(r,rq);
! 615: mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
! 616: }
! 617: for ( i = 0, s2 = 0; i < m; i++ ) {
! 618: r = random(); UTOQ(r,rq);
! 619: mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
! 620: }
! 621: ezgcdp(vl,s1,s2,&gcd);
! 622: if ( DP_Print > 2 )
! 623: { fprintf(asir_out,"(%d)",nmonop(gcd)); fflush(asir_out); }
! 624: for ( i = 0; i < m; i++ ) {
! 625: if ( !divtpz(vl,pl[i],gcd,&t) )
! 626: break;
! 627: }
! 628: if ( i == m )
! 629: break;
! 630: }
! 631: *pr = gcd;
1.5 noro 632: }
633:
1.20 noro 634: void dp_prim_mod(DP p,int mod,DP *rp)
1.5 noro 635: {
1.67 ! noro 636: P t,g;
! 637: MP m,mr,mr0;
1.5 noro 638:
1.67 ! noro 639: if ( !p )
! 640: *rp = 0;
! 641: else if ( NoGCD )
! 642: *rp = p;
! 643: else {
! 644: for ( m = BDY(p), g = (P)m->c, m = NEXT(m); m; m = NEXT(m) ) {
! 645: gcdprsmp(CO,mod,g,(P)m->c,&t); g = t;
! 646: }
! 647: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 648: NEXTMP(mr0,mr); divsmp(CO,mod,(P)m->c,g,(P *)&mr->c); mr->dl = m->dl;
! 649: }
! 650: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 651: }
1.5 noro 652: }
653:
1.20 noro 654: void dp_cont(DP p,Q *rp)
1.5 noro 655: {
1.67 ! noro 656: VECT v;
1.5 noro 657:
1.67 ! noro 658: dp_dtov(p,&v); igcdv(v,rp);
1.5 noro 659: }
660:
1.20 noro 661: void dp_dtov(DP dp,VECT *rp)
1.5 noro 662: {
1.67 ! noro 663: MP m,t;
! 664: int i,n;
! 665: VECT v;
! 666: pointer *p;
! 667:
! 668: m = BDY(dp);
! 669: for ( t = m, n = 0; t; t = NEXT(t), n++ );
! 670: MKVECT(v,n);
! 671: for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
! 672: p[i] = (pointer)(t->c);
! 673: *rp = v;
1.5 noro 674: }
675:
1.7 noro 676: /*
677: * s-poly computation
678: *
679: */
1.5 noro 680:
1.20 noro 681: void dp_sp(DP p1,DP p2,DP *rp)
1.5 noro 682: {
1.67 ! noro 683: int i,n,td;
! 684: int *w;
! 685: DL d1,d2,d;
! 686: MP m;
! 687: DP t,s1,s2,u;
! 688: Q c,c1,c2;
! 689: N gn,tn;
! 690:
! 691: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 692: w = (int *)ALLOCA(n*sizeof(int));
! 693: for ( i = 0, td = 0; i < n; i++ ) {
! 694: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
! 695: }
! 696:
! 697: NEWDL(d,n); d->td = td - d1->td;
! 698: for ( i = 0; i < n; i++ )
! 699: d->d[i] = w[i] - d1->d[i];
! 700: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 701: if ( INT(c1) && INT(c2) ) {
! 702: gcdn(NM(c1),NM(c2),&gn);
! 703: if ( !UNIN(gn) ) {
! 704: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 705: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 706: }
! 707: }
! 708:
! 709: NEWMP(m); m->dl = d; m->c = (Obj)c2; NEXT(m) = 0;
! 710: MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
! 711:
! 712: NEWDL(d,n); d->td = td - d2->td;
! 713: for ( i = 0; i < n; i++ )
! 714: d->d[i] = w[i] - d2->d[i];
! 715: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0;
! 716: MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
! 717:
! 718: subd(CO,t,u,rp);
! 719: if ( GenTrace ) {
! 720: LIST hist;
! 721: NODE node;
! 722:
! 723: node = mknode(4,ONE,NULLP,s1,ONE);
! 724: MKLIST(hist,node);
! 725: MKNODE(TraceList,hist,0);
! 726:
! 727: node = mknode(4,ONE,NULLP,NULLP,ONE);
! 728: chsgnd(s2,(DP *)&ARG2(node));
! 729: MKLIST(hist,node);
! 730: MKNODE(node,hist,TraceList); TraceList = node;
! 731: }
1.14 noro 732: }
733:
1.66 noro 734: void dpm_sp(DPM p1,DPM p2,DPM *rp)
735: {
1.67 ! noro 736: int i,n,td;
! 737: int *w;
! 738: DL d1,d2,d;
! 739: MP m;
! 740: DP s1,s2;
1.66 noro 741: DPM t,u;
1.67 ! noro 742: Q c,c1,c2;
! 743: N gn,tn;
1.66 noro 744:
1.67 ! noro 745: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.66 noro 746: if ( BDY(p1)->pos != BDY(p2)->pos ) {
747: *rp = 0;
748: return;
749: }
1.67 ! noro 750: w = (int *)ALLOCA(n*sizeof(int));
! 751: for ( i = 0, td = 0; i < n; i++ ) {
! 752: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
! 753: }
! 754:
! 755: NEWDL(d,n); d->td = td - d1->td;
! 756: for ( i = 0; i < n; i++ )
! 757: d->d[i] = w[i] - d1->d[i];
! 758: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 759: if ( INT(c1) && INT(c2) ) {
! 760: gcdn(NM(c1),NM(c2),&gn);
! 761: if ( !UNIN(gn) ) {
! 762: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 763: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 764: }
! 765: }
! 766:
! 767: NEWMP(m); m->dl = d; m->c = (Obj)c2; NEXT(m) = 0;
! 768: MKDP(n,m,s1); s1->sugar = d->td; mulobjdpm(CO,(Obj)s1,p1,&t);
! 769:
! 770: NEWDL(d,n); d->td = td - d2->td;
! 771: for ( i = 0; i < n; i++ )
! 772: d->d[i] = w[i] - d2->d[i];
! 773: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0;
! 774: MKDP(n,m,s2); s2->sugar = d->td; mulobjdpm(CO,(Obj)s2,p2,&u);
! 775:
! 776: subdpm(CO,t,u,rp);
! 777: if ( GenTrace ) {
! 778: LIST hist;
! 779: NODE node;
! 780:
! 781: node = mknode(4,ONE,NULLP,s1,ONE);
! 782: MKLIST(hist,node);
! 783: MKNODE(TraceList,hist,0);
! 784:
! 785: node = mknode(4,ONE,NULLP,NULLP,ONE);
! 786: chsgnd(s2,(DP *)&ARG2(node));
! 787: MKLIST(hist,node);
! 788: MKNODE(node,hist,TraceList); TraceList = node;
! 789: }
1.66 noro 790: }
791:
1.20 noro 792: void _dp_sp_dup(DP p1,DP p2,DP *rp)
1.14 noro 793: {
1.67 ! noro 794: int i,n,td;
! 795: int *w;
! 796: DL d1,d2,d;
! 797: MP m;
! 798: DP t,s1,s2,u;
! 799: Q c,c1,c2;
! 800: N gn,tn;
! 801:
! 802: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 803: w = (int *)ALLOCA(n*sizeof(int));
! 804: for ( i = 0, td = 0; i < n; i++ ) {
! 805: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
! 806: }
! 807:
! 808: _NEWDL(d,n); d->td = td - d1->td;
! 809: for ( i = 0; i < n; i++ )
! 810: d->d[i] = w[i] - d1->d[i];
! 811: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 812: if ( INT(c1) && INT(c2) ) {
! 813: gcdn(NM(c1),NM(c2),&gn);
! 814: if ( !UNIN(gn) ) {
! 815: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 816: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 817: }
! 818: }
! 819:
! 820: _NEWMP(m); m->dl = d; m->c = (Obj)c2; NEXT(m) = 0;
! 821: _MKDP(n,m,s1); s1->sugar = d->td; _muld_dup(CO,s1,p1,&t); _free_dp(s1);
! 822:
! 823: _NEWDL(d,n); d->td = td - d2->td;
! 824: for ( i = 0; i < n; i++ )
! 825: d->d[i] = w[i] - d2->d[i];
! 826: _NEWMP(m); m->dl = d; chsgnp((P)c1,(P *)&m->c); NEXT(m) = 0;
! 827: _MKDP(n,m,s2); s2->sugar = d->td; _muld_dup(CO,s2,p2,&u); _free_dp(s2);
! 828:
! 829: _addd_destructive(CO,t,u,rp);
! 830: if ( GenTrace ) {
! 831: LIST hist;
! 832: NODE node;
! 833:
! 834: node = mknode(4,ONE,NULLP,s1,ONE);
! 835: MKLIST(hist,node);
! 836: MKNODE(TraceList,hist,0);
! 837:
! 838: node = mknode(4,ONE,NULLP,NULLP,ONE);
! 839: chsgnd(s2,(DP *)&ARG2(node));
! 840: MKLIST(hist,node);
! 841: MKNODE(node,hist,TraceList); TraceList = node;
! 842: }
1.7 noro 843: }
844:
1.20 noro 845: void dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7 noro 846: {
1.67 ! noro 847: int i,n,td;
! 848: int *w;
! 849: DL d1,d2,d;
! 850: MP m;
! 851: DP t,s,u;
! 852:
! 853: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 854: w = (int *)ALLOCA(n*sizeof(int));
! 855: for ( i = 0, td = 0; i < n; i++ ) {
! 856: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
! 857: }
! 858: NEWDL_NOINIT(d,n); d->td = td - d1->td;
! 859: for ( i = 0; i < n; i++ )
! 860: d->d[i] = w[i] - d1->d[i];
! 861: NEWMP(m); m->dl = d; m->c = (Obj)BDY(p2)->c; NEXT(m) = 0;
! 862: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
! 863: NEWDL_NOINIT(d,n); d->td = td - d2->td;
! 864: for ( i = 0; i < n; i++ )
! 865: d->d[i] = w[i] - d2->d[i];
! 866: NEWMP(m); m->dl = d; m->c = (Obj)BDY(p1)->c; NEXT(m) = 0;
! 867: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
! 868: submd(CO,mod,t,u,rp);
1.7 noro 869: }
870:
1.20 noro 871: void _dp_sp_mod_dup(DP p1,DP p2,int mod,DP *rp)
1.7 noro 872: {
1.67 ! noro 873: int i,n,td;
! 874: int *w;
! 875: DL d1,d2,d;
! 876: MP m;
! 877: DP t,s,u;
! 878:
! 879: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 880: w = (int *)ALLOCA(n*sizeof(int));
! 881: for ( i = 0, td = 0; i < n; i++ ) {
! 882: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
! 883: }
! 884: _NEWDL(d,n); d->td = td - d1->td;
! 885: for ( i = 0; i < n; i++ )
! 886: d->d[i] = w[i] - d1->d[i];
! 887: _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
! 888: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
! 889: _NEWDL(d,n); d->td = td - d2->td;
! 890: for ( i = 0; i < n; i++ )
! 891: d->d[i] = w[i] - d2->d[i];
! 892: _NEWMP(m); m->dl = d; m->c = (Obj)STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
! 893: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
! 894: _addmd_destructive(mod,t,u,rp);
1.7 noro 895: }
896:
1.20 noro 897: void _dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7 noro 898: {
1.67 ! noro 899: int i,n,td;
! 900: int *w;
! 901: DL d1,d2,d;
! 902: MP m;
! 903: DP t,s,u;
! 904:
! 905: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 906: w = (int *)ALLOCA(n*sizeof(int));
! 907: for ( i = 0, td = 0; i < n; i++ ) {
! 908: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
! 909: }
! 910: NEWDL(d,n); d->td = td - d1->td;
! 911: for ( i = 0; i < n; i++ )
! 912: d->d[i] = w[i] - d1->d[i];
! 913: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
! 914: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
! 915: NEWDL(d,n); d->td = td - d2->td;
! 916: for ( i = 0; i < n; i++ )
! 917: d->d[i] = w[i] - d2->d[i];
! 918: NEWMP(m); m->dl = d; m->c = (Obj)STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
! 919: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
! 920: addmd_destructive(mod,t,u,rp);
1.7 noro 921: }
922:
923: /*
924: * m-reduction
1.13 noro 925: * do content reduction over Z or Q(x,...)
926: * do nothing over finite fields
1.7 noro 927: *
928: */
929:
1.20 noro 930: void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp,DP *multp)
1.7 noro 931: {
1.67 ! noro 932: int i,n;
! 933: DL d1,d2,d;
! 934: MP m;
! 935: DP t,s,r,h;
! 936: Q c,c1,c2;
! 937: N gn,tn;
! 938: P g,a;
! 939: P p[2];
! 940:
! 941: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 942: NEWDL(d,n); d->td = d1->td - d2->td;
! 943: for ( i = 0; i < n; i++ )
! 944: d->d[i] = d1->d[i]-d2->d[i];
! 945: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 946: if ( dp_fcoeffs == N_GFS ) {
! 947: p[0] = (P)c1; p[1] = (P)c2;
! 948: gcdsf(CO,p,2,&g);
! 949: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 950: } else if ( dp_fcoeffs ) {
! 951: /* do nothing */
! 952: } else if ( INT(c1) && INT(c2) ) {
! 953: gcdn(NM(c1),NM(c2),&gn);
! 954: if ( !UNIN(gn) ) {
! 955: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 956: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 957: }
! 958: } else {
! 959: ezgcdpz(CO,(P)c1,(P)c2,&g);
! 960: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 961: add_denomlist(g);
! 962: }
! 963: NEWMP(m); m->dl = d; chsgnp((P)c1,(P *)&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
! 964: *multp = s;
! 965: muld(CO,s,p2,&t); muldc(CO,p1,(Obj)c2,&s); addd(CO,s,t,&r);
! 966: muldc(CO,p0,(Obj)c2,&h);
! 967: *head = h; *rest = r; *dnp = (P)c2;
1.7 noro 968: }
969:
1.66 noro 970: void dpm_red(DPM p0,DPM p1,DPM p2,DPM *head,DPM *rest,P *dnp,DP *multp)
971: {
1.67 ! noro 972: int i,n,pos;
! 973: DL d1,d2,d;
! 974: MP m;
1.66 noro 975: DP s;
1.67 ! noro 976: DPM t,r,h,u,w;
! 977: Q c,c1,c2;
! 978: N gn,tn;
! 979: P g,a;
! 980: P p[2];
1.66 noro 981:
1.67 ! noro 982: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; pos = BDY(p1)->pos;
1.66 noro 983: if ( pos != BDY(p2)->pos )
984: error("dpm_red : cannot happen");
1.67 ! noro 985: NEWDL(d,n); d->td = d1->td - d2->td;
! 986: for ( i = 0; i < n; i++ )
! 987: d->d[i] = d1->d[i]-d2->d[i];
! 988: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 989: if ( dp_fcoeffs == N_GFS ) {
! 990: p[0] = (P)c1; p[1] = (P)c2;
! 991: gcdsf(CO,p,2,&g);
! 992: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 993: } else if ( dp_fcoeffs ) {
! 994: /* do nothing */
! 995: } else if ( INT(c1) && INT(c2) ) {
! 996: gcdn(NM(c1),NM(c2),&gn);
! 997: if ( !UNIN(gn) ) {
! 998: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 999: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 1000: }
! 1001: } else {
! 1002: ezgcdpz(CO,(P)c1,(P)c2,&g);
! 1003: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 1004: add_denomlist(g);
! 1005: }
! 1006: NEWMP(m); m->dl = d; chsgnp((P)c1,(P *)&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
! 1007: *multp = s;
! 1008: mulobjdpm(CO,(Obj)s,p2,&u); mulobjdpm(CO,(Obj)c2,p1,&w); adddpm(CO,u,w,&r);
! 1009: mulobjdpm(CO,(Obj)c2,p0,&h);
! 1010: *head = h; *rest = r; *dnp = (P)c2;
1.66 noro 1011: }
1012:
1013:
1.41 noro 1014: /*
1015: * m-reduction by a marked poly
1016: * do content reduction over Z or Q(x,...)
1017: * do nothing over finite fields
1018: *
1019: */
1020:
1021:
1022: void dp_red_marked(DP p0,DP p1,DP p2,DP hp2,DP *head,DP *rest,P *dnp,DP *multp)
1023: {
1.67 ! noro 1024: int i,n;
! 1025: DL d1,d2,d;
! 1026: MP m;
! 1027: DP t,s,r,h;
! 1028: Q c,c1,c2;
! 1029: N gn,tn;
! 1030: P g,a;
! 1031: P p[2];
! 1032:
! 1033: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl;
! 1034: NEWDL(d,n); d->td = d1->td - d2->td;
! 1035: for ( i = 0; i < n; i++ )
! 1036: d->d[i] = d1->d[i]-d2->d[i];
! 1037: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(hp2)->c;
! 1038: if ( dp_fcoeffs == N_GFS ) {
! 1039: p[0] = (P)c1; p[1] = (P)c2;
! 1040: gcdsf(CO,p,2,&g);
! 1041: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 1042: } else if ( dp_fcoeffs ) {
! 1043: /* do nothing */
! 1044: } else if ( INT(c1) && INT(c2) ) {
! 1045: gcdn(NM(c1),NM(c2),&gn);
! 1046: if ( !UNIN(gn) ) {
! 1047: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 1048: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 1049: }
! 1050: } else {
! 1051: ezgcdpz(CO,(P)c1,(P)c2,&g);
! 1052: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 1053: }
! 1054: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
! 1055: *multp = s;
! 1056: muld(CO,s,p2,&t); muldc(CO,p1,(Obj)c2,&s); subd(CO,s,t,&r);
! 1057: muldc(CO,p0,(Obj)c2,&h);
! 1058: *head = h; *rest = r; *dnp = (P)c2;
1.41 noro 1059: }
1060:
1.55 noro 1061: void dp_red_marked_mod(DP p0,DP p1,DP p2,DP hp2,int mod,DP *head,DP *rest,P *dnp,DP *multp)
1.44 noro 1062: {
1.67 ! noro 1063: int i,n;
! 1064: DL d1,d2,d;
! 1065: MP m;
! 1066: DP t,s,r,h;
! 1067: P c1,c2,g,u;
! 1068:
! 1069: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl;
! 1070: NEWDL(d,n); d->td = d1->td - d2->td;
! 1071: for ( i = 0; i < n; i++ )
! 1072: d->d[i] = d1->d[i]-d2->d[i];
! 1073: c1 = (P)BDY(p1)->c; c2 = (P)BDY(hp2)->c;
! 1074: gcdprsmp(CO,mod,c1,c2,&g);
! 1075: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
! 1076: if ( NUM(c2) ) {
! 1077: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
! 1078: }
! 1079: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0;
! 1080: MKDP(n,m,s); s->sugar = d->td;
! 1081: *multp = s;
! 1082: mulmd(CO,mod,s,p2,&t);
! 1083: if ( NUM(c2) ) {
! 1084: submd(CO,mod,p1,t,&r); h = p0;
! 1085: } else {
! 1086: mulmdc(CO,mod,p1,c2,&s); submd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
! 1087: }
! 1088: *head = h; *rest = r; *dnp = c2;
1.44 noro 1089: }
1090:
1.13 noro 1091: /* m-reduction over a field */
1092:
1.20 noro 1093: void dp_red_f(DP p1,DP p2,DP *rest)
1.13 noro 1094: {
1.67 ! noro 1095: int i,n;
! 1096: DL d1,d2,d;
! 1097: MP m;
! 1098: DP t,s;
! 1099: Obj a,b;
! 1100:
! 1101: n = p1->nv;
! 1102: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 1103:
! 1104: NEWDL(d,n); d->td = d1->td - d2->td;
! 1105: for ( i = 0; i < n; i++ )
! 1106: d->d[i] = d1->d[i]-d2->d[i];
! 1107:
! 1108: NEWMP(m); m->dl = d;
! 1109: divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
! 1110: C(m) = (Obj)b;
! 1111: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
1.13 noro 1112:
1.67 ! noro 1113: muld(CO,s,p2,&t); addd(CO,p1,t,rest);
1.13 noro 1114: }
1115:
1.66 noro 1116: void dpm_red_f(DPM p1,DPM p2,DPM *rest)
1117: {
1.67 ! noro 1118: int i,n;
! 1119: DL d1,d2,d;
! 1120: MP m;
! 1121: DPM t;
1.66 noro 1122: DP s;
1.67 ! noro 1123: Obj a,b;
1.66 noro 1124:
1.67 ! noro 1125: n = p1->nv;
! 1126: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.66 noro 1127:
1.67 ! noro 1128: NEWDL(d,n); d->td = d1->td - d2->td;
! 1129: for ( i = 0; i < n; i++ )
! 1130: d->d[i] = d1->d[i]-d2->d[i];
1.66 noro 1131:
1.67 ! noro 1132: NEWMP(m); m->dl = d;
! 1133: arf_div(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); arf_chsgn(a,&b);
! 1134: C(m) = b;
! 1135: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
1.66 noro 1136:
1.67 ! noro 1137: mulobjdpm(CO,(Obj)s,p2,&t); adddpm(CO,p1,t,rest);
1.66 noro 1138: }
1139:
1140:
1.20 noro 1141: void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp)
1.7 noro 1142: {
1.67 ! noro 1143: int i,n;
! 1144: DL d1,d2,d;
! 1145: MP m;
! 1146: DP t,s,r,h;
! 1147: P c1,c2,g,u;
! 1148:
! 1149: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 1150: NEWDL(d,n); d->td = d1->td - d2->td;
! 1151: for ( i = 0; i < n; i++ )
! 1152: d->d[i] = d1->d[i]-d2->d[i];
! 1153: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
! 1154: gcdprsmp(CO,mod,c1,c2,&g);
! 1155: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
! 1156: if ( NUM(c2) ) {
! 1157: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
! 1158: }
! 1159: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,(P *)&m->c); NEXT(m) = 0;
! 1160: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
! 1161: if ( NUM(c2) ) {
! 1162: addmd(CO,mod,p1,t,&r); h = p0;
! 1163: } else {
! 1164: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
! 1165: }
! 1166: *head = h; *rest = r; *dnp = c2;
1.7 noro 1167: }
1168:
1.10 noro 1169: struct oEGT eg_red_mod;
1170:
1.20 noro 1171: void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp)
1.7 noro 1172: {
1.67 ! noro 1173: int i,n;
! 1174: DL d1,d2,d;
! 1175: MP m;
! 1176: DP t,s;
! 1177: int c,c1,c2;
! 1178: extern int do_weyl;
! 1179:
! 1180: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 1181: _NEWDL(d,n); d->td = d1->td - d2->td;
! 1182: for ( i = 0; i < n; i++ )
! 1183: d->d[i] = d1->d[i]-d2->d[i];
! 1184: c = invm(ITOS(BDY(p2)->c),mod);
! 1185: c2 = ITOS(BDY(p1)->c);
! 1186: DMAR(c,c2,0,mod,c1);
! 1187: _NEWMP(m); m->dl = d; m->c = (Obj)STOI(mod-c1); NEXT(m) = 0;
1.16 noro 1188: #if 0
1.67 ! noro 1189: _MKDP(n,m,s); s->sugar = d->td;
! 1190: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 1191: #else
1.67 ! noro 1192: if ( do_weyl ) {
! 1193: _MKDP(n,m,s); s->sugar = d->td;
! 1194: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
! 1195: } else {
! 1196: _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
! 1197: }
1.16 noro 1198: #endif
1.10 noro 1199: /* get_eg(&t0); */
1.67 ! noro 1200: _addmd_destructive(mod,p1,t,rp);
1.10 noro 1201: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 1202: }
1203:
1204: /*
1205: * normal form computation
1206: *
1207: */
1.5 noro 1208:
1.20 noro 1209: void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp)
1.5 noro 1210: {
1.67 ! noro 1211: DP u,p,d,s,t,dmy;
! 1212: NODE l;
! 1213: MP m,mr;
! 1214: int i,n;
! 1215: int *wb;
! 1216: int sugar,psugar;
! 1217: P dn,tdn,tdn1;
! 1218:
! 1219: dn = (P)ONE;
! 1220: if ( !g ) {
! 1221: *rp = 0; *dnp = dn; return;
! 1222: }
! 1223: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1224: wb = (int *)ALLOCA(n*sizeof(int));
! 1225: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1226: wb[i] = QTOS((Q)BDY(l));
! 1227: sugar = g->sugar;
! 1228: for ( d = 0; g; ) {
! 1229: for ( u = 0, i = 0; i < n; i++ ) {
! 1230: if ( dp_redble(g,p = ps[wb[i]]) ) {
! 1231: dp_red(d,g,p,&t,&u,&tdn,&dmy);
! 1232: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1233: sugar = MAX(sugar,psugar);
! 1234: if ( !u ) {
! 1235: if ( d )
! 1236: d->sugar = sugar;
! 1237: *rp = d; *dnp = dn; return;
! 1238: } else {
! 1239: d = t;
! 1240: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
! 1241: }
! 1242: break;
! 1243: }
! 1244: }
! 1245: if ( u )
! 1246: g = u;
! 1247: else if ( !full ) {
! 1248: if ( g ) {
! 1249: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1250: }
! 1251: *rp = g; *dnp = dn; return;
! 1252: } else {
! 1253: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1254: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1255: addd(CO,d,t,&s); d = s;
! 1256: dp_rest(g,&t); g = t;
! 1257: }
! 1258: }
! 1259: if ( d )
! 1260: d->sugar = sugar;
! 1261: *rp = d; *dnp = dn;
1.5 noro 1262: }
1263:
1.43 noro 1264: void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Q *contp)
1265: {
1.67 ! noro 1266: struct oVECT v;
! 1267: int i,n1,n2,n;
! 1268: MP m,m0,t;
! 1269: Q *w;
! 1270: Q h;
! 1271:
! 1272: if ( p1 ) {
! 1273: for ( i = 0, m = BDY(p1); m; m = NEXT(m), i++ );
! 1274: n1 = i;
! 1275: } else
! 1276: n1 = 0;
! 1277: if ( p2 ) {
! 1278: for ( i = 0, m = BDY(p2); m; m = NEXT(m), i++ );
! 1279: n2 = i;
! 1280: } else
! 1281: n2 = 0;
! 1282: n = n1+n2;
! 1283: if ( !n ) {
! 1284: *r1p = 0; *r2p = 0; *contp = ONE; return;
! 1285: }
! 1286: w = (Q *)ALLOCA(n*sizeof(Q));
! 1287: v.len = n;
! 1288: v.body = (pointer *)w;
! 1289: i = 0;
! 1290: if ( p1 )
! 1291: for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) w[i] = (Q)m->c;
! 1292: if ( p2 )
! 1293: for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = (Q)m->c;
! 1294: h = w[0]; removecont_array((P *)w,n,1); divq(h,w[0],contp);
! 1295: i = 0;
! 1296: if ( p1 ) {
! 1297: for ( m0 = 0, t = BDY(p1); i < n1; i++, t = NEXT(t) ) {
! 1298: NEXTMP(m0,m); m->c = (Obj)w[i]; m->dl = t->dl;
! 1299: }
! 1300: NEXT(m) = 0;
! 1301: MKDP(p1->nv,m0,*r1p); (*r1p)->sugar = p1->sugar;
! 1302: } else
! 1303: *r1p = 0;
! 1304: if ( p2 ) {
! 1305: for ( m0 = 0, t = BDY(p2); i < n; i++, t = NEXT(t) ) {
! 1306: NEXTMP(m0,m); m->c = (Obj)w[i]; m->dl = t->dl;
! 1307: }
! 1308: NEXT(m) = 0;
! 1309: MKDP(p2->nv,m0,*r2p); (*r2p)->sugar = p2->sugar;
! 1310: } else
! 1311: *r2p = 0;
1.43 noro 1312: }
1313:
1.41 noro 1314: /* true nf by a marked GB */
1315:
1.43 noro 1316: void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP *rp,P *nmp,P *dnp)
1.41 noro 1317: {
1.67 ! noro 1318: DP u,p,d,s,t,dmy,hp;
! 1319: NODE l;
! 1320: MP m,mr;
! 1321: int i,n,hmag;
! 1322: int *wb;
! 1323: int sugar,psugar,multiple;
! 1324: P nm,tnm1,dn,tdn,tdn1;
! 1325: Q cont;
! 1326:
! 1327: multiple = 0;
! 1328: hmag = multiple*HMAG(g);
! 1329: nm = (P)ONE;
! 1330: dn = (P)ONE;
! 1331: if ( !g ) {
! 1332: *rp = 0; *dnp = dn; return;
! 1333: }
! 1334: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1335: wb = (int *)ALLOCA(n*sizeof(int));
! 1336: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1337: wb[i] = QTOS((Q)BDY(l));
! 1338: sugar = g->sugar;
! 1339: for ( d = 0; g; ) {
! 1340: for ( u = 0, i = 0; i < n; i++ ) {
! 1341: if ( dp_redble(g,hp = hps[wb[i]]) ) {
! 1342: p = ps[wb[i]];
! 1343: dp_red_marked(d,g,p,hp,&t,&u,&tdn,&dmy);
! 1344: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1345: sugar = MAX(sugar,psugar);
! 1346: if ( !u ) {
! 1347: goto last;
! 1348: } else {
! 1349: d = t;
! 1350: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
! 1351: }
! 1352: break;
! 1353: }
! 1354: }
! 1355: if ( u ) {
! 1356: g = u;
! 1357: if ( multiple && ((d && HMAG(d)>hmag) || (HMAG(g)>hmag)) ) {
! 1358: dp_removecont2(d,g,&t,&u,&cont); d = t; g = u;
! 1359: mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
! 1360: if ( d )
! 1361: hmag = multiple*HMAG(d);
! 1362: else
! 1363: hmag = multiple*HMAG(g);
! 1364: }
! 1365: } else {
! 1366: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1367: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1368: addd(CO,d,t,&s); d = s;
! 1369: dp_rest(g,&t); g = t;
! 1370: }
! 1371: }
1.43 noro 1372: last:
1.67 ! noro 1373: if ( d ) {
! 1374: dp_removecont2(d,0,&t,&u,&cont); d = t;
! 1375: mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
! 1376: d->sugar = sugar;
! 1377: }
! 1378: *rp = d; *nmp = nm; *dnp = dn;
1.41 noro 1379: }
1380:
1.44 noro 1381: void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
1382: {
1.67 ! noro 1383: DP hp,u,p,d,s,t,dmy;
! 1384: NODE l;
! 1385: MP m,mr;
! 1386: int i,n;
! 1387: int *wb;
! 1388: int sugar,psugar;
! 1389: P dn,tdn,tdn1;
! 1390:
! 1391: dn = (P)ONEM;
! 1392: if ( !g ) {
! 1393: *rp = 0; *dnp = dn; return;
! 1394: }
! 1395: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1396: wb = (int *)ALLOCA(n*sizeof(int));
! 1397: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1398: wb[i] = QTOS((Q)BDY(l));
! 1399: sugar = g->sugar;
! 1400: for ( d = 0; g; ) {
! 1401: for ( u = 0, i = 0; i < n; i++ ) {
! 1402: if ( dp_redble(g,hp = hps[wb[i]]) ) {
! 1403: p = ps[wb[i]];
! 1404: dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&dmy);
! 1405: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1406: sugar = MAX(sugar,psugar);
! 1407: if ( !u ) {
! 1408: if ( d )
! 1409: d->sugar = sugar;
! 1410: *rp = d; *dnp = dn; return;
! 1411: } else {
! 1412: d = t;
! 1413: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
! 1414: }
! 1415: break;
! 1416: }
! 1417: }
! 1418: if ( u )
! 1419: g = u;
! 1420: else {
! 1421: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1422: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1423: addmd(CO,mod,d,t,&s); d = s;
! 1424: dp_rest(g,&t); g = t;
! 1425: }
! 1426: }
! 1427: if ( d )
! 1428: d->sugar = sugar;
! 1429: *rp = d; *dnp = dn;
1.44 noro 1430: }
1431:
1.47 noro 1432: /* true nf by a marked GB and collect quotients */
1433:
1434: DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps,DP *hps,DP *rp,P *dnp)
1435: {
1.67 ! noro 1436: DP u,p,d,s,t,dmy,hp,mult;
! 1437: DP *q;
! 1438: NODE l;
! 1439: MP m,mr;
! 1440: int i,n,j;
! 1441: int *wb;
! 1442: int sugar,psugar,multiple;
! 1443: P nm,tnm1,dn,tdn,tdn1;
! 1444: Q cont;
! 1445:
! 1446: dn = (P)ONE;
! 1447: if ( !g ) {
! 1448: *rp = 0; *dnp = dn; return 0;
! 1449: }
! 1450: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1451: wb = (int *)ALLOCA(n*sizeof(int));
! 1452: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1453: wb[i] = QTOS((Q)BDY(l));
! 1454: q = (DP *)MALLOC(n*sizeof(DP));
! 1455: for ( i = 0; i < n; i++ ) q[i] = 0;
! 1456: sugar = g->sugar;
! 1457: for ( d = 0; g; ) {
! 1458: for ( u = 0, i = 0; i < n; i++ ) {
! 1459: if ( dp_redble(g,hp = hps[wb[i]]) ) {
! 1460: p = ps[wb[i]];
! 1461: dp_red_marked(d,g,p,hp,&t,&u,&tdn,&mult);
! 1462: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1463: sugar = MAX(sugar,psugar);
! 1464: for ( j = 0; j < n; j++ ) {
! 1465: muldc(CO,q[j],(Obj)tdn,&dmy); q[j] = dmy;
! 1466: }
! 1467: addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
! 1468: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
! 1469: d = t;
! 1470: if ( !u ) goto last;
! 1471: break;
! 1472: }
! 1473: }
! 1474: if ( u ) {
! 1475: g = u;
! 1476: } else {
! 1477: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1478: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1479: addd(CO,d,t,&s); d = s;
! 1480: dp_rest(g,&t); g = t;
! 1481: }
! 1482: }
1.47 noro 1483: last:
1.67 ! noro 1484: if ( d ) d->sugar = sugar;
! 1485: *rp = d; *dnp = dn;
! 1486: return q;
1.47 noro 1487: }
1488:
1.55 noro 1489: DP *dp_true_nf_and_quotient_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
1490: {
1.67 ! noro 1491: DP u,p,d,s,t,dmy,hp,mult;
! 1492: DP *q;
! 1493: NODE l;
! 1494: MP m,mr;
! 1495: int i,n,j;
! 1496: int *wb;
! 1497: int sugar,psugar;
! 1498: P dn,tdn,tdn1;
! 1499:
! 1500: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1501: q = (DP *)MALLOC(n*sizeof(DP));
! 1502: for ( i = 0; i < n; i++ ) q[i] = 0;
! 1503: dn = (P)ONEM;
! 1504: if ( !g ) {
! 1505: *rp = 0; *dnp = dn; return 0;
! 1506: }
! 1507: wb = (int *)ALLOCA(n*sizeof(int));
! 1508: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1509: wb[i] = QTOS((Q)BDY(l));
! 1510: sugar = g->sugar;
! 1511: for ( d = 0; g; ) {
! 1512: for ( u = 0, i = 0; i < n; i++ ) {
! 1513: if ( dp_redble(g,hp = hps[wb[i]]) ) {
! 1514: p = ps[wb[i]];
! 1515: dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&mult);
! 1516: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1517: sugar = MAX(sugar,psugar);
! 1518: for ( j = 0; j < n; j++ ) {
! 1519: mulmdc(CO,mod,q[j],(P)tdn,&dmy); q[j] = dmy;
! 1520: }
! 1521: addmd(CO,mod,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
! 1522: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
! 1523: d = t;
! 1524: if ( !u ) goto last;
! 1525: break;
! 1526: }
! 1527: }
! 1528: if ( u )
! 1529: g = u;
! 1530: else {
! 1531: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1532: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1533: addmd(CO,mod,d,t,&s); d = s;
! 1534: dp_rest(g,&t); g = t;
! 1535: }
! 1536: }
1.55 noro 1537: last:
1.67 ! noro 1538: if ( d )
! 1539: d->sugar = sugar;
! 1540: *rp = d; *dnp = dn;
! 1541: return q;
1.55 noro 1542: }
1543:
1.13 noro 1544: /* nf computation over Z */
1545:
1.20 noro 1546: void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
1.5 noro 1547: {
1.67 ! noro 1548: DP u,p,d,s,t,dmy1;
! 1549: P dmy;
! 1550: NODE l;
! 1551: MP m,mr;
! 1552: int i,n;
! 1553: int *wb;
! 1554: int hmag;
! 1555: int sugar,psugar;
! 1556:
! 1557: if ( !g ) {
! 1558: *rp = 0; return;
! 1559: }
! 1560: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1561: wb = (int *)ALLOCA(n*sizeof(int));
! 1562: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1563: wb[i] = QTOS((Q)BDY(l));
! 1564:
! 1565: hmag = multiple*HMAG(g);
! 1566: sugar = g->sugar;
! 1567:
! 1568: for ( d = 0; g; ) {
! 1569: for ( u = 0, i = 0; i < n; i++ ) {
! 1570: if ( dp_redble(g,p = ps[wb[i]]) ) {
! 1571: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
! 1572: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1573: sugar = MAX(sugar,psugar);
! 1574: if ( !u ) {
! 1575: if ( d )
! 1576: d->sugar = sugar;
! 1577: *rp = d; return;
! 1578: }
! 1579: d = t;
! 1580: break;
! 1581: }
! 1582: }
! 1583: if ( u ) {
! 1584: g = u;
! 1585: if ( d ) {
! 1586: if ( multiple && HMAG(d) > hmag ) {
! 1587: dp_ptozp2(d,g,&t,&u); d = t; g = u;
! 1588: hmag = multiple*HMAG(d);
! 1589: }
! 1590: } else {
! 1591: if ( multiple && HMAG(g) > hmag ) {
! 1592: dp_ptozp(g,&t); g = t;
! 1593: hmag = multiple*HMAG(g);
! 1594: }
! 1595: }
! 1596: }
! 1597: else if ( !full ) {
! 1598: if ( g ) {
! 1599: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1600: }
! 1601: *rp = g; return;
! 1602: } else {
! 1603: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1604: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1605: addd(CO,d,t,&s); d = s;
! 1606: dp_rest(g,&t); g = t;
! 1607:
! 1608: }
! 1609: }
! 1610: if ( d )
! 1611: d->sugar = sugar;
! 1612: *rp = d;
1.5 noro 1613: }
1614:
1.66 noro 1615: void dpm_nf_z(NODE b,DPM g,DPM *ps,int full,int multiple,DPM *rp)
1616: {
1.67 ! noro 1617: DPM u,p,d,s,t;
1.66 noro 1618: DP dmy1;
1.67 ! noro 1619: P dmy;
! 1620: NODE l;
! 1621: DMM m,mr;
! 1622: int i,n;
! 1623: int *wb;
! 1624: int hmag;
! 1625: int sugar,psugar;
! 1626:
! 1627: if ( !g ) {
! 1628: *rp = 0; return;
! 1629: }
! 1630: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1631: wb = (int *)ALLOCA(n*sizeof(int));
! 1632: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1633: wb[i] = QTOS((Q)BDY(l));
! 1634:
! 1635: hmag = multiple*HMAG(g);
! 1636: sugar = g->sugar;
! 1637:
! 1638: for ( d = 0; g; ) {
! 1639: for ( u = 0, i = 0; i < n; i++ ) {
! 1640: if ( dpm_redble(g,p = ps[wb[i]]) ) {
! 1641: dpm_red(d,g,p,&t,&u,&dmy,&dmy1);
! 1642: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1643: sugar = MAX(sugar,psugar);
! 1644: if ( !u ) {
! 1645: if ( d )
! 1646: d->sugar = sugar;
! 1647: *rp = d; return;
! 1648: }
! 1649: d = t;
! 1650: break;
! 1651: }
! 1652: }
! 1653: if ( u ) {
! 1654: g = u;
! 1655: if ( d ) {
! 1656: if ( multiple && HMAG(d) > hmag ) {
! 1657: dpm_ptozp2(d,g,&t,&u); d = t; g = u;
! 1658: hmag = multiple*HMAG(d);
! 1659: }
! 1660: } else {
! 1661: if ( multiple && HMAG(g) > hmag ) {
! 1662: dpm_ptozp(g,&t); g = t;
! 1663: hmag = multiple*HMAG(g);
! 1664: }
! 1665: }
! 1666: }
! 1667: else if ( !full ) {
! 1668: if ( g ) {
! 1669: MKDPM(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1670: }
! 1671: *rp = g; return;
! 1672: } else {
! 1673: m = BDY(g); NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos;
! 1674: NEXT(mr) = 0; MKDPM(g->nv,mr,t); t->sugar = mr->dl->td;
! 1675: adddpm(CO,d,t,&s); d = s;
! 1676: dpm_rest(g,&t); g = t;
! 1677: }
! 1678: }
! 1679: if ( d )
! 1680: d->sugar = sugar;
! 1681: *rp = d;
1.66 noro 1682: }
1683:
1.13 noro 1684: /* nf computation over a field */
1685:
1.20 noro 1686: void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp)
1.13 noro 1687: {
1.67 ! noro 1688: DP u,p,d,s,t;
! 1689: NODE l;
! 1690: MP m,mr;
! 1691: int i,n;
! 1692: int *wb;
! 1693: int sugar,psugar;
! 1694:
! 1695: if ( !g ) {
! 1696: *rp = 0; return;
! 1697: }
! 1698: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1699: wb = (int *)ALLOCA(n*sizeof(int));
! 1700: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1701: wb[i] = QTOS((Q)BDY(l));
! 1702:
! 1703: sugar = g->sugar;
! 1704: for ( d = 0; g; ) {
! 1705: for ( u = 0, i = 0; i < n; i++ ) {
! 1706: if ( dp_redble(g,p = ps[wb[i]]) ) {
! 1707: dp_red_f(g,p,&u);
! 1708: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1709: sugar = MAX(sugar,psugar);
! 1710: if ( !u ) {
! 1711: if ( d )
! 1712: d->sugar = sugar;
! 1713: *rp = d; return;
! 1714: }
! 1715: break;
! 1716: }
! 1717: }
! 1718: if ( u )
! 1719: g = u;
! 1720: else if ( !full ) {
! 1721: if ( g ) {
! 1722: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1723: }
! 1724: *rp = g; return;
! 1725: } else {
! 1726: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1727: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1728: addd(CO,d,t,&s); d = s;
! 1729: dp_rest(g,&t); g = t;
! 1730: }
! 1731: }
! 1732: if ( d )
! 1733: d->sugar = sugar;
! 1734: *rp = d;
1.13 noro 1735: }
1736:
1.66 noro 1737: void dpm_nf_f(NODE b,DPM g,DPM *ps,int full,DPM *rp)
1738: {
1.67 ! noro 1739: DPM u,p,d,s,t;
! 1740: NODE l;
! 1741: DMM m,mr;
! 1742: int i,n;
! 1743: int *wb;
! 1744: int sugar,psugar;
! 1745:
! 1746: if ( !g ) {
! 1747: *rp = 0; return;
! 1748: }
! 1749: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1750: wb = (int *)ALLOCA(n*sizeof(int));
! 1751: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1752: wb[i] = QTOS((Q)BDY(l));
! 1753:
! 1754: sugar = g->sugar;
! 1755: for ( d = 0; g; ) {
! 1756: for ( u = 0, i = 0; i < n; i++ ) {
! 1757: if ( dpm_redble(g,p = ps[wb[i]]) ) {
! 1758: dpm_red_f(g,p,&u);
! 1759: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1760: sugar = MAX(sugar,psugar);
! 1761: if ( !u ) {
! 1762: if ( d )
! 1763: d->sugar = sugar;
! 1764: *rp = d; return;
! 1765: }
! 1766: break;
! 1767: }
! 1768: }
! 1769: if ( u )
! 1770: g = u;
! 1771: else if ( !full ) {
! 1772: if ( g ) {
! 1773: MKDPM(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1774: }
! 1775: *rp = g; return;
! 1776: } else {
! 1777: m = BDY(g); NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos;
! 1778: NEXT(mr) = 0; MKDPM(g->nv,mr,t); t->sugar = mr->dl->td;
! 1779: adddpm(CO,d,t,&s); d = s;
! 1780: dpm_rest(g,&t); g = t;
! 1781: }
! 1782: }
! 1783: if ( d )
! 1784: d->sugar = sugar;
! 1785: *rp = d;
1.66 noro 1786: }
1787:
1.13 noro 1788: /* nf computation over GF(mod) (only for internal use) */
1789:
1.20 noro 1790: void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1791: {
1.67 ! noro 1792: DP u,p,d,s,t;
! 1793: P dmy;
! 1794: NODE l;
! 1795: MP m,mr;
! 1796: int sugar,psugar;
! 1797:
! 1798: if ( !g ) {
! 1799: *rp = 0; return;
! 1800: }
! 1801: sugar = g->sugar;
! 1802: for ( d = 0; g; ) {
! 1803: for ( u = 0, l = b; l; l = NEXT(l) ) {
! 1804: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
! 1805: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
! 1806: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1807: sugar = MAX(sugar,psugar);
! 1808: if ( !u ) {
! 1809: if ( d )
! 1810: d->sugar = sugar;
! 1811: *rp = d; return;
! 1812: }
! 1813: d = t;
! 1814: break;
! 1815: }
! 1816: }
! 1817: if ( u )
! 1818: g = u;
! 1819: else if ( !full ) {
! 1820: if ( g ) {
! 1821: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1822: }
! 1823: *rp = g; return;
! 1824: } else {
! 1825: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1826: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1827: addmd(CO,mod,d,t,&s); d = s;
! 1828: dp_rest(g,&t); g = t;
! 1829: }
! 1830: }
! 1831: if ( d )
! 1832: d->sugar = sugar;
! 1833: *rp = d;
1.5 noro 1834: }
1835:
1.20 noro 1836: void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp)
1.5 noro 1837: {
1.67 ! noro 1838: DP u,p,d,s,t;
! 1839: NODE l;
! 1840: MP m,mr;
! 1841: int i,n;
! 1842: int *wb;
! 1843: int sugar,psugar;
! 1844: P dn,tdn,tdn1;
! 1845:
! 1846: dn = (P)ONEM;
! 1847: if ( !g ) {
! 1848: *rp = 0; *dnp = dn; return;
! 1849: }
! 1850: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1851: wb = (int *)ALLOCA(n*sizeof(int));
! 1852: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1853: wb[i] = QTOS((Q)BDY(l));
! 1854: sugar = g->sugar;
! 1855: for ( d = 0; g; ) {
! 1856: for ( u = 0, i = 0; i < n; i++ ) {
! 1857: if ( dp_redble(g,p = ps[wb[i]]) ) {
! 1858: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
! 1859: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1860: sugar = MAX(sugar,psugar);
! 1861: if ( !u ) {
! 1862: if ( d )
! 1863: d->sugar = sugar;
! 1864: *rp = d; *dnp = dn; return;
! 1865: } else {
! 1866: d = t;
! 1867: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
! 1868: }
! 1869: break;
! 1870: }
! 1871: }
! 1872: if ( u )
! 1873: g = u;
! 1874: else if ( !full ) {
! 1875: if ( g ) {
! 1876: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1877: }
! 1878: *rp = g; *dnp = dn; return;
! 1879: } else {
! 1880: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 1881: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 1882: addmd(CO,mod,d,t,&s); d = s;
! 1883: dp_rest(g,&t); g = t;
! 1884: }
! 1885: }
! 1886: if ( d )
! 1887: d->sugar = sugar;
! 1888: *rp = d; *dnp = dn;
1.5 noro 1889: }
1890:
1.20 noro 1891: void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1892: {
1.67 ! noro 1893: DP u,p,d;
! 1894: NODE l;
! 1895: MP m,mrd;
! 1896: int sugar,psugar,n,h_reducible;
! 1897:
! 1898: if ( !g ) {
! 1899: *rp = 0; return;
! 1900: }
! 1901: sugar = g->sugar;
! 1902: n = g->nv;
! 1903: for ( d = 0; g; ) {
! 1904: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
! 1905: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
! 1906: h_reducible = 1;
! 1907: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1908: _dp_red_mod_destructive(g,p,mod,&u); g = u;
! 1909: sugar = MAX(sugar,psugar);
! 1910: if ( !g ) {
! 1911: if ( d )
! 1912: d->sugar = sugar;
! 1913: _dptodp(d,rp); _free_dp(d); return;
! 1914: }
! 1915: break;
! 1916: }
! 1917: }
! 1918: if ( !h_reducible ) {
! 1919: /* head term is not reducible */
! 1920: if ( !full ) {
! 1921: if ( g )
! 1922: g->sugar = sugar;
! 1923: _dptodp(g,rp); _free_dp(g); return;
! 1924: } else {
! 1925: m = BDY(g);
! 1926: if ( NEXT(m) ) {
! 1927: BDY(g) = NEXT(m); NEXT(m) = 0;
! 1928: } else {
! 1929: _FREEDP(g); g = 0;
! 1930: }
! 1931: if ( d ) {
! 1932: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
! 1933: NEXT(mrd) = m;
! 1934: } else {
! 1935: _MKDP(n,m,d);
! 1936: }
! 1937: }
! 1938: }
! 1939: }
! 1940: if ( d )
! 1941: d->sugar = sugar;
! 1942: _dptodp(d,rp); _free_dp(d);
1.5 noro 1943: }
1.13 noro 1944:
1945: /* reduction by linear base over a field */
1946:
1.20 noro 1947: void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p)
1.13 noro 1948: {
1.67 ! noro 1949: DP r1,r2,b1,b2,t,s;
! 1950: Obj c,c1,c2;
! 1951: NODE l,b;
! 1952: int n;
! 1953:
! 1954: if ( !p1 ) {
! 1955: *r1p = p1; *r2p = p2; return;
! 1956: }
! 1957: n = p1->nv;
! 1958: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
! 1959: if ( !r1 ) {
! 1960: *r1p = r1; *r2p = r2; return;
! 1961: }
! 1962: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
! 1963: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
! 1964: b2 = (DP)BDY(NEXT(b));
! 1965: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
! 1966: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
! 1967: muldc(CO,b1,(Obj)c,&t); addd(CO,r1,t,&s); r1 = s;
! 1968: muldc(CO,b2,(Obj)c,&t); addd(CO,r2,t,&s); r2 = s;
! 1969: }
! 1970: }
! 1971: *r1p = r1; *r2p = r2;
1.13 noro 1972: }
1973:
1974: /* reduction by linear base over GF(mod) */
1.5 noro 1975:
1.20 noro 1976: void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p)
1.5 noro 1977: {
1.67 ! noro 1978: DP r1,r2,b1,b2,t,s;
! 1979: P c;
! 1980: MQ c1,c2;
! 1981: NODE l,b;
! 1982: int n;
! 1983:
! 1984: if ( !p1 ) {
! 1985: *r1p = p1; *r2p = p2; return;
! 1986: }
! 1987: n = p1->nv;
! 1988: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
! 1989: if ( !r1 ) {
! 1990: *r1p = r1; *r2p = r2; return;
! 1991: }
! 1992: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
! 1993: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
! 1994: b2 = (DP)BDY(NEXT(b));
! 1995: invmq(mod,(MQ)BDY(b1)->c,&c1);
! 1996: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
! 1997: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
! 1998: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
! 1999: }
! 2000: }
! 2001: *r1p = r1; *r2p = r2;
1.5 noro 2002: }
2003:
1.20 noro 2004: void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp)
1.5 noro 2005: {
1.67 ! noro 2006: DP s,t,u;
! 2007: MP m;
! 2008: DL h;
! 2009: int i,n;
! 2010:
! 2011: if ( !p ) {
! 2012: *rp = p; return;
! 2013: }
! 2014: n = p->nv;
! 2015: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2016: h = m->dl;
! 2017: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
! 2018: i++;
! 2019: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),(P)m->c,&t);
! 2020: addmd(CO,mod,s,t,&u); s = u;
! 2021: }
! 2022: *rp = s;
1.24 noro 2023: }
2024:
2025: void dp_nf_tab_f(DP p,LIST *tab,DP *rp)
2026: {
1.67 ! noro 2027: DP s,t,u;
! 2028: MP m;
! 2029: DL h;
! 2030: int i,n;
! 2031:
! 2032: if ( !p ) {
! 2033: *rp = p; return;
! 2034: }
! 2035: n = p->nv;
! 2036: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2037: h = m->dl;
! 2038: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
! 2039: i++;
! 2040: muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
! 2041: addd(CO,s,t,&u); s = u;
! 2042: }
! 2043: *rp = s;
1.5 noro 2044: }
2045:
1.7 noro 2046: /*
2047: * setting flags
1.30 noro 2048: * call create_order_spec with vl=0 to set old type order.
1.7 noro 2049: *
2050: */
2051:
1.27 noro 2052: int create_order_spec(VL vl,Obj obj,struct order_spec **specp)
1.5 noro 2053: {
1.67 ! noro 2054: int i,j,n,s,row,col,ret,wlen;
! 2055: struct order_spec *spec;
! 2056: struct order_pair *l;
1.62 noro 2057: Obj wp,wm;
1.67 ! noro 2058: NODE node,t,tn,wpair;
! 2059: MAT m;
! 2060: VECT v;
! 2061: pointer **b,*bv;
! 2062: int **w;
! 2063:
! 2064: if ( vl && obj && OID(obj) == O_LIST ) {
! 2065: ret = create_composite_order_spec(vl,(LIST)obj,specp);
! 2066: if ( show_orderspec )
! 2067: print_composite_order_spec(*specp);
! 2068: return ret;
! 2069: }
! 2070:
! 2071: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
! 2072: if ( !obj || NUM(obj) ) {
! 2073: spec->id = 0; spec->obj = obj;
! 2074: spec->ord.simple = QTOS((Q)obj);
! 2075: return 1;
! 2076: } else if ( OID(obj) == O_LIST ) {
1.62 noro 2077: /* module order; obj = [0|1,w,ord] or [0|1,ord] */
1.67 ! noro 2078: node = BDY((LIST)obj);
! 2079: if ( !BDY(node) || NUM(BDY(node)) ) {
1.62 noro 2080: switch ( length(node) ) {
2081: case 2:
1.67 ! noro 2082: create_order_spec(0,(Obj)BDY(NEXT(node)),&spec);
! 2083: spec->id += 256; spec->obj = obj;
1.62 noro 2084: spec->top_weight = 0;
2085: spec->module_rank = 0;
2086: spec->module_top_weight = 0;
1.67 ! noro 2087: spec->ispot = (BDY(node)!=0);
! 2088: if ( spec->ispot ) {
! 2089: n = QTOS((Q)BDY(node));
! 2090: if ( n < 0 )
! 2091: spec->pot_nelim = -n;
! 2092: else
! 2093: spec->pot_nelim = 0;
! 2094: }
1.62 noro 2095: break;
2096:
2097: case 3:
1.67 ! noro 2098: create_order_spec(0,(Obj)BDY(NEXT(NEXT(node))),&spec);
! 2099: spec->id += 256; spec->obj = obj;
! 2100: spec->ispot = (BDY(node)!=0);
1.62 noro 2101: node = NEXT(node);
2102: if ( !BDY(node) || OID(BDY(node)) != O_LIST )
1.67 ! noro 2103: error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
1.62 noro 2104: wpair = BDY((LIST)BDY(node));
2105: if ( length(wpair) != 2 )
1.67 ! noro 2106: error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
1.62 noro 2107:
2108: wp = BDY(wpair);
2109: wm = BDY(NEXT(wpair));
2110: if ( !wp || OID(wp) != O_LIST || !wm || OID(wm) != O_LIST )
1.67 ! noro 2111: error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
1.62 noro 2112: spec->nv = length(BDY((LIST)wp));
2113: spec->top_weight = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
1.67 ! noro 2114: for ( i = 0, t = BDY((LIST)wp); i < spec->nv; t = NEXT(t), i++ )
1.62 noro 2115: spec->top_weight[i] = QTOS((Q)BDY(t));
2116:
2117: spec->module_rank = length(BDY((LIST)wm));
2118: spec->module_top_weight = (int *)MALLOC_ATOMIC(spec->module_rank*sizeof(int));
1.67 ! noro 2119: for ( i = 0, t = BDY((LIST)wm); i < spec->module_rank; t = NEXT(t), i++ )
1.62 noro 2120: spec->module_top_weight[i] = QTOS((Q)BDY(t));
2121: break;
2122: default:
1.67 ! noro 2123: error("create_order_spec : invalid arguments for module order");
1.62 noro 2124: }
2125:
1.67 ! noro 2126: *specp = spec;
! 2127: return 1;
! 2128: } else {
1.62 noro 2129: /* block order in polynomial ring */
1.67 ! noro 2130: for ( n = 0, t = node; t; t = NEXT(t), n++ );
! 2131: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
! 2132: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
! 2133: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
! 2134: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
! 2135: s += l[i].length;
! 2136: }
! 2137: spec->id = 1; spec->obj = obj;
! 2138: spec->ord.block.order_pair = l;
! 2139: spec->ord.block.length = n; spec->nv = s;
! 2140: return 1;
! 2141: }
! 2142: } else if ( OID(obj) == O_MAT ) {
! 2143: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
! 2144: w = almat(row,col);
! 2145: for ( i = 0; i < row; i++ )
! 2146: for ( j = 0; j < col; j++ )
! 2147: w[i][j] = QTOS((Q)b[i][j]);
! 2148: spec->id = 2; spec->obj = obj;
! 2149: spec->nv = col; spec->ord.matrix.row = row;
! 2150: spec->ord.matrix.matrix = w;
! 2151: return 1;
! 2152: } else
! 2153: return 0;
1.5 noro 2154: }
2155:
1.28 noro 2156: void print_composite_order_spec(struct order_spec *spec)
2157: {
1.67 ! noro 2158: int nv,n,len,i,j,k,start;
! 2159: struct weight_or_block *worb;
1.28 noro 2160:
1.67 ! noro 2161: nv = spec->nv;
! 2162: n = spec->ord.composite.length;
! 2163: worb = spec->ord.composite.w_or_b;
! 2164: for ( i = 0; i < n; i++, worb++ ) {
! 2165: len = worb->length;
! 2166: printf("[ ");
! 2167: switch ( worb->type ) {
! 2168: case IS_DENSE_WEIGHT:
! 2169: for ( j = 0; j < len; j++ )
! 2170: printf("%d ",worb->body.dense_weight[j]);
! 2171: for ( ; j < nv; j++ )
! 2172: printf("0 ");
! 2173: break;
! 2174: case IS_SPARSE_WEIGHT:
! 2175: for ( j = 0, k = 0; j < nv; j++ )
! 2176: if ( j == worb->body.sparse_weight[k].pos )
! 2177: printf("%d ",worb->body.sparse_weight[k++].value);
! 2178: else
! 2179: printf("0 ");
! 2180: break;
! 2181: case IS_BLOCK:
! 2182: start = worb->body.block.start;
! 2183: for ( j = 0; j < start; j++ ) printf("0 ");
! 2184: switch ( worb->body.block.order ) {
! 2185: case 0:
! 2186: for ( k = 0; k < len; k++, j++ ) printf("R ");
! 2187: break;
! 2188: case 1:
! 2189: for ( k = 0; k < len; k++, j++ ) printf("G ");
! 2190: break;
! 2191: case 2:
! 2192: for ( k = 0; k < len; k++, j++ ) printf("L ");
! 2193: break;
! 2194: }
! 2195: for ( ; j < nv; j++ ) printf("0 ");
! 2196: break;
! 2197: }
! 2198: printf("]\n");
! 2199: }
1.38 noro 2200: }
2201:
2202: struct order_spec *append_block(struct order_spec *spec,
1.67 ! noro 2203: int nv,int nalg,int ord)
1.38 noro 2204: {
1.67 ! noro 2205: MAT m,mat;
! 2206: int i,j,row,col,n;
! 2207: Q **b,**wp;
! 2208: int **w;
! 2209: NODE t,s,s0;
! 2210: struct order_pair *l,*l0;
! 2211: int n0,nv0;
! 2212: LIST list0,list1,list;
! 2213: Q oq,nq;
! 2214: struct order_spec *r;
! 2215:
! 2216: r = (struct order_spec *)MALLOC(sizeof(struct order_spec));
! 2217: switch ( spec->id ) {
! 2218: case 0:
! 2219: STOQ(spec->ord.simple,oq); STOQ(nv,nq);
! 2220: t = mknode(2,oq,nq); MKLIST(list0,t);
! 2221: STOQ(ord,oq); STOQ(nalg,nq);
! 2222: t = mknode(2,oq,nq); MKLIST(list1,t);
! 2223: t = mknode(2,list0,list1); MKLIST(list,t);
! 2224: l = (struct order_pair *)MALLOC_ATOMIC(2*sizeof(struct order_pair));
! 2225: l[0].order = spec->ord.simple; l[0].length = nv;
! 2226: l[1].order = ord; l[1].length = nalg;
! 2227: r->id = 1; r->obj = (Obj)list;
! 2228: r->ord.block.order_pair = l;
! 2229: r->ord.block.length = 2;
! 2230: r->nv = nv+nalg;
! 2231: break;
! 2232: case 1:
! 2233: if ( spec->nv != nv )
! 2234: error("append_block : number of variables mismatch");
! 2235: l0 = spec->ord.block.order_pair;
! 2236: n0 = spec->ord.block.length;
! 2237: nv0 = spec->nv;
! 2238: list0 = (LIST)spec->obj;
! 2239: n = n0+1;
! 2240: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
! 2241: for ( i = 0; i < n0; i++ )
! 2242: l[i] = l0[i];
! 2243: l[i].order = ord; l[i].length = nalg;
! 2244: for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) {
! 2245: NEXTNODE(s0,s); BDY(s) = BDY(t);
! 2246: }
! 2247: STOQ(ord,oq); STOQ(nalg,nq);
! 2248: t = mknode(2,oq,nq); MKLIST(list,t);
! 2249: NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0;
! 2250: MKLIST(list,s0);
! 2251: r->id = 1; r->obj = (Obj)list;
! 2252: r->ord.block.order_pair = l;
! 2253: r->ord.block.length = n;
! 2254: r->nv = nv+nalg;
! 2255: break;
! 2256: case 2:
! 2257: if ( spec->nv != nv )
! 2258: error("append_block : number of variables mismatch");
! 2259: m = (MAT)spec->obj;
! 2260: row = m->row; col = m->col; b = (Q **)BDY(m);
! 2261: w = almat(row+nalg,col+nalg);
! 2262: MKMAT(mat,row+nalg,col+nalg); wp = (Q **)BDY(mat);
! 2263: for ( i = 0; i < row; i++ )
! 2264: for ( j = 0; j < col; j++ ) {
! 2265: w[i][j] = QTOS(b[i][j]);
! 2266: wp[i][j] = b[i][j];
! 2267: }
! 2268: for ( i = 0; i < nalg; i++ ) {
! 2269: w[i+row][i+col] = 1;
! 2270: wp[i+row][i+col] = ONE;
! 2271: }
! 2272: r->id = 2; r->obj = (Obj)mat;
! 2273: r->nv = col+nalg; r->ord.matrix.row = row+nalg;
! 2274: r->ord.matrix.matrix = w;
! 2275: break;
! 2276: case 3:
! 2277: default:
! 2278: /* XXX */
! 2279: error("append_block : not implemented yet");
! 2280: }
! 2281: return r;
1.28 noro 2282: }
2283:
1.37 noro 2284: int comp_sw(struct sparse_weight *a, struct sparse_weight *b)
2285: {
1.67 ! noro 2286: if ( a->pos > b->pos ) return 1;
! 2287: else if ( a->pos < b->pos ) return -1;
! 2288: else return 0;
1.37 noro 2289: }
2290:
1.27 noro 2291: /* order = [w_or_b, w_or_b, ... ] */
2292: /* w_or_b = w or b */
2293: /* w = [1,2,...] or [x,1,y,2,...] */
2294: /* b = [@lex,x,y,...,z] etc */
2295:
2296: int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp)
2297: {
1.67 ! noro 2298: NODE wb,t,p;
! 2299: struct order_spec *spec;
! 2300: VL tvl;
! 2301: int n,i,j,k,l,start,end,len,w;
! 2302: int *dw;
! 2303: struct sparse_weight *sw;
! 2304: struct weight_or_block *w_or_b;
! 2305: Obj a0;
! 2306: NODE a;
! 2307: V v,sv,ev;
! 2308: SYMBOL sym;
! 2309: int *top;
! 2310:
! 2311: /* l = number of vars in vl */
! 2312: for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ );
! 2313: /* n = number of primitives in order */
! 2314: wb = BDY(order);
! 2315: n = length(wb);
! 2316: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
! 2317: spec->id = 3;
! 2318: spec->obj = (Obj)order;
! 2319: spec->nv = l;
! 2320: spec->ord.composite.length = n;
! 2321: w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *)
! 2322: MALLOC(sizeof(struct weight_or_block)*(n+1));
! 2323:
! 2324: /* top : register the top variable in each w_or_b specification */
! 2325: top = (int *)ALLOCA(l*sizeof(int));
! 2326: for ( i = 0; i < l; i++ ) top[i] = 0;
! 2327:
! 2328: for ( t = wb, i = 0; t; t = NEXT(t), i++ ) {
! 2329: if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST )
! 2330: error("a list of lists must be specified for the key \"order\"");
! 2331: a = BDY((LIST)BDY(t));
! 2332: len = length(a);
! 2333: a0 = (Obj)BDY(a);
! 2334: if ( !a0 || OID(a0) == O_N ) {
! 2335: /* a is a dense weight vector */
! 2336: dw = (int *)MALLOC(sizeof(int)*len);
! 2337: for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) {
! 2338: if ( !INT((Q)BDY(p)) )
! 2339: error("a dense weight vector must be specified as a list of integers");
! 2340: dw[j] = QTOS((Q)BDY(p));
! 2341: }
! 2342: w_or_b[i].type = IS_DENSE_WEIGHT;
! 2343: w_or_b[i].length = len;
! 2344: w_or_b[i].body.dense_weight = dw;
! 2345:
! 2346: /* find the top */
! 2347: for ( k = 0; k < len && !dw[k]; k++ );
! 2348: if ( k < len ) top[k] = 1;
! 2349:
! 2350: } else if ( OID(a0) == O_P ) {
! 2351: /* a is a sparse weight vector */
! 2352: len >>= 1;
! 2353: sw = (struct sparse_weight *)
! 2354: MALLOC(sizeof(struct sparse_weight)*len);
! 2355: for ( j = 0, p = a; j < len; j++ ) {
! 2356: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
! 2357: error("a sparse weight vector must be specified as [var1,weight1,...]");
! 2358: v = VR((P)BDY(p)); p = NEXT(p);
! 2359: for ( tvl = vl, k = 0; tvl && tvl->v != v;
! 2360: k++, tvl = NEXT(tvl) );
! 2361: if ( !tvl )
! 2362: error("invalid variable name in a sparse weight vector");
! 2363: sw[j].pos = k;
! 2364: if ( !INT((Q)BDY(p)) )
! 2365: error("a sparse weight vector must be specified as [var1,weight1,...]");
! 2366: sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p);
! 2367: }
! 2368: qsort(sw,len,sizeof(struct sparse_weight),
! 2369: (int (*)(const void *,const void *))comp_sw);
! 2370: w_or_b[i].type = IS_SPARSE_WEIGHT;
! 2371: w_or_b[i].length = len;
! 2372: w_or_b[i].body.sparse_weight = sw;
! 2373:
! 2374: /* find the top */
! 2375: for ( k = 0; k < len && !sw[k].value; k++ );
! 2376: if ( k < len ) top[sw[k].pos] = 1;
! 2377: } else if ( OID(a0) == O_RANGE ) {
! 2378: /* [range(v1,v2),w] */
! 2379: sv = VR((P)(((RANGE)a0)->start));
! 2380: ev = VR((P)(((RANGE)a0)->end));
! 2381: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
! 2382: if ( !tvl )
! 2383: error("invalid range");
! 2384: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
! 2385: if ( !tvl )
! 2386: error("invalid range");
! 2387: len = end-start+1;
! 2388: sw = (struct sparse_weight *)
! 2389: MALLOC(sizeof(struct sparse_weight)*len);
! 2390: w = QTOS((Q)BDY(NEXT(a)));
! 2391: for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) );
! 2392: for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) {
! 2393: sw[j].pos = k;
! 2394: sw[j].value = w;
! 2395: }
! 2396: w_or_b[i].type = IS_SPARSE_WEIGHT;
! 2397: w_or_b[i].length = len;
! 2398: w_or_b[i].body.sparse_weight = sw;
! 2399:
! 2400: /* register the top */
! 2401: if ( w ) top[start] = 1;
! 2402: } else if ( OID(a0) == O_SYMBOL ) {
! 2403: /* a is a block */
! 2404: sym = (SYMBOL)a0; a = NEXT(a); len--;
! 2405: if ( OID((Obj)BDY(a)) == O_RANGE ) {
! 2406: sv = VR((P)(((RANGE)BDY(a))->start));
! 2407: ev = VR((P)(((RANGE)BDY(a))->end));
! 2408: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
! 2409: if ( !tvl )
! 2410: error("invalid range");
! 2411: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
! 2412: if ( !tvl )
! 2413: error("invalid range");
! 2414: len = end-start+1;
! 2415: } else {
! 2416: for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a));
! 2417: tvl = NEXT(tvl), start++ );
! 2418: for ( p = NEXT(a), tvl = NEXT(tvl); p;
! 2419: p = NEXT(p), tvl = NEXT(tvl) ) {
! 2420: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
! 2421: error("a block must be specified as [ordsymbol,var1,var2,...]");
! 2422: if ( tvl->v != VR((P)BDY(p)) ) break;
! 2423: }
! 2424: if ( p )
! 2425: error("a block must be contiguous in the variable list");
! 2426: }
! 2427: w_or_b[i].type = IS_BLOCK;
! 2428: w_or_b[i].length = len;
! 2429: w_or_b[i].body.block.start = start;
! 2430: if ( !strcmp(sym->name,"@grlex") )
! 2431: w_or_b[i].body.block.order = 0;
! 2432: else if ( !strcmp(sym->name,"@glex") )
! 2433: w_or_b[i].body.block.order = 1;
! 2434: else if ( !strcmp(sym->name,"@lex") )
! 2435: w_or_b[i].body.block.order = 2;
! 2436: else
! 2437: error("invalid ordername");
! 2438: /* register the tops */
! 2439: for ( j = 0, k = start; j < len; j++, k++ )
! 2440: top[k] = 1;
! 2441: }
! 2442: }
! 2443: for ( k = 0; k < l && top[k]; k++ );
! 2444: if ( k < l ) {
! 2445: /* incomplete order specification; add @grlex */
! 2446: w_or_b[n].type = IS_BLOCK;
! 2447: w_or_b[n].length = l;
! 2448: w_or_b[n].body.block.start = 0;
! 2449: w_or_b[n].body.block.order = 0;
! 2450: spec->ord.composite.length = n+1;
! 2451: }
1.27 noro 2452: }
2453:
1.35 noro 2454: /* module order spec */
2455:
2456: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s)
2457: {
1.67 ! noro 2458: struct modorder_spec *spec;
! 2459: NODE n,t;
! 2460: LIST list;
! 2461: int *ds;
! 2462: int i,l;
! 2463: Q q;
! 2464:
! 2465: *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec));
! 2466: spec->id = id;
! 2467: if ( shift ) {
! 2468: n = BDY(shift);
! 2469: spec->len = l = length(n);
! 2470: spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int));
! 2471: for ( t = n, i = 0; t; t = NEXT(t), i++ )
! 2472: ds[i] = QTOS((Q)BDY(t));
! 2473: } else {
! 2474: spec->len = 0;
! 2475: spec->degree_shift = 0;
! 2476: }
! 2477: STOQ(id,q);
! 2478: n = mknode(2,q,shift);
! 2479: MKLIST(list,n);
! 2480: spec->obj = (Obj)list;
1.35 noro 2481: }
2482:
1.7 noro 2483: /*
2484: * converters
2485: *
2486: */
2487:
1.20 noro 2488: void dp_homo(DP p,DP *rp)
1.5 noro 2489: {
1.67 ! noro 2490: MP m,mr,mr0;
! 2491: int i,n,nv,td;
! 2492: DL dl,dlh;
! 2493:
! 2494: if ( !p )
! 2495: *rp = 0;
! 2496: else {
! 2497: n = p->nv; nv = n + 1;
! 2498: m = BDY(p); td = sugard(m);
! 2499: for ( mr0 = 0; m; m = NEXT(m) ) {
! 2500: NEXTMP(mr0,mr); mr->c = m->c;
! 2501: dl = m->dl;
! 2502: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
! 2503: dlh->td = td;
! 2504: for ( i = 0; i < n; i++ )
! 2505: dlh->d[i] = dl->d[i];
! 2506: dlh->d[n] = td - dl->td;
! 2507: }
! 2508: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 2509: }
1.5 noro 2510: }
2511:
1.20 noro 2512: void dp_dehomo(DP p,DP *rp)
1.5 noro 2513: {
1.67 ! noro 2514: MP m,mr,mr0;
! 2515: int i,n,nv;
! 2516: DL dl,dlh;
! 2517:
! 2518: if ( !p )
! 2519: *rp = 0;
! 2520: else {
! 2521: n = p->nv; nv = n - 1;
! 2522: m = BDY(p);
! 2523: for ( mr0 = 0; m; m = NEXT(m) ) {
! 2524: NEXTMP(mr0,mr); mr->c = m->c;
! 2525: dlh = m->dl;
! 2526: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
! 2527: dl->td = dlh->td - dlh->d[nv];
! 2528: for ( i = 0; i < nv; i++ )
! 2529: dl->d[i] = dlh->d[i];
! 2530: }
! 2531: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 2532: }
1.5 noro 2533: }
2534:
1.20 noro 2535: void dp_mod(DP p,int mod,NODE subst,DP *rp)
1.5 noro 2536: {
1.67 ! noro 2537: MP m,mr,mr0;
! 2538: P t,s,s1;
! 2539: V v;
! 2540: NODE tn;
! 2541:
! 2542: if ( !p )
! 2543: *rp = 0;
! 2544: else {
! 2545: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2546: for ( tn = subst, s = (P)m->c; tn; tn = NEXT(tn) ) {
! 2547: v = VR((P)BDY(tn)); tn = NEXT(tn);
! 2548: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
! 2549: }
! 2550: ptomp(mod,s,&t);
! 2551: if ( t ) {
! 2552: NEXTMP(mr0,mr); mr->c = (Obj)t; mr->dl = m->dl;
! 2553: }
! 2554: }
! 2555: if ( mr0 ) {
! 2556: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 2557: } else
! 2558: *rp = 0;
! 2559: }
1.5 noro 2560: }
2561:
1.20 noro 2562: void dp_rat(DP p,DP *rp)
1.5 noro 2563: {
1.67 ! noro 2564: MP m,mr,mr0;
1.5 noro 2565:
1.67 ! noro 2566: if ( !p )
! 2567: *rp = 0;
! 2568: else {
! 2569: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 2570: NEXTMP(mr0,mr); mptop((P)m->c,(P *)&mr->c); mr->dl = m->dl;
! 2571: }
! 2572: if ( mr0 ) {
! 2573: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 2574: } else
! 2575: *rp = 0;
! 2576: }
1.5 noro 2577: }
2578:
2579:
1.27 noro 2580: void homogenize_order(struct order_spec *old,int n,struct order_spec **newp)
1.5 noro 2581: {
1.67 ! noro 2582: struct order_pair *l;
! 2583: int length,nv,row,i,j;
! 2584: int **newm,**oldm;
! 2585: struct order_spec *new;
! 2586: int onv,nnv,nlen,olen,owlen;
! 2587: struct weight_or_block *owb,*nwb;
! 2588:
! 2589: *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec));
! 2590: switch ( old->id ) {
! 2591: case 0:
! 2592: switch ( old->ord.simple ) {
! 2593: case 0:
! 2594: new->id = 0; new->ord.simple = 0; break;
! 2595: case 1:
! 2596: l = (struct order_pair *)
! 2597: MALLOC_ATOMIC(2*sizeof(struct order_pair));
! 2598: l[0].length = n; l[0].order = 1;
! 2599: l[1].length = 1; l[1].order = 2;
! 2600: new->id = 1;
! 2601: new->ord.block.order_pair = l;
! 2602: new->ord.block.length = 2; new->nv = n+1;
! 2603: break;
! 2604: case 2:
! 2605: new->id = 0; new->ord.simple = 1; break;
! 2606: case 3: case 4: case 5:
! 2607: new->id = 0; new->ord.simple = old->ord.simple+3;
! 2608: dp_nelim = n-1; break;
! 2609: case 6: case 7: case 8: case 9:
! 2610: new->id = 0; new->ord.simple = old->ord.simple; break;
! 2611: default:
! 2612: error("homogenize_order : invalid input");
! 2613: }
! 2614: break;
! 2615: case 1: case 257:
! 2616: length = old->ord.block.length;
! 2617: l = (struct order_pair *)
! 2618: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
! 2619: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
! 2620: l[length].order = 2; l[length].length = 1;
! 2621: new->id = old->id; new->nv = n+1;
! 2622: new->ord.block.order_pair = l;
! 2623: new->ord.block.length = length+1;
! 2624: new->ispot = old->ispot;
! 2625: break;
! 2626: case 2: case 258:
! 2627: nv = old->nv; row = old->ord.matrix.row;
! 2628: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
! 2629: for ( i = 0; i <= nv; i++ )
! 2630: newm[0][i] = 1;
! 2631: for ( i = 0; i < row; i++ ) {
! 2632: for ( j = 0; j < nv; j++ )
! 2633: newm[i+1][j] = oldm[i][j];
! 2634: newm[i+1][j] = 0;
! 2635: }
! 2636: new->id = old->id; new->nv = nv+1;
! 2637: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
! 2638: new->ispot = old->ispot;
! 2639: break;
! 2640: case 3: case 259:
! 2641: onv = old->nv;
! 2642: nnv = onv+1;
! 2643: olen = old->ord.composite.length;
! 2644: nlen = olen+1;
! 2645: owb = old->ord.composite.w_or_b;
! 2646: nwb = (struct weight_or_block *)
! 2647: MALLOC(nlen*sizeof(struct weight_or_block));
! 2648: for ( i = 0; i < olen; i++ ) {
! 2649: nwb[i].type = owb[i].type;
! 2650: switch ( owb[i].type ) {
! 2651: case IS_DENSE_WEIGHT:
! 2652: owlen = owb[i].length;
! 2653: nwb[i].length = owlen+1;
! 2654: nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int));
! 2655: for ( j = 0; j < owlen; j++ )
! 2656: nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j];
! 2657: nwb[i].body.dense_weight[owlen] = 0;
! 2658: break;
! 2659: case IS_SPARSE_WEIGHT:
! 2660: nwb[i].length = owb[i].length;
! 2661: nwb[i].body.sparse_weight = owb[i].body.sparse_weight;
! 2662: break;
! 2663: case IS_BLOCK:
! 2664: nwb[i].length = owb[i].length;
! 2665: nwb[i].body.block = owb[i].body.block;
! 2666: break;
! 2667: }
! 2668: }
! 2669: nwb[i].type = IS_SPARSE_WEIGHT;
! 2670: nwb[i].body.sparse_weight =
! 2671: (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
! 2672: nwb[i].body.sparse_weight[0].pos = onv;
! 2673: nwb[i].body.sparse_weight[0].value = 1;
! 2674: new->id = old->id;
! 2675: new->nv = nnv;
! 2676: new->ord.composite.length = nlen;
! 2677: new->ord.composite.w_or_b = nwb;
! 2678: new->ispot = old->ispot;
! 2679: print_composite_order_spec(new);
! 2680: break;
! 2681: case 256: /* simple module order */
! 2682: switch ( old->ord.simple ) {
! 2683: case 0:
! 2684: new->id = 256; new->ord.simple = 0; break;
! 2685: case 1:
! 2686: l = (struct order_pair *)
! 2687: MALLOC_ATOMIC(2*sizeof(struct order_pair));
! 2688: l[0].length = n; l[0].order = old->ord.simple;
! 2689: l[1].length = 1; l[1].order = 2;
! 2690: new->id = 257;
! 2691: new->ord.block.order_pair = l;
! 2692: new->ord.block.length = 2; new->nv = n+1;
! 2693: break;
! 2694: case 2:
! 2695: new->id = 256; new->ord.simple = 1; break;
! 2696: default:
! 2697: error("homogenize_order : invalid input");
! 2698: }
! 2699: new->ispot = old->ispot;
! 2700: break;
! 2701: default:
! 2702: error("homogenize_order : invalid input");
! 2703: }
1.7 noro 2704: }
2705:
1.20 noro 2706: void qltozl(Q *w,int n,Q *dvr)
1.7 noro 2707: {
1.67 ! noro 2708: N nm,dn;
! 2709: N g,l1,l2,l3;
! 2710: Q c,d;
! 2711: int i;
! 2712: struct oVECT v;
! 2713:
! 2714: for ( i = 0; i < n; i++ )
! 2715: if ( w[i] && !INT(w[i]) )
! 2716: break;
! 2717: if ( i == n ) {
! 2718: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
! 2719: igcdv(&v,dvr); return;
! 2720: }
! 2721: for ( i = 0; !w[i]; i++ );
! 2722: c = w[i]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
! 2723: for ( i++; i < n; i++ ) {
! 2724: c = w[i];
! 2725: if ( !c ) continue;
! 2726: l1 = INT(c) ? ONEN : DN(c);
! 2727: gcdn(nm,NM(c),&g); nm = g;
! 2728: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
! 2729: }
! 2730: if ( UNIN(dn) )
! 2731: NTOQ(nm,1,d);
! 2732: else
! 2733: NDTOQ(nm,dn,1,d);
! 2734: *dvr = d;
1.7 noro 2735: }
1.5 noro 2736:
1.20 noro 2737: int comp_nm(Q *a,Q *b)
1.7 noro 2738: {
1.67 ! noro 2739: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
1.7 noro 2740: }
2741:
1.20 noro 2742: void sortbynm(Q *w,int n)
1.7 noro 2743: {
1.67 ! noro 2744: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
1.7 noro 2745: }
1.5 noro 2746:
2747:
1.7 noro 2748: /*
2749: * simple operations
2750: *
2751: */
1.5 noro 2752:
1.20 noro 2753: int dp_redble(DP p1,DP p2)
1.7 noro 2754: {
1.67 ! noro 2755: int i,n;
! 2756: DL d1,d2;
1.5 noro 2757:
1.67 ! noro 2758: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 2759: if ( d1->td < d2->td )
! 2760: return 0;
! 2761: else {
! 2762: for ( i = 0, n = p1->nv; i < n; i++ )
! 2763: if ( d1->d[i] < d2->d[i] )
! 2764: return 0;
! 2765: return 1;
! 2766: }
1.5 noro 2767: }
2768:
1.66 noro 2769: int dpm_redble(DPM p1,DPM p2)
2770: {
1.67 ! noro 2771: int i,n;
! 2772: DL d1,d2;
1.66 noro 2773:
2774: if ( BDY(p1)->pos != BDY(p2)->pos ) return 0;
1.67 ! noro 2775: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 2776: if ( d1->td < d2->td )
! 2777: return 0;
! 2778: else {
! 2779: for ( i = 0, n = p1->nv; i < n; i++ )
! 2780: if ( d1->d[i] < d2->d[i] )
! 2781: return 0;
! 2782: return 1;
! 2783: }
1.66 noro 2784: }
2785:
2786:
1.20 noro 2787: void dp_subd(DP p1,DP p2,DP *rp)
1.5 noro 2788: {
1.67 ! noro 2789: int i,n;
! 2790: DL d1,d2,d;
! 2791: MP m;
! 2792: DP s;
! 2793:
! 2794: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 2795: NEWDL(d,n); d->td = d1->td - d2->td;
! 2796: for ( i = 0; i < n; i++ )
! 2797: d->d[i] = d1->d[i]-d2->d[i];
! 2798: NEWMP(m); m->dl = d; m->c = (Obj)ONE; NEXT(m) = 0;
! 2799: MKDP(n,m,s); s->sugar = d->td;
! 2800: *rp = s;
1.7 noro 2801: }
2802:
1.20 noro 2803: void dltod(DL d,int n,DP *rp)
1.7 noro 2804: {
1.67 ! noro 2805: MP m;
! 2806: DP s;
1.7 noro 2807:
1.67 ! noro 2808: NEWMP(m); m->dl = d; m->c = (Obj)ONE; NEXT(m) = 0;
! 2809: MKDP(n,m,s); s->sugar = d->td;
! 2810: *rp = s;
1.5 noro 2811: }
2812:
1.20 noro 2813: void dp_hm(DP p,DP *rp)
1.5 noro 2814: {
1.67 ! noro 2815: MP m,mr;
1.5 noro 2816:
1.67 ! noro 2817: if ( !p )
! 2818: *rp = 0;
! 2819: else {
! 2820: m = BDY(p);
! 2821: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
! 2822: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
! 2823: }
1.5 noro 2824: }
2825:
1.35 noro 2826: void dp_ht(DP p,DP *rp)
2827: {
1.67 ! noro 2828: MP m,mr;
1.35 noro 2829:
1.67 ! noro 2830: if ( !p )
! 2831: *rp = 0;
! 2832: else {
! 2833: m = BDY(p);
! 2834: NEWMP(mr); mr->dl = m->dl; mr->c = (Obj)ONE; NEXT(mr) = 0;
! 2835: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
! 2836: }
1.35 noro 2837: }
2838:
1.66 noro 2839: void dpm_hm(DPM p,DPM *rp)
2840: {
1.67 ! noro 2841: DMM m,mr;
1.66 noro 2842:
1.67 ! noro 2843: if ( !p )
! 2844: *rp = 0;
! 2845: else {
! 2846: m = BDY(p);
! 2847: NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; NEXT(mr) = 0;
! 2848: MKDPM(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
! 2849: }
1.66 noro 2850: }
2851:
2852: void dpm_ht(DPM p,DPM *rp)
2853: {
1.67 ! noro 2854: DMM m,mr;
1.66 noro 2855:
1.67 ! noro 2856: if ( !p )
! 2857: *rp = 0;
! 2858: else {
! 2859: m = BDY(p);
! 2860: NEWDMM(mr); mr->dl = m->dl; mr->pos = m->pos; mr->c = (Obj)ONE; NEXT(mr) = 0;
! 2861: MKDPM(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
! 2862: }
1.66 noro 2863: }
2864:
2865:
1.20 noro 2866: void dp_rest(DP p,DP *rp)
1.5 noro 2867: {
1.67 ! noro 2868: MP m;
1.5 noro 2869:
1.67 ! noro 2870: m = BDY(p);
! 2871: if ( !NEXT(m) )
! 2872: *rp = 0;
! 2873: else {
! 2874: MKDP(p->nv,NEXT(m),*rp);
! 2875: if ( *rp )
! 2876: (*rp)->sugar = p->sugar;
! 2877: }
1.5 noro 2878: }
2879:
1.66 noro 2880: void dpm_rest(DPM p,DPM *rp)
2881: {
1.67 ! noro 2882: DMM m;
1.66 noro 2883:
1.67 ! noro 2884: m = BDY(p);
! 2885: if ( !NEXT(m) )
! 2886: *rp = 0;
! 2887: else {
! 2888: MKDPM(p->nv,NEXT(m),*rp);
! 2889: if ( *rp )
! 2890: (*rp)->sugar = p->sugar;
! 2891: }
1.66 noro 2892: }
2893:
1.20 noro 2894: DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl)
1.5 noro 2895: {
1.67 ! noro 2896: register int i, *d1, *d2, *d, td;
1.5 noro 2897:
1.67 ! noro 2898: if ( !dl ) NEWDL(dl,nv);
! 2899: d = dl->d, d1 = dl1->d, d2 = dl2->d;
! 2900: for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) {
! 2901: *d = *d1 > *d2 ? *d1 : *d2;
! 2902: td += MUL_WEIGHT(*d,i);
! 2903: }
! 2904: dl->td = td;
! 2905: return dl;
1.5 noro 2906: }
2907:
1.20 noro 2908: int dl_equal(int nv,DL dl1,DL dl2)
1.5 noro 2909: {
2910: register int *d1, *d2, n;
2911:
2912: if ( dl1->td != dl2->td ) return 0;
2913: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
2914: if ( *d1 != *d2 ) return 0;
2915: return 1;
2916: }
2917:
1.20 noro 2918: int dp_nt(DP p)
1.5 noro 2919: {
1.67 ! noro 2920: int i;
! 2921: MP m;
1.5 noro 2922:
1.67 ! noro 2923: if ( !p )
! 2924: return 0;
! 2925: else {
! 2926: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
! 2927: return i;
! 2928: }
1.5 noro 2929: }
2930:
1.20 noro 2931: int dp_homogeneous(DP p)
1.15 noro 2932: {
1.67 ! noro 2933: MP m;
! 2934: int d;
1.15 noro 2935:
1.67 ! noro 2936: if ( !p )
! 2937: return 1;
! 2938: else {
! 2939: m = BDY(p);
! 2940: d = m->dl->td;
! 2941: m = NEXT(m);
! 2942: for ( ; m; m = NEXT(m) ) {
! 2943: if ( m->dl->td != d )
! 2944: return 0;
! 2945: }
! 2946: return 1;
! 2947: }
1.16 noro 2948: }
2949:
1.20 noro 2950: void _print_mp(int nv,MP m)
1.16 noro 2951: {
1.67 ! noro 2952: int i;
1.16 noro 2953:
1.67 ! noro 2954: if ( !m )
! 2955: return;
! 2956: for ( ; m; m = NEXT(m) ) {
! 2957: fprintf(stderr,"%d<",ITOS(C(m)));
! 2958: for ( i = 0; i < nv; i++ ) {
! 2959: fprintf(stderr,"%d",m->dl->d[i]);
! 2960: if ( i != nv-1 )
! 2961: fprintf(stderr," ");
! 2962: }
! 2963: fprintf(stderr,">",C(m));
! 2964: }
! 2965: fprintf(stderr,"\n");
1.15 noro 2966: }
1.26 noro 2967:
2968: static int cmp_mp_nvar;
2969:
2970: int comp_mp(MP *a,MP *b)
2971: {
1.67 ! noro 2972: return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl);
1.26 noro 2973: }
2974:
2975: void dp_sort(DP p,DP *rp)
2976: {
1.67 ! noro 2977: MP t,mp,mp0;
! 2978: int i,n;
! 2979: DP r;
! 2980: MP *w;
! 2981:
! 2982: if ( !p ) {
! 2983: *rp = 0;
! 2984: return;
! 2985: }
! 2986: for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ );
! 2987: w = (MP *)ALLOCA(n*sizeof(MP));
! 2988: for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ )
! 2989: w[i] = t;
! 2990: cmp_mp_nvar = NV(p);
! 2991: qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp);
! 2992: mp0 = 0;
! 2993: for ( i = n-1; i >= 0; i-- ) {
! 2994: NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]);
! 2995: NEXT(mp) = mp0; mp0 = mp;
! 2996: }
! 2997: MKDP(p->nv,mp0,r);
! 2998: r->sugar = p->sugar;
! 2999: *rp = r;
1.26 noro 3000: }
3001:
1.32 noro 3002: DP extract_initial_term_from_dp(DP p,int *weight,int n);
3003: LIST extract_initial_term(LIST f,int *weight,int n);
3004:
3005: DP extract_initial_term_from_dp(DP p,int *weight,int n)
3006: {
1.67 ! noro 3007: int w,t,i,top;
! 3008: MP m,r0,r;
! 3009: DP dp;
! 3010:
! 3011: if ( !p ) return 0;
! 3012: top = 1;
! 3013: for ( m = BDY(p); m; m = NEXT(m) ) {
! 3014: for ( i = 0, t = 0; i < n; i++ )
! 3015: t += weight[i]*m->dl->d[i];
! 3016: if ( top || t > w ) {
! 3017: r0 = 0;
! 3018: w = t;
! 3019: top = 0;
! 3020: }
! 3021: if ( t == w ) {
! 3022: NEXTMP(r0,r);
! 3023: r->dl = m->dl;
! 3024: r->c = m->c;
! 3025: }
! 3026: }
! 3027: NEXT(r) = 0;
! 3028: MKDP(p->nv,r0,dp);
! 3029: return dp;
1.32 noro 3030: }
3031:
3032: LIST extract_initial_term(LIST f,int *weight,int n)
3033: {
1.67 ! noro 3034: NODE nd,r0,r;
! 3035: Obj p;
! 3036: LIST l;
! 3037:
! 3038: nd = BDY(f);
! 3039: for ( r0 = 0; nd; nd = NEXT(nd) ) {
! 3040: NEXTNODE(r0,r);
! 3041: p = (Obj)BDY(nd);
! 3042: BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n);
! 3043: }
! 3044: if ( r0 ) NEXT(r) = 0;
! 3045: MKLIST(l,r0);
! 3046: return l;
1.32 noro 3047: }
3048:
3049: LIST dp_initial_term(LIST f,struct order_spec *ord)
3050: {
1.67 ! noro 3051: int n,l,i;
! 3052: struct weight_or_block *worb;
! 3053: int *weight;
! 3054:
! 3055: switch ( ord->id ) {
! 3056: case 2: /* matrix order */
! 3057: /* extract the first row */
! 3058: n = ord->nv;
! 3059: weight = ord->ord.matrix.matrix[0];
! 3060: return extract_initial_term(f,weight,n);
! 3061: case 3: /* composite order */
! 3062: /* the first w_or_b */
! 3063: worb = ord->ord.composite.w_or_b;
! 3064: switch ( worb->type ) {
! 3065: case IS_DENSE_WEIGHT:
! 3066: n = worb->length;
! 3067: weight = worb->body.dense_weight;
! 3068: return extract_initial_term(f,weight,n);
! 3069: case IS_SPARSE_WEIGHT:
! 3070: n = ord->nv;
! 3071: weight = (int *)ALLOCA(n*sizeof(int));
! 3072: for ( i = 0; i < n; i++ ) weight[i] = 0;
! 3073: l = worb->length;
! 3074: for ( i = 0; i < l; i++ )
! 3075: weight[worb->body.sparse_weight[i].pos]
! 3076: = worb->body.sparse_weight[i].value;
! 3077: return extract_initial_term(f,weight,n);
! 3078: default:
! 3079: error("dp_initial_term : unsupported order");
! 3080: }
! 3081: default:
! 3082: error("dp_initial_term : unsupported order");
! 3083: }
1.32 noro 3084: }
3085:
3086: int highest_order_dp(DP p,int *weight,int n);
3087: LIST highest_order(LIST f,int *weight,int n);
3088:
3089: int highest_order_dp(DP p,int *weight,int n)
3090: {
1.67 ! noro 3091: int w,t,i,top;
! 3092: MP m;
1.32 noro 3093:
1.67 ! noro 3094: if ( !p ) return -1;
! 3095: top = 1;
! 3096: for ( m = BDY(p); m; m = NEXT(m) ) {
! 3097: for ( i = 0, t = 0; i < n; i++ )
! 3098: t += weight[i]*m->dl->d[i];
! 3099: if ( top || t > w ) {
! 3100: w = t;
! 3101: top = 0;
! 3102: }
! 3103: }
! 3104: return w;
1.32 noro 3105: }
3106:
3107: LIST highest_order(LIST f,int *weight,int n)
3108: {
1.67 ! noro 3109: int h;
! 3110: NODE nd,r0,r;
! 3111: Obj p;
! 3112: LIST l;
! 3113: Q q;
! 3114:
! 3115: nd = BDY(f);
! 3116: for ( r0 = 0; nd; nd = NEXT(nd) ) {
! 3117: NEXTNODE(r0,r);
! 3118: p = (Obj)BDY(nd);
! 3119: h = highest_order_dp((DP)p,weight,n);
! 3120: STOQ(h,q);
! 3121: BDY(r) = (pointer)q;
! 3122: }
! 3123: if ( r0 ) NEXT(r) = 0;
! 3124: MKLIST(l,r0);
! 3125: return l;
1.32 noro 3126: }
3127:
3128: LIST dp_order(LIST f,struct order_spec *ord)
3129: {
1.67 ! noro 3130: int n,l,i;
! 3131: struct weight_or_block *worb;
! 3132: int *weight;
! 3133:
! 3134: switch ( ord->id ) {
! 3135: case 2: /* matrix order */
! 3136: /* extract the first row */
! 3137: n = ord->nv;
! 3138: weight = ord->ord.matrix.matrix[0];
! 3139: return highest_order(f,weight,n);
! 3140: case 3: /* composite order */
! 3141: /* the first w_or_b */
! 3142: worb = ord->ord.composite.w_or_b;
! 3143: switch ( worb->type ) {
! 3144: case IS_DENSE_WEIGHT:
! 3145: n = worb->length;
! 3146: weight = worb->body.dense_weight;
! 3147: return highest_order(f,weight,n);
! 3148: case IS_SPARSE_WEIGHT:
! 3149: n = ord->nv;
! 3150: weight = (int *)ALLOCA(n*sizeof(int));
! 3151: for ( i = 0; i < n; i++ ) weight[i] = 0;
! 3152: l = worb->length;
! 3153: for ( i = 0; i < l; i++ )
! 3154: weight[worb->body.sparse_weight[i].pos]
! 3155: = worb->body.sparse_weight[i].value;
! 3156: return highest_order(f,weight,n);
! 3157: default:
! 3158: error("dp_initial_term : unsupported order");
! 3159: }
! 3160: default:
! 3161: error("dp_initial_term : unsupported order");
! 3162: }
1.35 noro 3163: }
3164:
3165: int dpv_ht(DPV p,DP *h)
3166: {
1.67 ! noro 3167: int len,max,maxi,i,t;
! 3168: DP *e;
! 3169: MP m,mr;
! 3170:
! 3171: len = p->len;
! 3172: e = p->body;
! 3173: max = -1;
! 3174: maxi = -1;
! 3175: for ( i = 0; i < len; i++ )
! 3176: if ( e[i] && (t = BDY(e[i])->dl->td) > max ) {
! 3177: max = t;
! 3178: maxi = i;
! 3179: }
! 3180: if ( max < 0 ) {
! 3181: *h = 0;
! 3182: return -1;
! 3183: } else {
! 3184: m = BDY(e[maxi]);
! 3185: NEWMP(mr); mr->dl = m->dl; mr->c = (Obj)ONE; NEXT(mr) = 0;
! 3186: MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td; /* XXX */
! 3187: return maxi;
! 3188: }
1.32 noro 3189: }
1.42 noro 3190:
3191: /* return 1 if 0 <_w1 v && v <_w2 0 */
3192:
3193: int in_c12(int n,int *v,int row1,int **w1,int row2, int **w2)
3194: {
1.67 ! noro 3195: int t1,t2;
1.42 noro 3196:
1.67 ! noro 3197: t1 = compare_zero(n,v,row1,w1);
! 3198: t2 = compare_zero(n,v,row2,w2);
! 3199: if ( t1 > 0 && t2 < 0 ) return 1;
! 3200: else return 0;
1.42 noro 3201: }
3202:
3203: /* 0 < u => 1, 0 > u => -1 */
3204:
3205: int compare_zero(int n,int *u,int row,int **w)
3206: {
1.67 ! noro 3207: int i,j,t;
! 3208: int *wi;
1.42 noro 3209:
1.67 ! noro 3210: for ( i = 0; i < row; i++ ) {
! 3211: wi = w[i];
! 3212: for ( j = 0, t = 0; j < n; j++ ) t += u[j]*wi[j];
! 3213: if ( t > 0 ) return 1;
! 3214: else if ( t < 0 ) return -1;
! 3215: }
! 3216: return 0;
1.42 noro 3217: }
3218:
3219: /* functions for generic groebner walk */
3220: /* u=0 means u=-infty */
3221:
3222: int compare_facet_preorder(int n,int *u,int *v,
1.67 ! noro 3223: int row1,int **w1,int row2,int **w2)
1.42 noro 3224: {
1.67 ! noro 3225: int i,j,s,t,tu,tv;
! 3226: int *w2i,*uv;
1.42 noro 3227:
1.67 ! noro 3228: if ( !u ) return 1;
! 3229: uv = W_ALLOC(n);
! 3230: for ( i = 0; i < row2; i++ ) {
! 3231: w2i = w2[i];
! 3232: for ( j = 0, tu = tv = 0; j < n; j++ )
! 3233: if ( s = w2i[j] ) {
! 3234: tu += s*u[j]; tv += s*v[j];
! 3235: }
! 3236: for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu;
! 3237: t = compare_zero(n,uv,row1,w1);
! 3238: if ( t > 0 ) return 1;
! 3239: else if ( t < 0 ) return 0;
! 3240: }
! 3241: return 1;
1.42 noro 3242: }
3243:
1.48 noro 3244: Q inner_product_with_small_vector(VECT w,int *v)
3245: {
1.67 ! noro 3246: int n,i;
! 3247: Q q,s,t,u;
1.48 noro 3248:
1.67 ! noro 3249: n = w->len;
! 3250: s = 0;
! 3251: for ( i = 0; i < n; i++ ) {
! 3252: STOQ(v[i],q); mulq((Q)w->body[i],q,&t); addq(t,s,&u); s = u;
! 3253: }
! 3254: return s;
1.48 noro 3255: }
3256:
3257: Q compute_last_t(NODE g,NODE gh,Q t,VECT w1,VECT w2,NODE *homo,VECT *wp)
3258: {
1.67 ! noro 3259: int n,i;
! 3260: int *wt;
! 3261: Q last,d1,d2,dn,nm,s,t1;
! 3262: VECT wd,wt1,wt2,w;
! 3263: NODE tg,tgh;
! 3264: MP f;
! 3265: int *h;
! 3266: NODE r0,r;
! 3267: MP m0,m;
! 3268: DP d;
! 3269:
! 3270: n = w1->len;
! 3271: wt = W_ALLOC(n);
! 3272: last = ONE;
! 3273: /* t1 = 1-t */
! 3274: for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
! 3275: f = BDY((DP)BDY(tg));
! 3276: h = BDY((DP)BDY(tgh))->dl->d;
! 3277: for ( ; f; f = NEXT(f) ) {
! 3278: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
! 3279: for ( i = 0; i < n && !wt[i]; i++ );
! 3280: if ( i == n ) continue;
! 3281: d1 = inner_product_with_small_vector(w1,wt);
! 3282: d2 = inner_product_with_small_vector(w2,wt);
! 3283: nm = d1; subq(d1,d2,&dn);
! 3284: /* if d1=d2 then nothing happens */
! 3285: if ( !dn ) continue;
! 3286: /* s satisfies ds = 0*/
! 3287: divq(nm,dn,&s);
! 3288:
! 3289: if ( cmpq(s,t) > 0 && cmpq(s,last) < 0 )
! 3290: last = s;
! 3291: else if ( !cmpq(s,t) ) {
! 3292: if ( cmpq(d2,0) < 0 ) {
! 3293: last = t;
! 3294: break;
! 3295: }
! 3296: }
! 3297: }
! 3298: }
! 3299: if ( !last ) {
! 3300: dn = ONE; nm = 0;
! 3301: } else {
! 3302: NTOQ(NM(last),1,nm);
! 3303: if ( INT(last) ) dn = ONE;
! 3304: else {
! 3305: NTOQ(DN(last),1,dn);
! 3306: }
! 3307: }
! 3308: /* (1-n/d)*w1+n/d*w2 -> w=(d-n)*w1+n*w2 */
! 3309: subq(dn,nm,&t1); mulvect(CO,(Obj)w1,(Obj)t1,(Obj *)&wt1);
! 3310: mulvect(CO,(Obj)w2,(Obj)nm,(Obj *)&wt2); addvect(CO,wt1,wt2,&w);
! 3311:
! 3312: r0 = 0;
! 3313: for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
! 3314: f = BDY((DP)BDY(tg));
! 3315: h = BDY((DP)BDY(tgh))->dl->d;
! 3316: for ( m0 = 0; f; f = NEXT(f) ) {
! 3317: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
! 3318: for ( i = 0; i < n && !wt[i]; i++ );
! 3319: if ( !inner_product_with_small_vector(w,wt) ) {
! 3320: NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
! 3321: }
! 3322: }
! 3323: NEXT(m) = 0;
! 3324: MKDP(((DP)BDY(tg))->nv,m0,d); d->sugar = ((DP)BDY(tg))->sugar;
! 3325: NEXTNODE(r0,r); BDY(r) = (pointer)d;
! 3326: }
! 3327: NEXT(r) = 0;
! 3328: *homo = r0;
! 3329: *wp = w;
! 3330: return last;
1.48 noro 3331: }
3332:
1.42 noro 3333: /* return 0 if last_w = infty */
3334:
3335: NODE compute_last_w(NODE g,NODE gh,int n,int **w,
1.67 ! noro 3336: int row1,int **w1,int row2,int **w2)
1.42 noro 3337: {
1.67 ! noro 3338: DP d;
! 3339: MP f,m0,m;
! 3340: int *wt,*v,*h;
! 3341: NODE t,s,n0,tn,n1,r0,r;
! 3342: int i;
! 3343:
! 3344: wt = W_ALLOC(n);
! 3345: n0 = 0;
! 3346: for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
! 3347: f = BDY((DP)BDY(t));
! 3348: h = BDY((DP)BDY(s))->dl->d;
! 3349: for ( ; f; f = NEXT(f) ) {
! 3350: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
! 3351: for ( i = 0; i < n && !wt[i]; i++ );
! 3352: if ( i == n ) continue;
! 3353:
! 3354: if ( in_c12(n,wt,row1,w1,row2,w2) &&
! 3355: compare_facet_preorder(n,*w,wt,row1,w1,row2,w2) ) {
! 3356: v = (int *)MALLOC_ATOMIC(n*sizeof(int));
! 3357: for ( i = 0; i < n; i++ ) v[i] = wt[i];
! 3358: MKNODE(n1,v,n0); n0 = n1;
! 3359: }
! 3360: }
! 3361: }
! 3362: if ( !n0 ) return 0;
! 3363: for ( t = n0; t; t = NEXT(t) ) {
! 3364: v = (int *)BDY(t);
! 3365: for ( s = n0; s; s = NEXT(s) )
! 3366: if ( !compare_facet_preorder(n,v,(int *)BDY(s),row1,w1,row2,w2) )
! 3367: break;
! 3368: if ( !s ) {
! 3369: *w = v;
! 3370: break;
! 3371: }
! 3372: }
! 3373: if ( !t )
! 3374: error("compute_last_w : cannot happen");
! 3375: r0 = 0;
! 3376: for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
! 3377: f = BDY((DP)BDY(t));
! 3378: h = BDY((DP)BDY(s))->dl->d;
! 3379: for ( m0 = 0; f; f = NEXT(f) ) {
! 3380: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
! 3381: for ( i = 0; i < n && !wt[i]; i++ );
! 3382: if ( i == n ||
! 3383: (compare_facet_preorder(n,wt,*w,row1,w1,row2,w2)
! 3384: && compare_facet_preorder(n,*w,wt,row1,w1,row2,w2)) ) {
! 3385: NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
! 3386: }
! 3387: }
! 3388: NEXT(m) = 0;
! 3389: MKDP(((DP)BDY(t))->nv,m0,d); d->sugar = ((DP)BDY(t))->sugar;
! 3390: NEXTNODE(r0,r); BDY(r) = (pointer)d;
! 3391: }
! 3392: NEXT(r) = 0;
! 3393: return r0;
1.42 noro 3394: }
1.44 noro 3395:
3396: /* compute a sufficient set of d(f)=u-v */
3397:
3398: NODE compute_essential_df(DP *g,DP *gh,int ng)
3399: {
1.67 ! noro 3400: int nv,i,j,k,t,lj;
! 3401: NODE r,r1,ri,rt,r0;
! 3402: MP m;
! 3403: MP *mj;
! 3404: DL di,hj,dl,dlt;
! 3405: int *d,*dt;
! 3406: LIST l;
! 3407: Q q;
! 3408:
! 3409: nv = g[0]->nv;
! 3410: r = 0;
! 3411: for ( j = 0; j < ng; j++ ) {
! 3412: for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ );
! 3413: mj = (MP *)ALLOCA(lj*sizeof(MP));
! 3414: for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ )
! 3415: mj[k] = m;
! 3416: for ( i = 0; i < lj; i++ ) {
! 3417: for ( di = mj[i]->dl, k = i+1; k < lj; k++ )
! 3418: if ( _dl_redble(di,mj[k]->dl,nv) ) break;
! 3419: if ( k < lj ) mj[i] = 0;
! 3420: }
! 3421: hj = BDY(gh[j])->dl;
! 3422: _NEWDL(dl,nv); d = dl->d;
! 3423: r0 = r;
! 3424: for ( i = 0; i < lj; i++ ) {
! 3425: if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) {
! 3426: for ( k = 0, t = 0; k < nv; k++ ) {
! 3427: d[k] = hj->d[k]-di->d[k];
! 3428: t += d[k];
! 3429: }
! 3430: dl->td = t;
1.45 noro 3431: #if 1
1.67 ! noro 3432: for ( rt = r0; rt; rt = NEXT(rt) ) {
! 3433: dlt = (DL)BDY(rt);
! 3434: if ( dlt->td != dl->td ) continue;
! 3435: for ( dt = dlt->d, k = 0; k < nv; k++ )
! 3436: if ( d[k] != dt[k] ) break;
! 3437: if ( k == nv ) break;
! 3438: }
1.45 noro 3439: #else
1.67 ! noro 3440: rt = 0;
1.45 noro 3441: #endif
1.67 ! noro 3442: if ( !rt ) {
! 3443: MKNODE(r1,dl,r); r = r1;
! 3444: _NEWDL(dl,nv); d = dl->d;
! 3445: }
! 3446: }
! 3447: }
! 3448: }
! 3449: for ( rt = r; rt; rt = NEXT(rt) ) {
! 3450: dl = (DL)BDY(rt); d = dl->d;
! 3451: ri = 0;
! 3452: for ( k = nv-1; k >= 0; k-- ) {
! 3453: STOQ(d[k],q);
! 3454: MKNODE(r1,q,ri); ri = r1;
! 3455: }
! 3456: MKNODE(r1,0,ri); MKLIST(l,r1);
! 3457: BDY(rt) = (pointer)l;
! 3458: }
! 3459: return r;
1.44 noro 3460: }
1.57 noro 3461:
3462: int comp_bits_divisible(int *a,int *b,int n)
3463: {
1.67 ! noro 3464: int bpi,i,wi,bi;
1.57 noro 3465:
1.67 ! noro 3466: bpi = (sizeof(int)/sizeof(char))*8;
! 3467: for ( i = 0; i < n; i++ ) {
! 3468: wi = i/bpi; bi = i%bpi;
! 3469: if ( !(a[wi]&(1<<bi)) && (b[wi]&(1<<bi)) ) return 0;
! 3470: }
! 3471: return 1;
1.57 noro 3472: }
3473:
3474: int comp_bits_lex(int *a,int *b,int n)
3475: {
1.67 ! noro 3476: int bpi,i,wi,ba,bb,bi;
1.57 noro 3477:
1.67 ! noro 3478: bpi = (sizeof(int)/sizeof(char))*8;
! 3479: for ( i = 0; i < n; i++ ) {
! 3480: wi = i/bpi; bi = i%bpi;
! 3481: ba = (a[wi]&(1<<bi))?1:0;
! 3482: bb = (b[wi]&(1<<bi))?1:0;
! 3483: if ( ba > bb ) return 1;
! 3484: else if ( ba < bb ) return -1;
! 3485: }
! 3486: return 0;
1.57 noro 3487: }
3488:
3489: NODE mono_raddec(NODE ideal)
3490: {
1.67 ! noro 3491: DP p;
! 3492: int nv,w,i,bpi,di,c,len;
! 3493: int *d,*s,*u,*new;
! 3494: NODE t,t1,v,r,rem,prev;
! 3495:
! 3496: if( !ideal ) return 0;
! 3497: p = (DP)BDY(ideal);
! 3498: nv = NV(p);
! 3499: bpi = (sizeof(int)/sizeof(char))*8;
! 3500: w = (nv+(bpi-1))/bpi;
! 3501: d = p->body->dl->d;
! 3502: if ( !NEXT(ideal) ) {
! 3503: for ( t = 0, i = nv-1; i >= 0; i-- ) {
! 3504: if ( d[i] ) {
! 3505: s = (int *)CALLOC(w,sizeof(int));
! 3506: s[i/bpi] |= 1<<(i%bpi);
! 3507: MKNODE(t1,s,t);
! 3508: t = t1;
! 3509: }
! 3510: }
! 3511: return t;
! 3512: }
! 3513: rem = mono_raddec(NEXT(ideal));
! 3514: r = 0;
! 3515: len = w*sizeof(int);
! 3516: u = (int *)CALLOC(w,sizeof(int));
! 3517: for ( i = nv-1; i >= 0; i-- ) {
! 3518: if ( d[i] ) {
! 3519: for ( t = rem; t; t = NEXT(t) ) {
! 3520: bcopy((char *)BDY(t),(char *)u,len);
! 3521: u[i/bpi] |= 1<<(i%bpi);
! 3522: for ( v = r; v; v = NEXT(v) ) {
! 3523: if ( comp_bits_divisible(u,(int *)BDY(v),nv) ) break;
! 3524: }
! 3525: if ( v ) continue;
! 3526: for ( v = r, prev = 0; v; v = NEXT(v) ) {
! 3527: if ( comp_bits_divisible((int *)BDY(v),u,nv) ) {
! 3528: if ( prev ) NEXT(prev) = NEXT(v);
! 3529: else r = NEXT(r);
! 3530: } else prev =v;
! 3531: }
! 3532: for ( v = r, prev = 0; v; prev = v, v = NEXT(v) ) {
! 3533: if ( comp_bits_lex(u,(int *)BDY(v),nv) < 0 ) break;
! 3534: }
! 3535: new = (int *)CALLOC(w,sizeof(int));
! 3536: bcopy((char *)u,(char *)new,len);
! 3537: MKNODE(t1,new,v);
! 3538: if ( prev ) NEXT(prev) = t1;
! 3539: else r = t1;
! 3540: }
! 3541: }
! 3542: }
! 3543: return r;
1.57 noro 3544: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dp-supp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin