![[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|
Return to dp-supp.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
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.37 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.36 2004/05/14 09:20:56 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.5 noro 56: #define HMAG(p) (p_mag(BDY(p)->c))
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);
! 68:
1.7 noro 69: /*
70: * content reduction
71: *
72: */
73:
1.20 noro 74: void dp_ptozp(DP p,DP *rp)
1.7 noro 75: {
76: MP m,mr,mr0;
77: int i,n;
78: Q *w;
79: Q dvr;
80: P t;
81:
82: if ( !p )
83: *rp = 0;
84: else {
85: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
86: w = (Q *)ALLOCA(n*sizeof(Q));
87: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
88: if ( NUM(m->c) )
89: w[i] = (Q)m->c;
90: else
91: ptozp(m->c,1,&w[i],&t);
92: sortbynm(w,n);
93: qltozl(w,n,&dvr);
94: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
95: NEXTMP(mr0,mr); divsp(CO,m->c,(P)dvr,&mr->c); mr->dl = m->dl;
96: }
97: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
98: }
99: }
100:
1.20 noro 101: void dp_ptozp2(DP p0,DP p1,DP *hp,DP *rp)
1.7 noro 102: {
103: DP t,s,h,r;
104: MP m,mr,mr0,m0;
105:
106: addd(CO,p0,p1,&t); dp_ptozp(t,&s);
107: if ( !p0 ) {
108: h = 0; r = s;
109: } else if ( !p1 ) {
110: h = s; r = 0;
111: } else {
112: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
113: m = NEXT(m), m0 = NEXT(m0) ) {
114: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
115: }
116: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
117: }
118: if ( h )
119: h->sugar = p0->sugar;
120: if ( r )
121: r->sugar = p1->sugar;
122: *hp = h; *rp = r;
123: }
1.1 noro 124:
1.20 noro 125: void dp_idiv(DP p,Q c,DP *rp)
1.1 noro 126: {
127: Q t;
128: N nm,q;
129: int sgn,s;
130: MP mr0,m,mr;
131:
132: if ( !p )
133: *rp = 0;
134: else if ( MUNIQ((Q)c) )
135: *rp = p;
136: else if ( MUNIQ((Q)c) )
137: chsgnd(p,rp);
138: else {
139: nm = NM(c); sgn = SGN(c);
140: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
141: NEXTMP(mr0,mr);
142:
143: divsn(NM((Q)(m->c)),nm,&q);
144: s = sgn*SGN((Q)(m->c));
145: NTOQ(q,s,t);
146: mr->c = (P)t;
147: mr->dl = m->dl;
148: }
149: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
150: if ( *rp )
151: (*rp)->sugar = p->sugar;
152: }
153: }
154:
1.20 noro 155: void dp_mbase(NODE hlist,NODE *mbase)
1.1 noro 156: {
157: DL *dl;
158: DL d;
159: int i,j,n,nvar,td;
160:
161: n = length(hlist); nvar = ((DP)BDY(hlist))->nv;
162: dl = (DL *)MALLOC(n*sizeof(DL));
163: for ( i = 0; i < n; i++, hlist = NEXT(hlist) )
164: dl[i] = BDY((DP)BDY(hlist))->dl;
165: NEWDL(d,nvar); *mbase = 0;
166: while ( 1 ) {
167: insert_to_node(d,mbase,nvar);
168: for ( i = nvar-1; i >= 0; ) {
1.21 noro 169: d->d[i]++;
170: d->td += MUL_WEIGHT(1,i);
1.1 noro 171: for ( j = 0; j < n; j++ ) {
172: if ( _dl_redble(dl[j],d,nvar) )
173: break;
174: }
175: if ( j < n ) {
176: for ( j = nvar-1; j >= i; j-- )
177: d->d[j] = 0;
178: for ( j = 0, td = 0; j < i; j++ )
1.21 noro 179: td += MUL_WEIGHT(d->d[j],j);
1.1 noro 180: d->td = td;
181: i--;
182: } else
183: break;
184: }
185: if ( i < 0 )
186: break;
187: }
188: }
189:
1.20 noro 190: int _dl_redble(DL d1,DL d2,int nvar)
1.1 noro 191: {
192: int i;
193:
194: if ( d1->td > d2->td )
195: return 0;
196: for ( i = 0; i < nvar; i++ )
197: if ( d1->d[i] > d2->d[i] )
198: break;
199: if ( i < nvar )
200: return 0;
201: else
202: return 1;
203: }
204:
1.20 noro 205: void insert_to_node(DL d,NODE *n,int nvar)
1.1 noro 206: {
207: DL d1;
208: MP m;
209: DP dp;
210: NODE n0,n1,n2;
211:
212: NEWDL(d1,nvar); d1->td = d->td;
213: bcopy((char *)d->d,(char *)d1->d,nvar*sizeof(int));
214: NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0;
215: MKDP(nvar,m,dp); dp->sugar = d->td;
216: if ( !(*n) ) {
217: MKNODE(n1,dp,0); *n = n1;
218: } else {
219: for ( n1 = *n, n0 = 0; n1; n0 = n1, n1 = NEXT(n1) )
220: if ( (*cmpdl)(nvar,d,BDY((DP)BDY(n1))->dl) > 0 ) {
221: MKNODE(n2,dp,n1);
222: if ( !n0 )
223: *n = n2;
224: else
225: NEXT(n0) = n2;
226: break;
227: }
228: if ( !n1 ) {
229: MKNODE(n2,dp,0); NEXT(n0) = n2;
230: }
231: }
232: }
233:
1.20 noro 234: void dp_vtod(Q *c,DP p,DP *rp)
1.1 noro 235: {
236: MP mr0,m,mr;
237: int i;
238:
239: if ( !p )
240: *rp = 0;
241: else {
242: for ( mr0 = 0, m = BDY(p), i = 0; m; m = NEXT(m), i++ ) {
243: NEXTMP(mr0,mr); mr->c = (P)c[i]; mr->dl = m->dl;
244: }
245: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
246: (*rp)->sugar = p->sugar;
247: }
248: }
249:
1.8 noro 250: extern int mpi_mag;
251: extern int PCoeffs;
252:
1.20 noro 253: void dp_ptozp_d(DP p,DP *rp)
1.1 noro 254: {
255: int i,j,k,l,n,nsep;
256: MP m;
257: NODE tn,n0,n1,n2,n3;
258: struct oVECT v;
259: VECT c,cs;
260: VECT qi,ri;
261: LIST *qr;
262: Obj dmy;
263: Q d0,d1,gcd,a,u,u1;
264: Q *q,*r;
265: STRING iqr_v;
266: pointer *b;
267: N qn,gn;
268: double get_rtime();
269: int blen;
1.8 noro 270: NODE dist;
271: int ndist;
1.1 noro 272: double t0;
273: double t_e,t_d,t_d1,t_c;
1.8 noro 274: extern int DP_NFStat;
275: extern LIST Dist;
1.20 noro 276: void Pox_rpc();
277: void Pox_pop_local();
1.1 noro 278:
279: if ( !p )
280: *rp = 0;
281: else {
1.8 noro 282: if ( PCoeffs ) {
283: dp_ptozp(p,rp); return;
284: }
1.9 noro 285: if ( !Dist || p_mag(BDY(p)->c) <= mpi_mag ) {
1.8 noro 286: dist = 0; ndist = 0;
287: if ( DP_NFStat ) fprintf(asir_out,"L");
288: } else {
289: dist = BDY(Dist); ndist = length(dist);
290: if ( DP_NFStat ) fprintf(asir_out,"D");
291: }
1.1 noro 292: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
293: nsep = ndist + 1;
294: if ( n <= nsep ) {
295: dp_ptozp(p,rp); return;
296: }
297: t0 = get_rtime();
298: dp_dtov(p,&c);
299: igcdv_estimate(c,&d0);
300: t_e = get_rtime()-t0;
301: t0 = get_rtime();
302: dp_dtov(p,&c);
303: sepvect(c,nsep,&cs);
304: MKSTR(iqr_v,"iqr");
305: qr = (LIST *)CALLOC(nsep,sizeof(LIST));
306: q = (Q *)CALLOC(n,sizeof(Q));
307: r = (Q *)CALLOC(n,sizeof(Q));
308: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
309: MKNODE(n3,d0,0); MKNODE(n2,b[i],n3);
310: MKNODE(n1,iqr_v,n2); MKNODE(n0,BDY(tn),n1);
311: Pox_rpc(n0,&dmy);
312: }
313: iqrv(b[i],d0,&qr[i]);
314: dp_dtov(p,&c);
315: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
316: Pox_pop_local(tn,&qr[i]);
317: if ( OID(qr[i]) == O_ERR ) {
318: printexpr(CO,(Obj)qr[i]);
319: error("dp_ptozp_d : aborted");
320: }
321: }
322: t_d = get_rtime()-t0;
323: t_d1 = t_d/n;
324: t0 = get_rtime();
325: for ( i = j = 0; i < nsep; i++ ) {
326: tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
327: for ( k = 0, l = qi->len; k < l; k++, j++ ) {
328: q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
329: }
330: }
331: v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
332: if ( d1 ) {
333: gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
334: divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
335: for ( i = 0; i < n; i++ ) {
336: mulq(a,q[i],&u);
337: if ( r[i] ) {
338: divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
339: addq(u,u1,&q[i]);
340: } else
341: q[i] = u;
342: }
343: } else
344: gcd = d0;
345: dp_vtod(q,p,rp);
346: t_c = get_rtime()-t0;
347: blen=p_mag((P)gcd);
348: pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
349: if ( 0 )
350: fprintf(stderr,"(%d,%d)",p_mag((P)d0)-blen,blen);
351: }
352: }
353:
1.20 noro 354: void dp_ptozp2_d(DP p0,DP p1,DP *hp,DP *rp)
1.1 noro 355: {
356: DP t,s,h,r;
357: MP m,mr,mr0,m0;
358:
1.8 noro 359: addd(CO,p0,p1,&t); dp_ptozp_d(t,&s);
1.1 noro 360: if ( !p0 ) {
361: h = 0; r = s;
362: } else if ( !p1 ) {
363: h = s; r = 0;
364: } else {
365: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
366: m = NEXT(m), m0 = NEXT(m0) ) {
367: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
368: }
369: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
370: }
371: if ( h )
372: h->sugar = p0->sugar;
373: if ( r )
374: r->sugar = p1->sugar;
375: *hp = h; *rp = r;
1.5 noro 376: }
377:
1.22 noro 378: int have_sf_coef(P p)
379: {
380: DCP dc;
381:
382: if ( !p )
383: return 0;
384: else if ( NUM(p) )
385: return NID((Num)p) == N_GFS ? 1 : 0;
386: else {
387: for ( dc = DC(p); dc; dc = NEXT(dc) )
388: if ( have_sf_coef(COEF(dc)) )
389: return 1;
390: return 0;
391: }
392: }
393:
1.25 noro 394: void head_coef(P p,Num *c)
395: {
396: if ( !p )
397: *c = 0;
398: else if ( NUM(p) )
399: *c = (Num)p;
400: else
401: head_coef(COEF(DC(p)),c);
402: }
403:
404: void dp_monic_sf(DP p,DP *rp)
405: {
406: Num c;
407:
408: if ( !p )
409: *rp = 0;
410: else {
411: head_coef(BDY(p)->c,&c);
412: divsdc(CO,p,(P)c,rp);
413: }
414: }
415:
1.20 noro 416: void dp_prim(DP p,DP *rp)
1.5 noro 417: {
1.7 noro 418: P t,g;
419: DP p1;
420: MP m,mr,mr0;
421: int i,n;
422: P *w;
423: Q *c;
424: Q dvr;
1.5 noro 425:
1.7 noro 426: if ( !p )
427: *rp = 0;
1.23 noro 428: else if ( dp_fcoeffs == N_GFS ) {
429: for ( m = BDY(p); m; m = NEXT(m) )
1.22 noro 430: if ( OID(m->c) == O_N ) {
431: /* GCD of coeffs = 1 */
1.25 noro 432: dp_monic_sf(p,rp);
1.22 noro 433: return;
1.23 noro 434: } else break;
435: /* compute GCD over the finite fieid */
436: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
437: w = (P *)ALLOCA(n*sizeof(P));
438: for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ )
439: w[i] = m->c;
440: gcdsf(CO,w,n,&g);
441: if ( NUM(g) )
1.25 noro 442: dp_monic_sf(p,rp);
1.23 noro 443: else {
444: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
445: NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
1.22 noro 446: }
1.25 noro 447: NEXT(mr) = 0; MKDP(p->nv,mr0,p1); p1->sugar = p->sugar;
448: dp_monic_sf(p1,rp);
1.22 noro 449: }
1.23 noro 450: return;
451: } else if ( dp_fcoeffs )
1.7 noro 452: *rp = p;
1.23 noro 453: else if ( NoGCD )
1.7 noro 454: dp_ptozp(p,rp);
455: else {
456: dp_ptozp(p,&p1); p = p1;
457: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
458: if ( n == 1 ) {
459: m = BDY(p);
460: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
461: MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
462: return;
463: }
464: w = (P *)ALLOCA(n*sizeof(P));
465: c = (Q *)ALLOCA(n*sizeof(Q));
466: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
467: if ( NUM(m->c) ) {
468: c[i] = (Q)m->c; w[i] = (P)ONE;
469: } else
470: ptozp(m->c,1,&c[i],&w[i]);
471: qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
472: if ( NUM(g) )
473: *rp = p;
474: else {
475: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
476: NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
477: }
478: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 479: }
1.7 noro 480: }
1.5 noro 481: }
482:
1.20 noro 483: void heu_nezgcdnpz(VL vl,P *pl,int m,P *pr)
1.5 noro 484: {
485: int i,r;
486: P gcd,t,s1,s2,u;
487: Q rq;
488:
489: while ( 1 ) {
490: for ( i = 0, s1 = 0; i < m; i++ ) {
491: r = random(); UTOQ(r,rq);
492: mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
493: }
494: for ( i = 0, s2 = 0; i < m; i++ ) {
495: r = random(); UTOQ(r,rq);
496: mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
497: }
498: ezgcdp(vl,s1,s2,&gcd);
499: for ( i = 0; i < m; i++ ) {
500: if ( !divtpz(vl,pl[i],gcd,&t) )
501: break;
502: }
503: if ( i == m )
504: break;
505: }
506: *pr = gcd;
507: }
508:
1.20 noro 509: void dp_prim_mod(DP p,int mod,DP *rp)
1.5 noro 510: {
511: P t,g;
512: MP m,mr,mr0;
513:
514: if ( !p )
515: *rp = 0;
516: else if ( NoGCD )
517: *rp = p;
518: else {
519: for ( m = BDY(p), g = m->c, m = NEXT(m); m; m = NEXT(m) ) {
520: gcdprsmp(CO,mod,g,m->c,&t); g = t;
521: }
522: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
523: NEXTMP(mr0,mr); divsmp(CO,mod,m->c,g,&mr->c); mr->dl = m->dl;
524: }
525: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
526: }
527: }
528:
1.20 noro 529: void dp_cont(DP p,Q *rp)
1.5 noro 530: {
1.7 noro 531: VECT v;
1.5 noro 532:
1.7 noro 533: dp_dtov(p,&v); igcdv(v,rp);
1.5 noro 534: }
535:
1.20 noro 536: void dp_dtov(DP dp,VECT *rp)
1.5 noro 537: {
1.7 noro 538: MP m,t;
539: int i,n;
540: VECT v;
541: pointer *p;
1.5 noro 542:
1.7 noro 543: m = BDY(dp);
544: for ( t = m, n = 0; t; t = NEXT(t), n++ );
545: MKVECT(v,n);
546: for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
547: p[i] = (pointer)(t->c);
548: *rp = v;
1.5 noro 549: }
550:
1.7 noro 551: /*
552: * s-poly computation
553: *
554: */
1.5 noro 555:
1.20 noro 556: void dp_sp(DP p1,DP p2,DP *rp)
1.5 noro 557: {
1.7 noro 558: int i,n,td;
559: int *w;
560: DL d1,d2,d;
561: MP m;
562: DP t,s1,s2,u;
563: Q c,c1,c2;
564: N gn,tn;
1.5 noro 565:
1.7 noro 566: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
567: w = (int *)ALLOCA(n*sizeof(int));
568: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 569: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.5 noro 570: }
1.7 noro 571:
572: NEWDL(d,n); d->td = td - d1->td;
573: for ( i = 0; i < n; i++ )
574: d->d[i] = w[i] - d1->d[i];
575: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
576: if ( INT(c1) && INT(c2) ) {
577: gcdn(NM(c1),NM(c2),&gn);
578: if ( !UNIN(gn) ) {
579: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
580: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1.5 noro 581: }
582: }
1.7 noro 583:
584: NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
585: MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
586:
587: NEWDL(d,n); d->td = td - d2->td;
588: for ( i = 0; i < n; i++ )
589: d->d[i] = w[i] - d2->d[i];
590: NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0;
591: MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
592:
593: subd(CO,t,u,rp);
1.14 noro 594: if ( GenTrace ) {
595: LIST hist;
596: NODE node;
597:
598: node = mknode(4,ONE,0,s1,ONE);
599: MKLIST(hist,node);
600: MKNODE(TraceList,hist,0);
601:
602: node = mknode(4,ONE,0,0,ONE);
603: chsgnd(s2,(DP *)&ARG2(node));
604: MKLIST(hist,node);
605: MKNODE(node,hist,TraceList); TraceList = node;
606: }
607: }
608:
1.20 noro 609: void _dp_sp_dup(DP p1,DP p2,DP *rp)
1.14 noro 610: {
611: int i,n,td;
612: int *w;
613: DL d1,d2,d;
614: MP m;
615: DP t,s1,s2,u;
616: Q c,c1,c2;
617: N gn,tn;
618:
619: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
620: w = (int *)ALLOCA(n*sizeof(int));
621: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 622: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.14 noro 623: }
624:
625: _NEWDL(d,n); d->td = td - d1->td;
626: for ( i = 0; i < n; i++ )
627: d->d[i] = w[i] - d1->d[i];
628: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
629: if ( INT(c1) && INT(c2) ) {
630: gcdn(NM(c1),NM(c2),&gn);
631: if ( !UNIN(gn) ) {
632: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
633: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
634: }
635: }
636:
637: _NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
638: _MKDP(n,m,s1); s1->sugar = d->td; _muld_dup(CO,s1,p1,&t); _free_dp(s1);
639:
640: _NEWDL(d,n); d->td = td - d2->td;
641: for ( i = 0; i < n; i++ )
642: d->d[i] = w[i] - d2->d[i];
643: _NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0;
644: _MKDP(n,m,s2); s2->sugar = d->td; _muld_dup(CO,s2,p2,&u); _free_dp(s2);
645:
646: _addd_destructive(CO,t,u,rp);
1.7 noro 647: if ( GenTrace ) {
648: LIST hist;
649: NODE node;
650:
651: node = mknode(4,ONE,0,s1,ONE);
652: MKLIST(hist,node);
653: MKNODE(TraceList,hist,0);
654:
655: node = mknode(4,ONE,0,0,ONE);
656: chsgnd(s2,(DP *)&ARG2(node));
657: MKLIST(hist,node);
658: MKNODE(node,hist,TraceList); TraceList = node;
659: }
660: }
661:
1.20 noro 662: void dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7 noro 663: {
664: int i,n,td;
665: int *w;
666: DL d1,d2,d;
667: MP m;
668: DP t,s,u;
669:
670: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
671: w = (int *)ALLOCA(n*sizeof(int));
672: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 673: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7 noro 674: }
1.18 noro 675: NEWDL_NOINIT(d,n); d->td = td - d1->td;
1.7 noro 676: for ( i = 0; i < n; i++ )
677: d->d[i] = w[i] - d1->d[i];
678: NEWMP(m); m->dl = d; m->c = (P)BDY(p2)->c; NEXT(m) = 0;
679: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
1.18 noro 680: NEWDL_NOINIT(d,n); d->td = td - d2->td;
1.7 noro 681: for ( i = 0; i < n; i++ )
682: d->d[i] = w[i] - d2->d[i];
683: NEWMP(m); m->dl = d; m->c = (P)BDY(p1)->c; NEXT(m) = 0;
684: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
685: submd(CO,mod,t,u,rp);
686: }
687:
1.20 noro 688: void _dp_sp_mod_dup(DP p1,DP p2,int mod,DP *rp)
1.7 noro 689: {
690: int i,n,td;
691: int *w;
692: DL d1,d2,d;
693: MP m;
694: DP t,s,u;
695:
696: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
697: w = (int *)ALLOCA(n*sizeof(int));
698: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 699: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7 noro 700: }
701: _NEWDL(d,n); d->td = td - d1->td;
702: for ( i = 0; i < n; i++ )
703: d->d[i] = w[i] - d1->d[i];
704: _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
705: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
706: _NEWDL(d,n); d->td = td - d2->td;
707: for ( i = 0; i < n; i++ )
708: d->d[i] = w[i] - d2->d[i];
709: _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
710: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
711: _addmd_destructive(mod,t,u,rp);
712: }
713:
1.20 noro 714: void _dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7 noro 715: {
716: int i,n,td;
717: int *w;
718: DL d1,d2,d;
719: MP m;
720: DP t,s,u;
721:
722: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
723: w = (int *)ALLOCA(n*sizeof(int));
724: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 725: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7 noro 726: }
727: NEWDL(d,n); d->td = td - d1->td;
728: for ( i = 0; i < n; i++ )
729: d->d[i] = w[i] - d1->d[i];
730: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
731: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
732: NEWDL(d,n); d->td = td - d2->td;
733: for ( i = 0; i < n; i++ )
734: d->d[i] = w[i] - d2->d[i];
735: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
736: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
737: addmd_destructive(mod,t,u,rp);
738: }
739:
740: /*
741: * m-reduction
1.13 noro 742: * do content reduction over Z or Q(x,...)
743: * do nothing over finite fields
1.7 noro 744: *
745: */
746:
1.20 noro 747: void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp,DP *multp)
1.7 noro 748: {
749: int i,n;
750: DL d1,d2,d;
751: MP m;
752: DP t,s,r,h;
753: Q c,c1,c2;
754: N gn,tn;
755: P g,a;
1.23 noro 756: P p[2];
1.7 noro 757:
758: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
759: NEWDL(d,n); d->td = d1->td - d2->td;
760: for ( i = 0; i < n; i++ )
761: d->d[i] = d1->d[i]-d2->d[i];
762: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
1.23 noro 763: if ( dp_fcoeffs == N_GFS ) {
764: p[0] = (P)c1; p[1] = (P)c2;
765: gcdsf(CO,p,2,&g);
766: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
767: } else if ( dp_fcoeffs ) {
1.7 noro 768: /* do nothing */
769: } else if ( INT(c1) && INT(c2) ) {
770: gcdn(NM(c1),NM(c2),&gn);
771: if ( !UNIN(gn) ) {
772: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
773: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
774: }
775: } else {
776: ezgcdpz(CO,(P)c1,(P)c2,&g);
777: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
778: }
779: NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
780: *multp = s;
781: muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
782: muldc(CO,p0,(P)c2,&h);
783: *head = h; *rest = r; *dnp = (P)c2;
784: }
785:
1.13 noro 786: /* m-reduction over a field */
787:
1.20 noro 788: void dp_red_f(DP p1,DP p2,DP *rest)
1.13 noro 789: {
790: int i,n;
791: DL d1,d2,d;
792: MP m;
1.20 noro 793: DP t,s;
1.13 noro 794: Obj a,b;
795:
796: n = p1->nv;
797: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
798:
799: NEWDL(d,n); d->td = d1->td - d2->td;
800: for ( i = 0; i < n; i++ )
801: d->d[i] = d1->d[i]-d2->d[i];
802:
803: NEWMP(m); m->dl = d;
804: divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
805: C(m) = (P)b;
806: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
807:
808: muld(CO,s,p2,&t); addd(CO,p1,t,rest);
809: }
810:
1.20 noro 811: void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp)
1.7 noro 812: {
813: int i,n;
814: DL d1,d2,d;
815: MP m;
816: DP t,s,r,h;
817: P c1,c2,g,u;
818:
819: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
820: NEWDL(d,n); d->td = d1->td - d2->td;
821: for ( i = 0; i < n; i++ )
822: d->d[i] = d1->d[i]-d2->d[i];
823: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
824: gcdprsmp(CO,mod,c1,c2,&g);
825: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
826: if ( NUM(c2) ) {
827: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
828: }
829: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1.11 noro 830: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7 noro 831: if ( NUM(c2) ) {
832: addmd(CO,mod,p1,t,&r); h = p0;
833: } else {
834: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
835: }
836: *head = h; *rest = r; *dnp = c2;
837: }
838:
1.10 noro 839: struct oEGT eg_red_mod;
840:
1.20 noro 841: void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp)
1.7 noro 842: {
843: int i,n;
844: DL d1,d2,d;
845: MP m;
846: DP t,s;
1.16 noro 847: int c,c1,c2;
848: extern int do_weyl;
1.7 noro 849:
850: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
851: _NEWDL(d,n); d->td = d1->td - d2->td;
852: for ( i = 0; i < n; i++ )
853: d->d[i] = d1->d[i]-d2->d[i];
1.16 noro 854: c = invm(ITOS(BDY(p2)->c),mod);
855: c2 = ITOS(BDY(p1)->c);
856: DMAR(c,c2,0,mod,c1);
1.7 noro 857: _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
1.16 noro 858: #if 0
1.7 noro 859: _MKDP(n,m,s); s->sugar = d->td;
860: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 861: #else
862: if ( do_weyl ) {
1.19 noro 863: _MKDP(n,m,s); s->sugar = d->td;
864: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 865: } else {
866: _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
867: }
868: #endif
1.10 noro 869: /* get_eg(&t0); */
1.7 noro 870: _addmd_destructive(mod,p1,t,rp);
1.10 noro 871: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 872: }
873:
874: /*
875: * normal form computation
876: *
877: */
1.5 noro 878:
1.20 noro 879: void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp)
1.5 noro 880: {
881: DP u,p,d,s,t,dmy;
882: NODE l;
883: MP m,mr;
884: int i,n;
885: int *wb;
886: int sugar,psugar;
887: P dn,tdn,tdn1;
888:
889: dn = (P)ONE;
890: if ( !g ) {
891: *rp = 0; *dnp = dn; return;
892: }
893: for ( n = 0, l = b; l; l = NEXT(l), n++ );
894: wb = (int *)ALLOCA(n*sizeof(int));
895: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
896: wb[i] = QTOS((Q)BDY(l));
897: sugar = g->sugar;
898: for ( d = 0; g; ) {
899: for ( u = 0, i = 0; i < n; i++ ) {
900: if ( dp_redble(g,p = ps[wb[i]]) ) {
901: dp_red(d,g,p,&t,&u,&tdn,&dmy);
902: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
903: sugar = MAX(sugar,psugar);
904: if ( !u ) {
905: if ( d )
906: d->sugar = sugar;
907: *rp = d; *dnp = dn; return;
908: } else {
909: d = t;
910: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
911: }
912: break;
913: }
914: }
915: if ( u )
916: g = u;
917: else if ( !full ) {
918: if ( g ) {
919: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
920: }
921: *rp = g; *dnp = dn; return;
922: } else {
923: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
924: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
925: addd(CO,d,t,&s); d = s;
926: dp_rest(g,&t); g = t;
927: }
928: }
929: if ( d )
930: d->sugar = sugar;
931: *rp = d; *dnp = dn;
932: }
933:
1.13 noro 934: /* nf computation over Z */
935:
1.20 noro 936: void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
1.5 noro 937: {
938: DP u,p,d,s,t,dmy1;
939: P dmy;
940: NODE l;
941: MP m,mr;
942: int i,n;
943: int *wb;
944: int hmag;
945: int sugar,psugar;
946:
947: if ( !g ) {
948: *rp = 0; return;
949: }
950: for ( n = 0, l = b; l; l = NEXT(l), n++ );
951: wb = (int *)ALLOCA(n*sizeof(int));
952: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
953: wb[i] = QTOS((Q)BDY(l));
1.12 noro 954:
1.13 noro 955: hmag = multiple*HMAG(g);
1.5 noro 956: sugar = g->sugar;
1.12 noro 957:
1.5 noro 958: for ( d = 0; g; ) {
959: for ( u = 0, i = 0; i < n; i++ ) {
960: if ( dp_redble(g,p = ps[wb[i]]) ) {
961: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
962: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
963: sugar = MAX(sugar,psugar);
964: if ( !u ) {
965: if ( d )
966: d->sugar = sugar;
967: *rp = d; return;
968: }
969: d = t;
970: break;
971: }
972: }
973: if ( u ) {
974: g = u;
975: if ( d ) {
1.13 noro 976: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 977: dp_ptozp2(d,g,&t,&u); d = t; g = u;
978: hmag = multiple*HMAG(d);
979: }
980: } else {
1.13 noro 981: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 982: dp_ptozp(g,&t); g = t;
983: hmag = multiple*HMAG(g);
984: }
985: }
986: }
987: else if ( !full ) {
988: if ( g ) {
989: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
990: }
991: *rp = g; return;
992: } else {
993: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
994: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
995: addd(CO,d,t,&s); d = s;
996: dp_rest(g,&t); g = t;
997:
998: }
999: }
1000: if ( d )
1001: d->sugar = sugar;
1002: *rp = d;
1003: }
1004:
1.13 noro 1005: /* nf computation over a field */
1006:
1.20 noro 1007: void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp)
1.13 noro 1008: {
1009: DP u,p,d,s,t;
1010: NODE l;
1011: MP m,mr;
1012: int i,n;
1013: int *wb;
1014: int sugar,psugar;
1015:
1016: if ( !g ) {
1017: *rp = 0; return;
1018: }
1019: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1020: wb = (int *)ALLOCA(n*sizeof(int));
1021: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1022: wb[i] = QTOS((Q)BDY(l));
1023:
1024: sugar = g->sugar;
1025: for ( d = 0; g; ) {
1026: for ( u = 0, i = 0; i < n; i++ ) {
1027: if ( dp_redble(g,p = ps[wb[i]]) ) {
1028: dp_red_f(g,p,&u);
1029: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1030: sugar = MAX(sugar,psugar);
1031: if ( !u ) {
1032: if ( d )
1033: d->sugar = sugar;
1034: *rp = d; return;
1035: }
1036: break;
1037: }
1038: }
1039: if ( u )
1040: g = u;
1041: else if ( !full ) {
1042: if ( g ) {
1043: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1044: }
1045: *rp = g; return;
1046: } else {
1047: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1048: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1049: addd(CO,d,t,&s); d = s;
1050: dp_rest(g,&t); g = t;
1051: }
1052: }
1053: if ( d )
1054: d->sugar = sugar;
1055: *rp = d;
1056: }
1057:
1058: /* nf computation over GF(mod) (only for internal use) */
1059:
1.20 noro 1060: void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1061: {
1062: DP u,p,d,s,t;
1063: P dmy;
1064: NODE l;
1065: MP m,mr;
1066: int sugar,psugar;
1067:
1068: if ( !g ) {
1069: *rp = 0; return;
1070: }
1071: sugar = g->sugar;
1072: for ( d = 0; g; ) {
1073: for ( u = 0, l = b; l; l = NEXT(l) ) {
1074: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1075: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
1076: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1077: sugar = MAX(sugar,psugar);
1078: if ( !u ) {
1079: if ( d )
1080: d->sugar = sugar;
1081: *rp = d; return;
1082: }
1083: d = t;
1084: break;
1085: }
1086: }
1087: if ( u )
1088: g = u;
1089: else if ( !full ) {
1090: if ( g ) {
1091: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1092: }
1093: *rp = g; return;
1094: } else {
1095: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1096: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1097: addmd(CO,mod,d,t,&s); d = s;
1098: dp_rest(g,&t); g = t;
1099: }
1100: }
1101: if ( d )
1102: d->sugar = sugar;
1103: *rp = d;
1104: }
1105:
1.20 noro 1106: void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp)
1.5 noro 1107: {
1108: DP u,p,d,s,t;
1109: NODE l;
1110: MP m,mr;
1111: int i,n;
1112: int *wb;
1113: int sugar,psugar;
1114: P dn,tdn,tdn1;
1115:
1116: dn = (P)ONEM;
1117: if ( !g ) {
1118: *rp = 0; *dnp = dn; return;
1119: }
1120: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1121: wb = (int *)ALLOCA(n*sizeof(int));
1122: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1123: wb[i] = QTOS((Q)BDY(l));
1124: sugar = g->sugar;
1125: for ( d = 0; g; ) {
1126: for ( u = 0, i = 0; i < n; i++ ) {
1127: if ( dp_redble(g,p = ps[wb[i]]) ) {
1128: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
1129: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1130: sugar = MAX(sugar,psugar);
1131: if ( !u ) {
1132: if ( d )
1133: d->sugar = sugar;
1134: *rp = d; *dnp = dn; return;
1135: } else {
1136: d = t;
1137: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1138: }
1139: break;
1140: }
1141: }
1142: if ( u )
1143: g = u;
1144: else if ( !full ) {
1145: if ( g ) {
1146: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1147: }
1148: *rp = g; *dnp = dn; return;
1149: } else {
1150: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1151: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1152: addmd(CO,mod,d,t,&s); d = s;
1153: dp_rest(g,&t); g = t;
1154: }
1155: }
1156: if ( d )
1157: d->sugar = sugar;
1158: *rp = d; *dnp = dn;
1159: }
1160:
1.20 noro 1161: void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1162: {
1.20 noro 1163: DP u,p,d;
1.7 noro 1164: NODE l;
1.20 noro 1165: MP m,mrd;
1166: int sugar,psugar,n,h_reducible;
1.5 noro 1167:
1.7 noro 1168: if ( !g ) {
1169: *rp = 0; return;
1.5 noro 1170: }
1.7 noro 1171: sugar = g->sugar;
1172: n = g->nv;
1173: for ( d = 0; g; ) {
1174: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
1175: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1176: h_reducible = 1;
1177: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1178: _dp_red_mod_destructive(g,p,mod,&u); g = u;
1179: sugar = MAX(sugar,psugar);
1180: if ( !g ) {
1181: if ( d )
1182: d->sugar = sugar;
1183: _dptodp(d,rp); _free_dp(d); return;
1184: }
1185: break;
1186: }
1187: }
1188: if ( !h_reducible ) {
1189: /* head term is not reducible */
1190: if ( !full ) {
1191: if ( g )
1192: g->sugar = sugar;
1193: _dptodp(g,rp); _free_dp(g); return;
1194: } else {
1195: m = BDY(g);
1196: if ( NEXT(m) ) {
1197: BDY(g) = NEXT(m); NEXT(m) = 0;
1198: } else {
1199: _FREEDP(g); g = 0;
1200: }
1201: if ( d ) {
1202: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1203: NEXT(mrd) = m;
1204: } else {
1205: _MKDP(n,m,d);
1206: }
1207: }
1208: }
1.5 noro 1209: }
1.7 noro 1210: if ( d )
1211: d->sugar = sugar;
1212: _dptodp(d,rp); _free_dp(d);
1.5 noro 1213: }
1.13 noro 1214:
1215: /* reduction by linear base over a field */
1216:
1.20 noro 1217: void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p)
1.13 noro 1218: {
1219: DP r1,r2,b1,b2,t,s;
1220: Obj c,c1,c2;
1221: NODE l,b;
1222: int n;
1223:
1224: if ( !p1 ) {
1225: *r1p = p1; *r2p = p2; return;
1226: }
1227: n = p1->nv;
1228: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1229: if ( !r1 ) {
1230: *r1p = r1; *r2p = r2; return;
1231: }
1232: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1233: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1234: b2 = (DP)BDY(NEXT(b));
1235: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
1236: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
1237: muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
1238: muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
1239: }
1240: }
1241: *r1p = r1; *r2p = r2;
1242: }
1243:
1244: /* reduction by linear base over GF(mod) */
1.5 noro 1245:
1.20 noro 1246: void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p)
1.5 noro 1247: {
1.7 noro 1248: DP r1,r2,b1,b2,t,s;
1249: P c;
1250: MQ c1,c2;
1251: NODE l,b;
1252: int n;
1253:
1254: if ( !p1 ) {
1255: *r1p = p1; *r2p = p2; return;
1256: }
1257: n = p1->nv;
1258: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1259: if ( !r1 ) {
1260: *r1p = r1; *r2p = r2; return;
1261: }
1262: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1263: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1264: b2 = (DP)BDY(NEXT(b));
1265: invmq(mod,(MQ)BDY(b1)->c,&c1);
1266: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
1267: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
1268: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
1269: }
1270: }
1271: *r1p = r1; *r2p = r2;
1.5 noro 1272: }
1273:
1.20 noro 1274: void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp)
1.5 noro 1275: {
1.7 noro 1276: DP s,t,u;
1277: MP m;
1278: DL h;
1279: int i,n;
1280:
1281: if ( !p ) {
1282: *rp = p; return;
1283: }
1284: n = p->nv;
1285: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1286: h = m->dl;
1287: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1288: i++;
1289: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1290: addmd(CO,mod,s,t,&u); s = u;
1.24 noro 1291: }
1292: *rp = s;
1293: }
1294:
1295: void dp_nf_tab_f(DP p,LIST *tab,DP *rp)
1296: {
1297: DP s,t,u;
1298: MP m;
1299: DL h;
1300: int i,n;
1301:
1302: if ( !p ) {
1303: *rp = p; return;
1304: }
1305: n = p->nv;
1306: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1307: h = m->dl;
1308: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1309: i++;
1310: muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1311: addd(CO,s,t,&u); s = u;
1.7 noro 1312: }
1313: *rp = s;
1.5 noro 1314: }
1315:
1.7 noro 1316: /*
1317: * setting flags
1.30 noro 1318: * call create_order_spec with vl=0 to set old type order.
1.7 noro 1319: *
1320: */
1321:
1.27 noro 1322: int create_order_spec(VL vl,Obj obj,struct order_spec **specp)
1.5 noro 1323: {
1.37 ! noro 1324: int i,j,n,s,row,col,ret;
1.27 noro 1325: struct order_spec *spec;
1.7 noro 1326: struct order_pair *l;
1327: NODE node,t,tn;
1328: MAT m;
1329: pointer **b;
1330: int **w;
1.5 noro 1331:
1.37 ! noro 1332: if ( vl && obj && OID(obj) == O_LIST ) {
! 1333: ret = create_composite_order_spec(vl,(LIST)obj,specp);
! 1334: if ( show_orderspec )
! 1335: print_composite_order_spec(*specp);
! 1336: return ret;
! 1337: }
1.27 noro 1338:
1339: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 1340: if ( !obj || NUM(obj) ) {
1341: spec->id = 0; spec->obj = obj;
1342: spec->ord.simple = QTOS((Q)obj);
1343: return 1;
1344: } else if ( OID(obj) == O_LIST ) {
1345: node = BDY((LIST)obj);
1346: for ( n = 0, t = node; t; t = NEXT(t), n++ );
1347: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1348: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
1349: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
1350: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
1351: s += l[i].length;
1352: }
1353: spec->id = 1; spec->obj = obj;
1354: spec->ord.block.order_pair = l;
1355: spec->ord.block.length = n; spec->nv = s;
1356: return 1;
1357: } else if ( OID(obj) == O_MAT ) {
1358: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
1359: w = almat(row,col);
1360: for ( i = 0; i < row; i++ )
1361: for ( j = 0; j < col; j++ )
1362: w[i][j] = QTOS((Q)b[i][j]);
1363: spec->id = 2; spec->obj = obj;
1364: spec->nv = col; spec->ord.matrix.row = row;
1365: spec->ord.matrix.matrix = w;
1366: return 1;
1367: } else
1.5 noro 1368: return 0;
1369: }
1370:
1.28 noro 1371: void print_composite_order_spec(struct order_spec *spec)
1372: {
1373: int nv,n,len,i,j,k,start;
1374: struct weight_or_block *worb;
1375:
1376: nv = spec->nv;
1377: n = spec->ord.composite.length;
1378: worb = spec->ord.composite.w_or_b;
1379: for ( i = 0; i < n; i++, worb++ ) {
1380: len = worb->length;
1381: printf("[ ");
1382: switch ( worb->type ) {
1383: case IS_DENSE_WEIGHT:
1384: for ( j = 0; j < len; j++ )
1385: printf("%d ",worb->body.dense_weight[j]);
1386: for ( ; j < nv; j++ )
1387: printf("0 ");
1388: break;
1389: case IS_SPARSE_WEIGHT:
1390: for ( j = 0, k = 0; j < nv; j++ )
1391: if ( j == worb->body.sparse_weight[k].pos )
1392: printf("%d ",worb->body.sparse_weight[k++].value);
1393: else
1394: printf("0 ");
1395: break;
1396: case IS_BLOCK:
1397: start = worb->body.block.start;
1398: for ( j = 0; j < start; j++ ) printf("0 ");
1399: switch ( worb->body.block.order ) {
1400: case 0:
1401: for ( k = 0; k < len; k++, j++ ) printf("R ");
1402: break;
1403: case 1:
1404: for ( k = 0; k < len; k++, j++ ) printf("G ");
1405: break;
1406: case 2:
1407: for ( k = 0; k < len; k++, j++ ) printf("L ");
1408: break;
1409: }
1410: for ( ; j < nv; j++ ) printf("0 ");
1411: break;
1412: }
1413: printf("]\n");
1414: }
1415: }
1416:
1.37 ! noro 1417: int comp_sw(struct sparse_weight *a, struct sparse_weight *b)
! 1418: {
! 1419: if ( a->pos > b->pos ) return 1;
! 1420: else if ( a->pos < b->pos ) return -1;
! 1421: else return 0;
! 1422: }
! 1423:
1.27 noro 1424: /* order = [w_or_b, w_or_b, ... ] */
1425: /* w_or_b = w or b */
1426: /* w = [1,2,...] or [x,1,y,2,...] */
1427: /* b = [@lex,x,y,...,z] etc */
1428:
1429: int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp)
1430: {
1431: NODE wb,t,p;
1432: struct order_spec *spec;
1433: VL tvl;
1.29 noro 1434: int n,i,j,k,l,start,end,len,w;
1.27 noro 1435: int *dw;
1436: struct sparse_weight *sw;
1437: struct weight_or_block *w_or_b;
1438: Obj a0;
1439: NODE a;
1.29 noro 1440: V v,sv,ev;
1441: SYMBOL sym;
1442: int *top;
1.27 noro 1443:
1444: /* l = number of vars in vl */
1445: for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ );
1446: /* n = number of primitives in order */
1447: wb = BDY(order);
1448: n = length(wb);
1449: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1450: spec->id = 3;
1451: spec->obj = (Obj)order;
1452: spec->nv = l;
1453: spec->ord.composite.length = n;
1.28 noro 1454: w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *)
1.29 noro 1455: MALLOC(sizeof(struct weight_or_block)*(n+1));
1456:
1457: /* top : register the top variable in each w_or_b specification */
1458: top = (int *)ALLOCA(l*sizeof(int));
1459: for ( i = 0; i < l; i++ ) top[i] = 0;
1460:
1.28 noro 1461: for ( t = wb, i = 0; t; t = NEXT(t), i++ ) {
1.30 noro 1462: if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST )
1463: error("a list of lists must be specified for the key \"order\"");
1.28 noro 1464: a = BDY((LIST)BDY(t));
1.27 noro 1465: len = length(a);
1466: a0 = (Obj)BDY(a);
1467: if ( !a0 || OID(a0) == O_N ) {
1.28 noro 1468: /* a is a dense weight vector */
1.27 noro 1469: dw = (int *)MALLOC(sizeof(int)*len);
1.30 noro 1470: for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) {
1471: if ( !INT((Q)BDY(p)) )
1472: error("a dense weight vector must be specified as a list of integers");
1.27 noro 1473: dw[j] = QTOS((Q)BDY(p));
1.30 noro 1474: }
1.27 noro 1475: w_or_b[i].type = IS_DENSE_WEIGHT;
1476: w_or_b[i].length = len;
1477: w_or_b[i].body.dense_weight = dw;
1.29 noro 1478:
1479: /* find the top */
1480: for ( k = 0; k < len && !dw[k]; k++ );
1481: if ( k < len ) top[k] = 1;
1482:
1.27 noro 1483: } else if ( OID(a0) == O_P ) {
1.28 noro 1484: /* a is a sparse weight vector */
1485: len >>= 1;
1.27 noro 1486: sw = (struct sparse_weight *)
1487: MALLOC(sizeof(struct sparse_weight)*len);
1488: for ( j = 0, p = a; j < len; j++ ) {
1.30 noro 1489: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
1490: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 1491: v = VR((P)BDY(p)); p = NEXT(p);
1.27 noro 1492: for ( tvl = vl, k = 0; tvl && tvl->v != v;
1493: k++, tvl = NEXT(tvl) );
1494: if ( !tvl )
1.30 noro 1495: error("invalid variable name in a sparse weight vector");
1.27 noro 1496: sw[j].pos = k;
1.30 noro 1497: if ( !INT((Q)BDY(p)) )
1498: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 1499: sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p);
1.27 noro 1500: }
1.37 ! noro 1501: qsort(sw,len,sizeof(struct sparse_weight),
! 1502: (int (*)(const void *,const void *))comp_sw);
1.27 noro 1503: w_or_b[i].type = IS_SPARSE_WEIGHT;
1504: w_or_b[i].length = len;
1505: w_or_b[i].body.sparse_weight = sw;
1.29 noro 1506:
1507: /* find the top */
1508: for ( k = 0; k < len && !sw[k].value; k++ );
1509: if ( k < len ) top[sw[k].pos] = 1;
1510: } else if ( OID(a0) == O_RANGE ) {
1511: /* [range(v1,v2),w] */
1512: sv = VR((P)(((RANGE)a0)->start));
1513: ev = VR((P)(((RANGE)a0)->end));
1514: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
1515: if ( !tvl )
1516: error("invalid range");
1517: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
1518: if ( !tvl )
1519: error("invalid range");
1520: len = end-start+1;
1521: sw = (struct sparse_weight *)
1522: MALLOC(sizeof(struct sparse_weight)*len);
1523: w = QTOS((Q)BDY(NEXT(a)));
1524: for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) );
1525: for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) {
1526: sw[j].pos = k;
1527: sw[j].value = w;
1528: }
1529: w_or_b[i].type = IS_SPARSE_WEIGHT;
1530: w_or_b[i].length = len;
1531: w_or_b[i].body.sparse_weight = sw;
1532:
1533: /* register the top */
1534: if ( w ) top[start] = 1;
1.28 noro 1535: } else if ( OID(a0) == O_SYMBOL ) {
1536: /* a is a block */
1.29 noro 1537: sym = (SYMBOL)a0; a = NEXT(a); len--;
1538: if ( OID((Obj)BDY(a)) == O_RANGE ) {
1539: sv = VR((P)(((RANGE)BDY(a))->start));
1540: ev = VR((P)(((RANGE)BDY(a))->end));
1541: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
1542: if ( !tvl )
1543: error("invalid range");
1544: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
1545: if ( !tvl )
1546: error("invalid range");
1547: len = end-start+1;
1548: } else {
1549: for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a));
1.28 noro 1550: tvl = NEXT(tvl), start++ );
1.29 noro 1551: for ( p = NEXT(a), tvl = NEXT(tvl); p;
1.30 noro 1552: p = NEXT(p), tvl = NEXT(tvl) ) {
1553: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
1554: error("a block must be specified as [ordsymbol,var1,var2,...]");
1.29 noro 1555: if ( tvl->v != VR((P)BDY(p)) ) break;
1.30 noro 1556: }
1.29 noro 1557: if ( p )
1.30 noro 1558: error("a block must be contiguous in the variable list");
1.29 noro 1559: }
1.28 noro 1560: w_or_b[i].type = IS_BLOCK;
1561: w_or_b[i].length = len;
1562: w_or_b[i].body.block.start = start;
1563: if ( !strcmp(sym->name,"@grlex") )
1564: w_or_b[i].body.block.order = 0;
1565: else if ( !strcmp(sym->name,"@glex") )
1566: w_or_b[i].body.block.order = 1;
1567: else if ( !strcmp(sym->name,"@lex") )
1568: w_or_b[i].body.block.order = 2;
1569: else
1.29 noro 1570: error("invalid ordername");
1571: /* register the tops */
1572: for ( j = 0, k = start; j < len; j++, k++ )
1573: top[k] = 1;
1.28 noro 1574: }
1.29 noro 1575: }
1576: for ( k = 0; k < l && top[k]; k++ );
1577: if ( k < l ) {
1578: /* incomplete order specification; add @grlex */
1579: w_or_b[n].type = IS_BLOCK;
1580: w_or_b[n].length = l;
1581: w_or_b[n].body.block.start = 0;
1582: w_or_b[n].body.block.order = 0;
1583: spec->ord.composite.length = n+1;
1.27 noro 1584: }
1585: }
1586:
1.35 noro 1587: /* module order spec */
1588:
1589: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s)
1590: {
1591: struct modorder_spec *spec;
1592: NODE n,t;
1593: LIST list;
1594: int *ds;
1595: int i,l;
1596: Q q;
1597:
1598: *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec));
1599: spec->id = id;
1600: if ( shift ) {
1601: n = BDY(shift);
1602: spec->len = l = length(n);
1603: spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int));
1604: for ( t = n, i = 0; t; t = NEXT(t), i++ )
1605: ds[i] = QTOS((Q)BDY(t));
1606: } else {
1607: spec->len = 0;
1608: spec->degree_shift = 0;
1609: }
1610: STOQ(id,q);
1611: n = mknode(2,q,shift);
1612: MKLIST(list,n);
1613: spec->obj = (Obj)list;
1614: }
1615:
1.7 noro 1616: /*
1617: * converters
1618: *
1619: */
1620:
1.20 noro 1621: void dp_homo(DP p,DP *rp)
1.5 noro 1622: {
1.7 noro 1623: MP m,mr,mr0;
1624: int i,n,nv,td;
1625: DL dl,dlh;
1.5 noro 1626:
1.7 noro 1627: if ( !p )
1628: *rp = 0;
1629: else {
1630: n = p->nv; nv = n + 1;
1631: m = BDY(p); td = sugard(m);
1632: for ( mr0 = 0; m; m = NEXT(m) ) {
1633: NEXTMP(mr0,mr); mr->c = m->c;
1634: dl = m->dl;
1635: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1636: dlh->td = td;
1637: for ( i = 0; i < n; i++ )
1638: dlh->d[i] = dl->d[i];
1639: dlh->d[n] = td - dl->td;
1640: }
1641: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 1642: }
1643: }
1644:
1.20 noro 1645: void dp_dehomo(DP p,DP *rp)
1.5 noro 1646: {
1.7 noro 1647: MP m,mr,mr0;
1648: int i,n,nv;
1649: DL dl,dlh;
1.5 noro 1650:
1.7 noro 1651: if ( !p )
1652: *rp = 0;
1653: else {
1654: n = p->nv; nv = n - 1;
1655: m = BDY(p);
1656: for ( mr0 = 0; m; m = NEXT(m) ) {
1657: NEXTMP(mr0,mr); mr->c = m->c;
1658: dlh = m->dl;
1659: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1660: dl->td = dlh->td - dlh->d[nv];
1661: for ( i = 0; i < nv; i++ )
1662: dl->d[i] = dlh->d[i];
1663: }
1664: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1665: }
1.5 noro 1666: }
1667:
1.20 noro 1668: void dp_mod(DP p,int mod,NODE subst,DP *rp)
1.5 noro 1669: {
1.7 noro 1670: MP m,mr,mr0;
1671: P t,s,s1;
1672: V v;
1673: NODE tn;
1.5 noro 1674:
1.7 noro 1675: if ( !p )
1676: *rp = 0;
1677: else {
1678: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1679: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
1680: v = VR((P)BDY(tn)); tn = NEXT(tn);
1681: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
1682: }
1683: ptomp(mod,s,&t);
1684: if ( t ) {
1685: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
1686: }
1687: }
1688: if ( mr0 ) {
1689: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1690: } else
1691: *rp = 0;
1692: }
1.5 noro 1693: }
1694:
1.20 noro 1695: void dp_rat(DP p,DP *rp)
1.5 noro 1696: {
1.7 noro 1697: MP m,mr,mr0;
1.5 noro 1698:
1.7 noro 1699: if ( !p )
1700: *rp = 0;
1701: else {
1702: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1703: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5 noro 1704: }
1.7 noro 1705: if ( mr0 ) {
1706: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1707: } else
1708: *rp = 0;
1.5 noro 1709: }
1710: }
1711:
1712:
1.27 noro 1713: void homogenize_order(struct order_spec *old,int n,struct order_spec **newp)
1.5 noro 1714: {
1.7 noro 1715: struct order_pair *l;
1716: int length,nv,row,i,j;
1717: int **newm,**oldm;
1.27 noro 1718: struct order_spec *new;
1.31 noro 1719: int onv,nnv,nlen,olen,owlen;
1720: struct weight_or_block *owb,*nwb;
1.5 noro 1721:
1.27 noro 1722: *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 1723: switch ( old->id ) {
1724: case 0:
1725: switch ( old->ord.simple ) {
1726: case 0:
1727: new->id = 0; new->ord.simple = 0; break;
1728: case 1:
1729: l = (struct order_pair *)
1730: MALLOC_ATOMIC(2*sizeof(struct order_pair));
1731: l[0].length = n; l[0].order = 1;
1732: l[1].length = 1; l[1].order = 2;
1733: new->id = 1;
1734: new->ord.block.order_pair = l;
1735: new->ord.block.length = 2; new->nv = n+1;
1736: break;
1737: case 2:
1738: new->id = 0; new->ord.simple = 1; break;
1739: case 3: case 4: case 5:
1740: new->id = 0; new->ord.simple = old->ord.simple+3;
1741: dp_nelim = n-1; break;
1742: case 6: case 7: case 8: case 9:
1743: new->id = 0; new->ord.simple = old->ord.simple; break;
1744: default:
1745: error("homogenize_order : invalid input");
1746: }
1747: break;
1748: case 1:
1749: length = old->ord.block.length;
1750: l = (struct order_pair *)
1751: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
1752: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
1753: l[length].order = 2; l[length].length = 1;
1754: new->id = 1; new->nv = n+1;
1755: new->ord.block.order_pair = l;
1756: new->ord.block.length = length+1;
1757: break;
1758: case 2:
1759: nv = old->nv; row = old->ord.matrix.row;
1760: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
1761: for ( i = 0; i <= nv; i++ )
1762: newm[0][i] = 1;
1763: for ( i = 0; i < row; i++ ) {
1764: for ( j = 0; j < nv; j++ )
1765: newm[i+1][j] = oldm[i][j];
1766: newm[i+1][j] = 0;
1767: }
1768: new->id = 2; new->nv = nv+1;
1769: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1.31 noro 1770: break;
1771: case 3:
1772: onv = old->nv;
1773: nnv = onv+1;
1774: olen = old->ord.composite.length;
1775: nlen = olen+1;
1776: owb = old->ord.composite.w_or_b;
1777: nwb = (struct weight_or_block *)
1778: MALLOC(nlen*sizeof(struct weight_or_block));
1779: for ( i = 0; i < olen; i++ ) {
1780: nwb[i].type = owb[i].type;
1781: switch ( owb[i].type ) {
1782: case IS_DENSE_WEIGHT:
1783: owlen = owb[i].length;
1784: nwb[i].length = owlen+1;
1785: nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int));
1786: for ( j = 0; j < owlen; j++ )
1787: nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j];
1788: nwb[i].body.dense_weight[owlen] = 0;
1789: break;
1790: case IS_SPARSE_WEIGHT:
1791: nwb[i].length = owb[i].length;
1792: nwb[i].body.sparse_weight = owb[i].body.sparse_weight;
1793: break;
1794: case IS_BLOCK:
1795: nwb[i].length = owb[i].length;
1796: nwb[i].body.block = owb[i].body.block;
1797: break;
1798: }
1799: }
1800: nwb[i].type = IS_SPARSE_WEIGHT;
1801: nwb[i].body.sparse_weight =
1802: (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
1803: nwb[i].body.sparse_weight[0].pos = onv;
1804: nwb[i].body.sparse_weight[0].value = 1;
1805: new->id = 3;
1806: new->nv = nnv;
1807: new->ord.composite.length = nlen;
1808: new->ord.composite.w_or_b = nwb;
1809: print_composite_order_spec(new);
1.7 noro 1810: break;
1811: default:
1812: error("homogenize_order : invalid input");
1.5 noro 1813: }
1.7 noro 1814: }
1815:
1.20 noro 1816: void qltozl(Q *w,int n,Q *dvr)
1.7 noro 1817: {
1818: N nm,dn;
1819: N g,l1,l2,l3;
1820: Q c,d;
1821: int i;
1822: struct oVECT v;
1.5 noro 1823:
1824: for ( i = 0; i < n; i++ )
1.7 noro 1825: if ( w[i] && !INT(w[i]) )
1826: break;
1827: if ( i == n ) {
1828: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
1829: igcdv(&v,dvr); return;
1830: }
1831: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
1832: for ( i = 1; i < n; i++ ) {
1833: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
1834: gcdn(nm,NM(c),&g); nm = g;
1835: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 1836: }
1.7 noro 1837: if ( UNIN(dn) )
1838: NTOQ(nm,1,d);
1839: else
1840: NDTOQ(nm,dn,1,d);
1841: *dvr = d;
1842: }
1.5 noro 1843:
1.20 noro 1844: int comp_nm(Q *a,Q *b)
1.7 noro 1845: {
1846: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
1847: }
1848:
1.20 noro 1849: void sortbynm(Q *w,int n)
1.7 noro 1850: {
1851: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
1852: }
1.5 noro 1853:
1854:
1.7 noro 1855: /*
1856: * simple operations
1857: *
1858: */
1.5 noro 1859:
1.20 noro 1860: int dp_redble(DP p1,DP p2)
1.7 noro 1861: {
1862: int i,n;
1863: DL d1,d2;
1.5 noro 1864:
1.7 noro 1865: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1866: if ( d1->td < d2->td )
1867: return 0;
1868: else {
1869: for ( i = 0, n = p1->nv; i < n; i++ )
1870: if ( d1->d[i] < d2->d[i] )
1871: return 0;
1872: return 1;
1.5 noro 1873: }
1874: }
1875:
1.20 noro 1876: void dp_subd(DP p1,DP p2,DP *rp)
1.5 noro 1877: {
1.7 noro 1878: int i,n;
1.5 noro 1879: DL d1,d2,d;
1880: MP m;
1.7 noro 1881: DP s;
1.5 noro 1882:
1883: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 noro 1884: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 1885: for ( i = 0; i < n; i++ )
1.7 noro 1886: d->d[i] = d1->d[i]-d2->d[i];
1887: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1888: MKDP(n,m,s); s->sugar = d->td;
1889: *rp = s;
1890: }
1891:
1.20 noro 1892: void dltod(DL d,int n,DP *rp)
1.7 noro 1893: {
1894: MP m;
1895: DP s;
1896:
1897: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1898: MKDP(n,m,s); s->sugar = d->td;
1899: *rp = s;
1.5 noro 1900: }
1901:
1.20 noro 1902: void dp_hm(DP p,DP *rp)
1.5 noro 1903: {
1904: MP m,mr;
1905:
1906: if ( !p )
1907: *rp = 0;
1908: else {
1909: m = BDY(p);
1910: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
1911: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1912: }
1913: }
1914:
1.35 noro 1915: void dp_ht(DP p,DP *rp)
1916: {
1917: MP m,mr;
1918:
1919: if ( !p )
1920: *rp = 0;
1921: else {
1922: m = BDY(p);
1923: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
1924: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1925: }
1926: }
1927:
1.20 noro 1928: void dp_rest(DP p,DP *rp)
1.5 noro 1929: {
1930: MP m;
1931:
1932: m = BDY(p);
1933: if ( !NEXT(m) )
1934: *rp = 0;
1935: else {
1936: MKDP(p->nv,NEXT(m),*rp);
1937: if ( *rp )
1938: (*rp)->sugar = p->sugar;
1939: }
1940: }
1941:
1.20 noro 1942: DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl)
1.5 noro 1943: {
1.21 noro 1944: register int i, *d1, *d2, *d, td;
1.5 noro 1945:
1946: if ( !dl ) NEWDL(dl,nv);
1947: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1.21 noro 1948: for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) {
1949: *d = *d1 > *d2 ? *d1 : *d2;
1950: td += MUL_WEIGHT(*d,i);
1951: }
1.5 noro 1952: dl->td = td;
1953: return dl;
1954: }
1955:
1.20 noro 1956: int dl_equal(int nv,DL dl1,DL dl2)
1.5 noro 1957: {
1958: register int *d1, *d2, n;
1959:
1960: if ( dl1->td != dl2->td ) return 0;
1961: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
1962: if ( *d1 != *d2 ) return 0;
1963: return 1;
1964: }
1965:
1.20 noro 1966: int dp_nt(DP p)
1.5 noro 1967: {
1968: int i;
1969: MP m;
1970:
1971: if ( !p )
1972: return 0;
1973: else {
1974: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
1975: return i;
1976: }
1977: }
1978:
1.20 noro 1979: int dp_homogeneous(DP p)
1.15 noro 1980: {
1981: MP m;
1982: int d;
1983:
1984: if ( !p )
1985: return 1;
1986: else {
1987: m = BDY(p);
1988: d = m->dl->td;
1989: m = NEXT(m);
1990: for ( ; m; m = NEXT(m) ) {
1991: if ( m->dl->td != d )
1992: return 0;
1993: }
1994: return 1;
1995: }
1.16 noro 1996: }
1997:
1.20 noro 1998: void _print_mp(int nv,MP m)
1.16 noro 1999: {
2000: int i;
2001:
1.17 noro 2002: if ( !m )
1.16 noro 2003: return;
2004: for ( ; m; m = NEXT(m) ) {
2005: fprintf(stderr,"%d<",ITOS(C(m)));
2006: for ( i = 0; i < nv; i++ ) {
2007: fprintf(stderr,"%d",m->dl->d[i]);
2008: if ( i != nv-1 )
2009: fprintf(stderr," ");
2010: }
2011: fprintf(stderr,">",C(m));
2012: }
2013: fprintf(stderr,"\n");
1.15 noro 2014: }
1.26 noro 2015:
2016: static int cmp_mp_nvar;
2017:
2018: int comp_mp(MP *a,MP *b)
2019: {
2020: return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl);
2021: }
2022:
2023: void dp_sort(DP p,DP *rp)
2024: {
2025: MP t,mp,mp0;
2026: int i,n;
2027: DP r;
2028: MP *w;
2029:
2030: if ( !p ) {
2031: *rp = 0;
2032: return;
2033: }
2034: for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ );
2035: w = (MP *)ALLOCA(n*sizeof(MP));
2036: for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ )
2037: w[i] = t;
2038: cmp_mp_nvar = NV(p);
2039: qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp);
2040: mp0 = 0;
2041: for ( i = n-1; i >= 0; i-- ) {
2042: NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]);
2043: NEXT(mp) = mp0; mp0 = mp;
2044: }
2045: MKDP(p->nv,mp0,r);
2046: r->sugar = p->sugar;
2047: *rp = r;
2048: }
2049:
1.32 noro 2050: DP extract_initial_term_from_dp(DP p,int *weight,int n);
2051: LIST extract_initial_term(LIST f,int *weight,int n);
2052:
2053: DP extract_initial_term_from_dp(DP p,int *weight,int n)
2054: {
1.34 noro 2055: int w,t,i,top;
1.32 noro 2056: MP m,r0,r;
2057: DP dp;
2058:
2059: if ( !p ) return 0;
1.34 noro 2060: top = 1;
1.32 noro 2061: for ( m = BDY(p); m; m = NEXT(m) ) {
2062: for ( i = 0, t = 0; i < n; i++ )
2063: t += weight[i]*m->dl->d[i];
1.34 noro 2064: if ( top || t > w ) {
1.32 noro 2065: r0 = 0;
2066: w = t;
1.34 noro 2067: top = 0;
1.32 noro 2068: }
2069: if ( t == w ) {
2070: NEXTMP(r0,r);
2071: r->dl = m->dl;
2072: r->c = m->c;
2073: }
2074: }
2075: NEXT(r) = 0;
2076: MKDP(p->nv,r0,dp);
2077: return dp;
2078: }
2079:
2080: LIST extract_initial_term(LIST f,int *weight,int n)
2081: {
2082: NODE nd,r0,r;
2083: Obj p;
2084: LIST l;
2085:
2086: nd = BDY(f);
2087: for ( r0 = 0; nd; nd = NEXT(nd) ) {
2088: NEXTNODE(r0,r);
2089: p = (Obj)BDY(nd);
2090: BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n);
2091: }
2092: if ( r0 ) NEXT(r) = 0;
2093: MKLIST(l,r0);
2094: return l;
2095: }
2096:
2097: LIST dp_initial_term(LIST f,struct order_spec *ord)
2098: {
2099: int n,l,i;
2100: struct weight_or_block *worb;
2101: int *weight;
2102:
2103: switch ( ord->id ) {
2104: case 2: /* matrix order */
2105: /* extract the first row */
2106: n = ord->nv;
2107: weight = ord->ord.matrix.matrix[0];
2108: return extract_initial_term(f,weight,n);
2109: case 3: /* composite order */
2110: /* the first w_or_b */
2111: worb = ord->ord.composite.w_or_b;
2112: switch ( worb->type ) {
2113: case IS_DENSE_WEIGHT:
2114: n = worb->length;
2115: weight = worb->body.dense_weight;
2116: return extract_initial_term(f,weight,n);
2117: case IS_SPARSE_WEIGHT:
2118: n = ord->nv;
2119: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 2120: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 2121: l = worb->length;
2122: for ( i = 0; i < l; i++ )
2123: weight[worb->body.sparse_weight[i].pos]
2124: = worb->body.sparse_weight[i].value;
2125: return extract_initial_term(f,weight,n);
2126: default:
2127: error("dp_initial_term : unsupported order");
2128: }
2129: default:
2130: error("dp_initial_term : unsupported order");
2131: }
2132: }
2133:
2134: int highest_order_dp(DP p,int *weight,int n);
2135: LIST highest_order(LIST f,int *weight,int n);
2136:
2137: int highest_order_dp(DP p,int *weight,int n)
2138: {
1.34 noro 2139: int w,t,i,top;
1.32 noro 2140: MP m;
2141:
2142: if ( !p ) return -1;
1.34 noro 2143: top = 1;
1.32 noro 2144: for ( m = BDY(p); m; m = NEXT(m) ) {
2145: for ( i = 0, t = 0; i < n; i++ )
2146: t += weight[i]*m->dl->d[i];
1.34 noro 2147: if ( top || t > w ) {
1.32 noro 2148: w = t;
1.34 noro 2149: top = 0;
2150: }
1.32 noro 2151: }
2152: return w;
2153: }
2154:
2155: LIST highest_order(LIST f,int *weight,int n)
2156: {
2157: int h;
2158: NODE nd,r0,r;
2159: Obj p;
2160: LIST l;
2161: Q q;
2162:
2163: nd = BDY(f);
2164: for ( r0 = 0; nd; nd = NEXT(nd) ) {
2165: NEXTNODE(r0,r);
2166: p = (Obj)BDY(nd);
2167: h = highest_order_dp((DP)p,weight,n);
2168: STOQ(h,q);
2169: BDY(r) = (pointer)q;
2170: }
2171: if ( r0 ) NEXT(r) = 0;
2172: MKLIST(l,r0);
2173: return l;
2174: }
2175:
2176: LIST dp_order(LIST f,struct order_spec *ord)
2177: {
2178: int n,l,i;
2179: struct weight_or_block *worb;
2180: int *weight;
2181:
2182: switch ( ord->id ) {
2183: case 2: /* matrix order */
2184: /* extract the first row */
2185: n = ord->nv;
2186: weight = ord->ord.matrix.matrix[0];
2187: return highest_order(f,weight,n);
2188: case 3: /* composite order */
2189: /* the first w_or_b */
2190: worb = ord->ord.composite.w_or_b;
2191: switch ( worb->type ) {
2192: case IS_DENSE_WEIGHT:
2193: n = worb->length;
2194: weight = worb->body.dense_weight;
2195: return highest_order(f,weight,n);
2196: case IS_SPARSE_WEIGHT:
2197: n = ord->nv;
2198: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 2199: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 2200: l = worb->length;
2201: for ( i = 0; i < l; i++ )
2202: weight[worb->body.sparse_weight[i].pos]
2203: = worb->body.sparse_weight[i].value;
2204: return highest_order(f,weight,n);
2205: default:
2206: error("dp_initial_term : unsupported order");
2207: }
2208: default:
2209: error("dp_initial_term : unsupported order");
1.35 noro 2210: }
2211: }
2212:
2213: int dpv_ht(DPV p,DP *h)
2214: {
2215: int len,max,maxi,i,t;
2216: DP *e;
2217: MP m,mr;
2218:
2219: len = p->len;
2220: e = p->body;
2221: max = -1;
2222: maxi = -1;
2223: for ( i = 0; i < len; i++ )
2224: if ( e[i] && (t = BDY(e[i])->dl->td) > max ) {
2225: max = t;
2226: maxi = i;
2227: }
2228: if ( max < 0 ) {
2229: *h = 0;
2230: return -1;
2231: } else {
2232: m = BDY(e[maxi]);
2233: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
2234: MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td; /* XXX */
2235: return maxi;
1.32 noro 2236: }
2237: }