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