![[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.40 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.39 2005/08/25 18:59:11 ohara 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.13 noro 816: /* m-reduction over a field */
817:
1.20 noro 818: void dp_red_f(DP p1,DP p2,DP *rest)
1.13 noro 819: {
820: int i,n;
821: DL d1,d2,d;
822: MP m;
1.20 noro 823: DP t,s;
1.13 noro 824: Obj a,b;
825:
826: n = p1->nv;
827: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
828:
829: NEWDL(d,n); d->td = d1->td - d2->td;
830: for ( i = 0; i < n; i++ )
831: d->d[i] = d1->d[i]-d2->d[i];
832:
833: NEWMP(m); m->dl = d;
834: divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
835: C(m) = (P)b;
836: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
837:
838: muld(CO,s,p2,&t); addd(CO,p1,t,rest);
839: }
840:
1.20 noro 841: void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp)
1.7 noro 842: {
843: int i,n;
844: DL d1,d2,d;
845: MP m;
846: DP t,s,r,h;
847: P c1,c2,g,u;
848:
849: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
850: NEWDL(d,n); d->td = d1->td - d2->td;
851: for ( i = 0; i < n; i++ )
852: d->d[i] = d1->d[i]-d2->d[i];
853: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
854: gcdprsmp(CO,mod,c1,c2,&g);
855: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
856: if ( NUM(c2) ) {
857: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
858: }
859: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1.11 noro 860: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7 noro 861: if ( NUM(c2) ) {
862: addmd(CO,mod,p1,t,&r); h = p0;
863: } else {
864: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
865: }
866: *head = h; *rest = r; *dnp = c2;
867: }
868:
1.10 noro 869: struct oEGT eg_red_mod;
870:
1.20 noro 871: void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp)
1.7 noro 872: {
873: int i,n;
874: DL d1,d2,d;
875: MP m;
876: DP t,s;
1.16 noro 877: int c,c1,c2;
878: extern int do_weyl;
1.7 noro 879:
880: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
881: _NEWDL(d,n); d->td = d1->td - d2->td;
882: for ( i = 0; i < n; i++ )
883: d->d[i] = d1->d[i]-d2->d[i];
1.16 noro 884: c = invm(ITOS(BDY(p2)->c),mod);
885: c2 = ITOS(BDY(p1)->c);
886: DMAR(c,c2,0,mod,c1);
1.7 noro 887: _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
1.16 noro 888: #if 0
1.7 noro 889: _MKDP(n,m,s); s->sugar = d->td;
890: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 891: #else
892: if ( do_weyl ) {
1.19 noro 893: _MKDP(n,m,s); s->sugar = d->td;
894: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 895: } else {
896: _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
897: }
898: #endif
1.10 noro 899: /* get_eg(&t0); */
1.7 noro 900: _addmd_destructive(mod,p1,t,rp);
1.10 noro 901: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 902: }
903:
904: /*
905: * normal form computation
906: *
907: */
1.5 noro 908:
1.20 noro 909: void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp)
1.5 noro 910: {
911: DP u,p,d,s,t,dmy;
912: NODE l;
913: MP m,mr;
914: int i,n;
915: int *wb;
916: int sugar,psugar;
917: P dn,tdn,tdn1;
918:
919: dn = (P)ONE;
920: if ( !g ) {
921: *rp = 0; *dnp = dn; return;
922: }
923: for ( n = 0, l = b; l; l = NEXT(l), n++ );
924: wb = (int *)ALLOCA(n*sizeof(int));
925: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
926: wb[i] = QTOS((Q)BDY(l));
927: sugar = g->sugar;
928: for ( d = 0; g; ) {
929: for ( u = 0, i = 0; i < n; i++ ) {
930: if ( dp_redble(g,p = ps[wb[i]]) ) {
931: dp_red(d,g,p,&t,&u,&tdn,&dmy);
932: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
933: sugar = MAX(sugar,psugar);
934: if ( !u ) {
935: if ( d )
936: d->sugar = sugar;
937: *rp = d; *dnp = dn; return;
938: } else {
939: d = t;
940: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
941: }
942: break;
943: }
944: }
945: if ( u )
946: g = u;
947: else if ( !full ) {
948: if ( g ) {
949: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
950: }
951: *rp = g; *dnp = dn; return;
952: } else {
953: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
954: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
955: addd(CO,d,t,&s); d = s;
956: dp_rest(g,&t); g = t;
957: }
958: }
959: if ( d )
960: d->sugar = sugar;
961: *rp = d; *dnp = dn;
962: }
963:
1.13 noro 964: /* nf computation over Z */
965:
1.20 noro 966: void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
1.5 noro 967: {
968: DP u,p,d,s,t,dmy1;
969: P dmy;
970: NODE l;
971: MP m,mr;
972: int i,n;
973: int *wb;
974: int hmag;
975: int sugar,psugar;
976:
977: if ( !g ) {
978: *rp = 0; return;
979: }
980: for ( n = 0, l = b; l; l = NEXT(l), n++ );
981: wb = (int *)ALLOCA(n*sizeof(int));
982: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
983: wb[i] = QTOS((Q)BDY(l));
1.12 noro 984:
1.13 noro 985: hmag = multiple*HMAG(g);
1.5 noro 986: sugar = g->sugar;
1.12 noro 987:
1.5 noro 988: for ( d = 0; g; ) {
989: for ( u = 0, i = 0; i < n; i++ ) {
990: if ( dp_redble(g,p = ps[wb[i]]) ) {
991: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
992: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
993: sugar = MAX(sugar,psugar);
994: if ( !u ) {
995: if ( d )
996: d->sugar = sugar;
997: *rp = d; return;
998: }
999: d = t;
1000: break;
1001: }
1002: }
1003: if ( u ) {
1004: g = u;
1005: if ( d ) {
1.13 noro 1006: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 1007: dp_ptozp2(d,g,&t,&u); d = t; g = u;
1008: hmag = multiple*HMAG(d);
1009: }
1010: } else {
1.13 noro 1011: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 1012: dp_ptozp(g,&t); g = t;
1013: hmag = multiple*HMAG(g);
1014: }
1015: }
1016: }
1017: else if ( !full ) {
1018: if ( g ) {
1019: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1020: }
1021: *rp = g; return;
1022: } else {
1023: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1024: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1025: addd(CO,d,t,&s); d = s;
1026: dp_rest(g,&t); g = t;
1027:
1028: }
1029: }
1030: if ( d )
1031: d->sugar = sugar;
1032: *rp = d;
1033: }
1034:
1.13 noro 1035: /* nf computation over a field */
1036:
1.20 noro 1037: void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp)
1.13 noro 1038: {
1039: DP u,p,d,s,t;
1040: NODE l;
1041: MP m,mr;
1042: int i,n;
1043: int *wb;
1044: int sugar,psugar;
1045:
1046: if ( !g ) {
1047: *rp = 0; return;
1048: }
1049: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1050: wb = (int *)ALLOCA(n*sizeof(int));
1051: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1052: wb[i] = QTOS((Q)BDY(l));
1053:
1054: sugar = g->sugar;
1055: for ( d = 0; g; ) {
1056: for ( u = 0, i = 0; i < n; i++ ) {
1057: if ( dp_redble(g,p = ps[wb[i]]) ) {
1058: dp_red_f(g,p,&u);
1059: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1060: sugar = MAX(sugar,psugar);
1061: if ( !u ) {
1062: if ( d )
1063: d->sugar = sugar;
1064: *rp = d; return;
1065: }
1066: break;
1067: }
1068: }
1069: if ( u )
1070: g = u;
1071: else if ( !full ) {
1072: if ( g ) {
1073: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1074: }
1075: *rp = g; return;
1076: } else {
1077: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1078: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1079: addd(CO,d,t,&s); d = s;
1080: dp_rest(g,&t); g = t;
1081: }
1082: }
1083: if ( d )
1084: d->sugar = sugar;
1085: *rp = d;
1086: }
1087:
1088: /* nf computation over GF(mod) (only for internal use) */
1089:
1.20 noro 1090: void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1091: {
1092: DP u,p,d,s,t;
1093: P dmy;
1094: NODE l;
1095: MP m,mr;
1096: int sugar,psugar;
1097:
1098: if ( !g ) {
1099: *rp = 0; return;
1100: }
1101: sugar = g->sugar;
1102: for ( d = 0; g; ) {
1103: for ( u = 0, l = b; l; l = NEXT(l) ) {
1104: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1105: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
1106: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1107: sugar = MAX(sugar,psugar);
1108: if ( !u ) {
1109: if ( d )
1110: d->sugar = sugar;
1111: *rp = d; return;
1112: }
1113: d = t;
1114: break;
1115: }
1116: }
1117: if ( u )
1118: g = u;
1119: else if ( !full ) {
1120: if ( g ) {
1121: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1122: }
1123: *rp = g; return;
1124: } else {
1125: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1126: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1127: addmd(CO,mod,d,t,&s); d = s;
1128: dp_rest(g,&t); g = t;
1129: }
1130: }
1131: if ( d )
1132: d->sugar = sugar;
1133: *rp = d;
1134: }
1135:
1.20 noro 1136: void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp)
1.5 noro 1137: {
1138: DP u,p,d,s,t;
1139: NODE l;
1140: MP m,mr;
1141: int i,n;
1142: int *wb;
1143: int sugar,psugar;
1144: P dn,tdn,tdn1;
1145:
1146: dn = (P)ONEM;
1147: if ( !g ) {
1148: *rp = 0; *dnp = dn; return;
1149: }
1150: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1151: wb = (int *)ALLOCA(n*sizeof(int));
1152: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1153: wb[i] = QTOS((Q)BDY(l));
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_mod(d,g,p,mod,&t,&u,&tdn);
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; *dnp = dn; return;
1165: } else {
1166: d = t;
1167: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1168: }
1169: break;
1170: }
1171: }
1172: if ( u )
1173: g = u;
1174: else if ( !full ) {
1175: if ( g ) {
1176: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1177: }
1178: *rp = g; *dnp = dn; return;
1179: } else {
1180: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1181: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1182: addmd(CO,mod,d,t,&s); d = s;
1183: dp_rest(g,&t); g = t;
1184: }
1185: }
1186: if ( d )
1187: d->sugar = sugar;
1188: *rp = d; *dnp = dn;
1189: }
1190:
1.20 noro 1191: void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1192: {
1.20 noro 1193: DP u,p,d;
1.7 noro 1194: NODE l;
1.20 noro 1195: MP m,mrd;
1196: int sugar,psugar,n,h_reducible;
1.5 noro 1197:
1.7 noro 1198: if ( !g ) {
1199: *rp = 0; return;
1.5 noro 1200: }
1.7 noro 1201: sugar = g->sugar;
1202: n = g->nv;
1203: for ( d = 0; g; ) {
1204: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
1205: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1206: h_reducible = 1;
1207: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1208: _dp_red_mod_destructive(g,p,mod,&u); g = u;
1209: sugar = MAX(sugar,psugar);
1210: if ( !g ) {
1211: if ( d )
1212: d->sugar = sugar;
1213: _dptodp(d,rp); _free_dp(d); return;
1214: }
1215: break;
1216: }
1217: }
1218: if ( !h_reducible ) {
1219: /* head term is not reducible */
1220: if ( !full ) {
1221: if ( g )
1222: g->sugar = sugar;
1223: _dptodp(g,rp); _free_dp(g); return;
1224: } else {
1225: m = BDY(g);
1226: if ( NEXT(m) ) {
1227: BDY(g) = NEXT(m); NEXT(m) = 0;
1228: } else {
1229: _FREEDP(g); g = 0;
1230: }
1231: if ( d ) {
1232: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1233: NEXT(mrd) = m;
1234: } else {
1235: _MKDP(n,m,d);
1236: }
1237: }
1238: }
1.5 noro 1239: }
1.7 noro 1240: if ( d )
1241: d->sugar = sugar;
1242: _dptodp(d,rp); _free_dp(d);
1.5 noro 1243: }
1.13 noro 1244:
1245: /* reduction by linear base over a field */
1246:
1.20 noro 1247: void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p)
1.13 noro 1248: {
1249: DP r1,r2,b1,b2,t,s;
1250: Obj c,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: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
1266: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
1267: muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
1268: muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
1269: }
1270: }
1271: *r1p = r1; *r2p = r2;
1272: }
1273:
1274: /* reduction by linear base over GF(mod) */
1.5 noro 1275:
1.20 noro 1276: void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p)
1.5 noro 1277: {
1.7 noro 1278: DP r1,r2,b1,b2,t,s;
1279: P c;
1280: MQ c1,c2;
1281: NODE l,b;
1282: int n;
1283:
1284: if ( !p1 ) {
1285: *r1p = p1; *r2p = p2; return;
1286: }
1287: n = p1->nv;
1288: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1289: if ( !r1 ) {
1290: *r1p = r1; *r2p = r2; return;
1291: }
1292: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1293: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1294: b2 = (DP)BDY(NEXT(b));
1295: invmq(mod,(MQ)BDY(b1)->c,&c1);
1296: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
1297: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
1298: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
1299: }
1300: }
1301: *r1p = r1; *r2p = r2;
1.5 noro 1302: }
1303:
1.20 noro 1304: void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp)
1.5 noro 1305: {
1.7 noro 1306: DP s,t,u;
1307: MP m;
1308: DL h;
1309: int i,n;
1310:
1311: if ( !p ) {
1312: *rp = p; return;
1313: }
1314: n = p->nv;
1315: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1316: h = m->dl;
1317: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1318: i++;
1319: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1320: addmd(CO,mod,s,t,&u); s = u;
1.24 noro 1321: }
1322: *rp = s;
1323: }
1324:
1325: void dp_nf_tab_f(DP p,LIST *tab,DP *rp)
1326: {
1327: DP s,t,u;
1328: MP m;
1329: DL h;
1330: int i,n;
1331:
1332: if ( !p ) {
1333: *rp = p; return;
1334: }
1335: n = p->nv;
1336: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1337: h = m->dl;
1338: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1339: i++;
1340: muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1341: addd(CO,s,t,&u); s = u;
1.7 noro 1342: }
1343: *rp = s;
1.5 noro 1344: }
1345:
1.7 noro 1346: /*
1347: * setting flags
1.30 noro 1348: * call create_order_spec with vl=0 to set old type order.
1.7 noro 1349: *
1350: */
1351:
1.27 noro 1352: int create_order_spec(VL vl,Obj obj,struct order_spec **specp)
1.5 noro 1353: {
1.37 noro 1354: int i,j,n,s,row,col,ret;
1.27 noro 1355: struct order_spec *spec;
1.7 noro 1356: struct order_pair *l;
1357: NODE node,t,tn;
1358: MAT m;
1359: pointer **b;
1360: int **w;
1.5 noro 1361:
1.37 noro 1362: if ( vl && obj && OID(obj) == O_LIST ) {
1363: ret = create_composite_order_spec(vl,(LIST)obj,specp);
1364: if ( show_orderspec )
1365: print_composite_order_spec(*specp);
1366: return ret;
1367: }
1.27 noro 1368:
1369: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 1370: if ( !obj || NUM(obj) ) {
1371: spec->id = 0; spec->obj = obj;
1372: spec->ord.simple = QTOS((Q)obj);
1373: return 1;
1374: } else if ( OID(obj) == O_LIST ) {
1375: node = BDY((LIST)obj);
1376: for ( n = 0, t = node; t; t = NEXT(t), n++ );
1377: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1378: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
1379: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
1380: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
1381: s += l[i].length;
1382: }
1383: spec->id = 1; spec->obj = obj;
1384: spec->ord.block.order_pair = l;
1385: spec->ord.block.length = n; spec->nv = s;
1386: return 1;
1387: } else if ( OID(obj) == O_MAT ) {
1388: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
1389: w = almat(row,col);
1390: for ( i = 0; i < row; i++ )
1391: for ( j = 0; j < col; j++ )
1392: w[i][j] = QTOS((Q)b[i][j]);
1393: spec->id = 2; spec->obj = obj;
1394: spec->nv = col; spec->ord.matrix.row = row;
1395: spec->ord.matrix.matrix = w;
1396: return 1;
1397: } else
1.5 noro 1398: return 0;
1399: }
1400:
1.28 noro 1401: void print_composite_order_spec(struct order_spec *spec)
1402: {
1403: int nv,n,len,i,j,k,start;
1404: struct weight_or_block *worb;
1405:
1406: nv = spec->nv;
1407: n = spec->ord.composite.length;
1408: worb = spec->ord.composite.w_or_b;
1409: for ( i = 0; i < n; i++, worb++ ) {
1410: len = worb->length;
1411: printf("[ ");
1412: switch ( worb->type ) {
1413: case IS_DENSE_WEIGHT:
1414: for ( j = 0; j < len; j++ )
1415: printf("%d ",worb->body.dense_weight[j]);
1416: for ( ; j < nv; j++ )
1417: printf("0 ");
1418: break;
1419: case IS_SPARSE_WEIGHT:
1420: for ( j = 0, k = 0; j < nv; j++ )
1421: if ( j == worb->body.sparse_weight[k].pos )
1422: printf("%d ",worb->body.sparse_weight[k++].value);
1423: else
1424: printf("0 ");
1425: break;
1426: case IS_BLOCK:
1427: start = worb->body.block.start;
1428: for ( j = 0; j < start; j++ ) printf("0 ");
1429: switch ( worb->body.block.order ) {
1430: case 0:
1431: for ( k = 0; k < len; k++, j++ ) printf("R ");
1432: break;
1433: case 1:
1434: for ( k = 0; k < len; k++, j++ ) printf("G ");
1435: break;
1436: case 2:
1437: for ( k = 0; k < len; k++, j++ ) printf("L ");
1438: break;
1439: }
1440: for ( ; j < nv; j++ ) printf("0 ");
1441: break;
1442: }
1443: printf("]\n");
1444: }
1.38 noro 1445: }
1446:
1447: struct order_spec *append_block(struct order_spec *spec,
1448: int nv,int nalg,int ord)
1449: {
1450: MAT m,mat;
1451: int i,j,row,col,n;
1452: Q **b,**wp;
1453: int **w;
1454: NODE t,s,s0;
1455: struct order_pair *l,*l0;
1456: int n0,nv0;
1457: LIST list0,list1,list;
1458: Q oq,nq;
1459: struct order_spec *r;
1460:
1461: r = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1462: switch ( spec->id ) {
1463: case 0:
1464: STOQ(spec->ord.simple,oq); STOQ(nv,nq);
1465: t = mknode(2,oq,nq); MKLIST(list0,t);
1466: STOQ(ord,oq); STOQ(nalg,nq);
1467: t = mknode(2,oq,nq); MKLIST(list1,t);
1468: t = mknode(2,list0,list1); MKLIST(list,t);
1469: l = (struct order_pair *)MALLOC_ATOMIC(2*sizeof(struct order_pair));
1470: l[0].order = spec->ord.simple; l[0].length = nv;
1471: l[1].order = ord; l[1].length = nalg;
1472: r->id = 1; r->obj = (Obj)list;
1473: r->ord.block.order_pair = l;
1474: r->ord.block.length = 2;
1475: r->nv = nv+nalg;
1476: break;
1477: case 1:
1478: if ( spec->nv != nv )
1479: error("append_block : number of variables mismatch");
1480: l0 = spec->ord.block.order_pair;
1481: n0 = spec->ord.block.length;
1482: nv0 = spec->nv;
1483: list0 = (LIST)spec->obj;
1484: n = n0+1;
1485: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1486: for ( i = 0; i < n0; i++ )
1487: l[i] = l0[i];
1488: l[i].order = ord; l[i].length = nalg;
1489: for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) {
1490: NEXTNODE(s0,s); BDY(s) = BDY(t);
1491: }
1492: STOQ(ord,oq); STOQ(nalg,nq);
1493: t = mknode(2,oq,nq); MKLIST(list,t);
1494: NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0;
1495: MKLIST(list,s0);
1496: r->id = 1; r->obj = (Obj)list;
1497: r->ord.block.order_pair = l;
1498: r->ord.block.length = n;
1499: r->nv = nv+nalg;
1500: break;
1501: case 2:
1502: if ( spec->nv != nv )
1503: error("append_block : number of variables mismatch");
1504: m = (MAT)spec->obj;
1505: row = m->row; col = m->col; b = (Q **)BDY(m);
1506: w = almat(row+nalg,col+nalg);
1507: MKMAT(mat,row+nalg,col+nalg); wp = (Q **)BDY(mat);
1508: for ( i = 0; i < row; i++ )
1509: for ( j = 0; j < col; j++ ) {
1510: w[i][j] = QTOS(b[i][j]);
1511: wp[i][j] = b[i][j];
1512: }
1513: for ( i = 0; i < nalg; i++ ) {
1514: w[i+row][i+col] = 1;
1515: wp[i+row][i+col] = ONE;
1516: }
1517: r->id = 2; r->obj = (Obj)mat;
1518: r->nv = col+nalg; r->ord.matrix.row = row+nalg;
1519: r->ord.matrix.matrix = w;
1520: break;
1521: case 3:
1522: default:
1523: /* XXX */
1524: error("append_block : not implemented yet");
1525: }
1526: return r;
1.28 noro 1527: }
1528:
1.37 noro 1529: int comp_sw(struct sparse_weight *a, struct sparse_weight *b)
1530: {
1531: if ( a->pos > b->pos ) return 1;
1532: else if ( a->pos < b->pos ) return -1;
1533: else return 0;
1534: }
1535:
1.27 noro 1536: /* order = [w_or_b, w_or_b, ... ] */
1537: /* w_or_b = w or b */
1538: /* w = [1,2,...] or [x,1,y,2,...] */
1539: /* b = [@lex,x,y,...,z] etc */
1540:
1541: int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp)
1542: {
1543: NODE wb,t,p;
1544: struct order_spec *spec;
1545: VL tvl;
1.29 noro 1546: int n,i,j,k,l,start,end,len,w;
1.27 noro 1547: int *dw;
1548: struct sparse_weight *sw;
1549: struct weight_or_block *w_or_b;
1550: Obj a0;
1551: NODE a;
1.29 noro 1552: V v,sv,ev;
1553: SYMBOL sym;
1554: int *top;
1.27 noro 1555:
1556: /* l = number of vars in vl */
1557: for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ );
1558: /* n = number of primitives in order */
1559: wb = BDY(order);
1560: n = length(wb);
1561: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1562: spec->id = 3;
1563: spec->obj = (Obj)order;
1564: spec->nv = l;
1565: spec->ord.composite.length = n;
1.28 noro 1566: w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *)
1.29 noro 1567: MALLOC(sizeof(struct weight_or_block)*(n+1));
1568:
1569: /* top : register the top variable in each w_or_b specification */
1570: top = (int *)ALLOCA(l*sizeof(int));
1571: for ( i = 0; i < l; i++ ) top[i] = 0;
1572:
1.28 noro 1573: for ( t = wb, i = 0; t; t = NEXT(t), i++ ) {
1.30 noro 1574: if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST )
1575: error("a list of lists must be specified for the key \"order\"");
1.28 noro 1576: a = BDY((LIST)BDY(t));
1.27 noro 1577: len = length(a);
1578: a0 = (Obj)BDY(a);
1579: if ( !a0 || OID(a0) == O_N ) {
1.28 noro 1580: /* a is a dense weight vector */
1.27 noro 1581: dw = (int *)MALLOC(sizeof(int)*len);
1.30 noro 1582: for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) {
1583: if ( !INT((Q)BDY(p)) )
1584: error("a dense weight vector must be specified as a list of integers");
1.27 noro 1585: dw[j] = QTOS((Q)BDY(p));
1.30 noro 1586: }
1.27 noro 1587: w_or_b[i].type = IS_DENSE_WEIGHT;
1588: w_or_b[i].length = len;
1589: w_or_b[i].body.dense_weight = dw;
1.29 noro 1590:
1591: /* find the top */
1592: for ( k = 0; k < len && !dw[k]; k++ );
1593: if ( k < len ) top[k] = 1;
1594:
1.27 noro 1595: } else if ( OID(a0) == O_P ) {
1.28 noro 1596: /* a is a sparse weight vector */
1597: len >>= 1;
1.27 noro 1598: sw = (struct sparse_weight *)
1599: MALLOC(sizeof(struct sparse_weight)*len);
1600: for ( j = 0, p = a; j < len; j++ ) {
1.30 noro 1601: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
1602: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 1603: v = VR((P)BDY(p)); p = NEXT(p);
1.27 noro 1604: for ( tvl = vl, k = 0; tvl && tvl->v != v;
1605: k++, tvl = NEXT(tvl) );
1606: if ( !tvl )
1.30 noro 1607: error("invalid variable name in a sparse weight vector");
1.27 noro 1608: sw[j].pos = k;
1.30 noro 1609: if ( !INT((Q)BDY(p)) )
1610: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 1611: sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p);
1.27 noro 1612: }
1.37 noro 1613: qsort(sw,len,sizeof(struct sparse_weight),
1614: (int (*)(const void *,const void *))comp_sw);
1.27 noro 1615: w_or_b[i].type = IS_SPARSE_WEIGHT;
1616: w_or_b[i].length = len;
1617: w_or_b[i].body.sparse_weight = sw;
1.29 noro 1618:
1619: /* find the top */
1620: for ( k = 0; k < len && !sw[k].value; k++ );
1621: if ( k < len ) top[sw[k].pos] = 1;
1622: } else if ( OID(a0) == O_RANGE ) {
1623: /* [range(v1,v2),w] */
1624: sv = VR((P)(((RANGE)a0)->start));
1625: ev = VR((P)(((RANGE)a0)->end));
1626: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
1627: if ( !tvl )
1628: error("invalid range");
1629: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
1630: if ( !tvl )
1631: error("invalid range");
1632: len = end-start+1;
1633: sw = (struct sparse_weight *)
1634: MALLOC(sizeof(struct sparse_weight)*len);
1635: w = QTOS((Q)BDY(NEXT(a)));
1636: for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) );
1637: for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) {
1638: sw[j].pos = k;
1639: sw[j].value = w;
1640: }
1641: w_or_b[i].type = IS_SPARSE_WEIGHT;
1642: w_or_b[i].length = len;
1643: w_or_b[i].body.sparse_weight = sw;
1644:
1645: /* register the top */
1646: if ( w ) top[start] = 1;
1.28 noro 1647: } else if ( OID(a0) == O_SYMBOL ) {
1648: /* a is a block */
1.29 noro 1649: sym = (SYMBOL)a0; a = NEXT(a); len--;
1650: if ( OID((Obj)BDY(a)) == O_RANGE ) {
1651: sv = VR((P)(((RANGE)BDY(a))->start));
1652: ev = VR((P)(((RANGE)BDY(a))->end));
1653: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
1654: if ( !tvl )
1655: error("invalid range");
1656: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
1657: if ( !tvl )
1658: error("invalid range");
1659: len = end-start+1;
1660: } else {
1661: for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a));
1.28 noro 1662: tvl = NEXT(tvl), start++ );
1.29 noro 1663: for ( p = NEXT(a), tvl = NEXT(tvl); p;
1.30 noro 1664: p = NEXT(p), tvl = NEXT(tvl) ) {
1665: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
1666: error("a block must be specified as [ordsymbol,var1,var2,...]");
1.29 noro 1667: if ( tvl->v != VR((P)BDY(p)) ) break;
1.30 noro 1668: }
1.29 noro 1669: if ( p )
1.30 noro 1670: error("a block must be contiguous in the variable list");
1.29 noro 1671: }
1.28 noro 1672: w_or_b[i].type = IS_BLOCK;
1673: w_or_b[i].length = len;
1674: w_or_b[i].body.block.start = start;
1675: if ( !strcmp(sym->name,"@grlex") )
1676: w_or_b[i].body.block.order = 0;
1677: else if ( !strcmp(sym->name,"@glex") )
1678: w_or_b[i].body.block.order = 1;
1679: else if ( !strcmp(sym->name,"@lex") )
1680: w_or_b[i].body.block.order = 2;
1681: else
1.29 noro 1682: error("invalid ordername");
1683: /* register the tops */
1684: for ( j = 0, k = start; j < len; j++, k++ )
1685: top[k] = 1;
1.28 noro 1686: }
1.29 noro 1687: }
1688: for ( k = 0; k < l && top[k]; k++ );
1689: if ( k < l ) {
1690: /* incomplete order specification; add @grlex */
1691: w_or_b[n].type = IS_BLOCK;
1692: w_or_b[n].length = l;
1693: w_or_b[n].body.block.start = 0;
1694: w_or_b[n].body.block.order = 0;
1695: spec->ord.composite.length = n+1;
1.27 noro 1696: }
1697: }
1698:
1.35 noro 1699: /* module order spec */
1700:
1701: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s)
1702: {
1703: struct modorder_spec *spec;
1704: NODE n,t;
1705: LIST list;
1706: int *ds;
1707: int i,l;
1708: Q q;
1709:
1710: *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec));
1711: spec->id = id;
1712: if ( shift ) {
1713: n = BDY(shift);
1714: spec->len = l = length(n);
1715: spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int));
1716: for ( t = n, i = 0; t; t = NEXT(t), i++ )
1717: ds[i] = QTOS((Q)BDY(t));
1718: } else {
1719: spec->len = 0;
1720: spec->degree_shift = 0;
1721: }
1722: STOQ(id,q);
1723: n = mknode(2,q,shift);
1724: MKLIST(list,n);
1725: spec->obj = (Obj)list;
1726: }
1727:
1.7 noro 1728: /*
1729: * converters
1730: *
1731: */
1732:
1.20 noro 1733: void dp_homo(DP p,DP *rp)
1.5 noro 1734: {
1.7 noro 1735: MP m,mr,mr0;
1736: int i,n,nv,td;
1737: DL dl,dlh;
1.5 noro 1738:
1.7 noro 1739: if ( !p )
1740: *rp = 0;
1741: else {
1742: n = p->nv; nv = n + 1;
1743: m = BDY(p); td = sugard(m);
1744: for ( mr0 = 0; m; m = NEXT(m) ) {
1745: NEXTMP(mr0,mr); mr->c = m->c;
1746: dl = m->dl;
1747: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1748: dlh->td = td;
1749: for ( i = 0; i < n; i++ )
1750: dlh->d[i] = dl->d[i];
1751: dlh->d[n] = td - dl->td;
1752: }
1753: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 1754: }
1755: }
1756:
1.20 noro 1757: void dp_dehomo(DP p,DP *rp)
1.5 noro 1758: {
1.7 noro 1759: MP m,mr,mr0;
1760: int i,n,nv;
1761: DL dl,dlh;
1.5 noro 1762:
1.7 noro 1763: if ( !p )
1764: *rp = 0;
1765: else {
1766: n = p->nv; nv = n - 1;
1767: m = BDY(p);
1768: for ( mr0 = 0; m; m = NEXT(m) ) {
1769: NEXTMP(mr0,mr); mr->c = m->c;
1770: dlh = m->dl;
1771: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1772: dl->td = dlh->td - dlh->d[nv];
1773: for ( i = 0; i < nv; i++ )
1774: dl->d[i] = dlh->d[i];
1775: }
1776: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1777: }
1.5 noro 1778: }
1779:
1.20 noro 1780: void dp_mod(DP p,int mod,NODE subst,DP *rp)
1.5 noro 1781: {
1.7 noro 1782: MP m,mr,mr0;
1783: P t,s,s1;
1784: V v;
1785: NODE tn;
1.5 noro 1786:
1.7 noro 1787: if ( !p )
1788: *rp = 0;
1789: else {
1790: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1791: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
1792: v = VR((P)BDY(tn)); tn = NEXT(tn);
1793: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
1794: }
1795: ptomp(mod,s,&t);
1796: if ( t ) {
1797: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
1798: }
1799: }
1800: if ( mr0 ) {
1801: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1802: } else
1803: *rp = 0;
1804: }
1.5 noro 1805: }
1806:
1.20 noro 1807: void dp_rat(DP p,DP *rp)
1.5 noro 1808: {
1.7 noro 1809: MP m,mr,mr0;
1.5 noro 1810:
1.7 noro 1811: if ( !p )
1812: *rp = 0;
1813: else {
1814: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1815: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5 noro 1816: }
1.7 noro 1817: if ( mr0 ) {
1818: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1819: } else
1820: *rp = 0;
1.5 noro 1821: }
1822: }
1823:
1824:
1.27 noro 1825: void homogenize_order(struct order_spec *old,int n,struct order_spec **newp)
1.5 noro 1826: {
1.7 noro 1827: struct order_pair *l;
1828: int length,nv,row,i,j;
1829: int **newm,**oldm;
1.27 noro 1830: struct order_spec *new;
1.31 noro 1831: int onv,nnv,nlen,olen,owlen;
1832: struct weight_or_block *owb,*nwb;
1.5 noro 1833:
1.27 noro 1834: *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 1835: switch ( old->id ) {
1836: case 0:
1837: switch ( old->ord.simple ) {
1838: case 0:
1839: new->id = 0; new->ord.simple = 0; break;
1840: case 1:
1841: l = (struct order_pair *)
1842: MALLOC_ATOMIC(2*sizeof(struct order_pair));
1843: l[0].length = n; l[0].order = 1;
1844: l[1].length = 1; l[1].order = 2;
1845: new->id = 1;
1846: new->ord.block.order_pair = l;
1847: new->ord.block.length = 2; new->nv = n+1;
1848: break;
1849: case 2:
1850: new->id = 0; new->ord.simple = 1; break;
1851: case 3: case 4: case 5:
1852: new->id = 0; new->ord.simple = old->ord.simple+3;
1853: dp_nelim = n-1; break;
1854: case 6: case 7: case 8: case 9:
1855: new->id = 0; new->ord.simple = old->ord.simple; break;
1856: default:
1857: error("homogenize_order : invalid input");
1858: }
1859: break;
1860: case 1:
1861: length = old->ord.block.length;
1862: l = (struct order_pair *)
1863: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
1864: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
1865: l[length].order = 2; l[length].length = 1;
1866: new->id = 1; new->nv = n+1;
1867: new->ord.block.order_pair = l;
1868: new->ord.block.length = length+1;
1869: break;
1870: case 2:
1871: nv = old->nv; row = old->ord.matrix.row;
1872: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
1873: for ( i = 0; i <= nv; i++ )
1874: newm[0][i] = 1;
1875: for ( i = 0; i < row; i++ ) {
1876: for ( j = 0; j < nv; j++ )
1877: newm[i+1][j] = oldm[i][j];
1878: newm[i+1][j] = 0;
1879: }
1880: new->id = 2; new->nv = nv+1;
1881: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1.31 noro 1882: break;
1883: case 3:
1884: onv = old->nv;
1885: nnv = onv+1;
1886: olen = old->ord.composite.length;
1887: nlen = olen+1;
1888: owb = old->ord.composite.w_or_b;
1889: nwb = (struct weight_or_block *)
1890: MALLOC(nlen*sizeof(struct weight_or_block));
1891: for ( i = 0; i < olen; i++ ) {
1892: nwb[i].type = owb[i].type;
1893: switch ( owb[i].type ) {
1894: case IS_DENSE_WEIGHT:
1895: owlen = owb[i].length;
1896: nwb[i].length = owlen+1;
1897: nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int));
1898: for ( j = 0; j < owlen; j++ )
1899: nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j];
1900: nwb[i].body.dense_weight[owlen] = 0;
1901: break;
1902: case IS_SPARSE_WEIGHT:
1903: nwb[i].length = owb[i].length;
1904: nwb[i].body.sparse_weight = owb[i].body.sparse_weight;
1905: break;
1906: case IS_BLOCK:
1907: nwb[i].length = owb[i].length;
1908: nwb[i].body.block = owb[i].body.block;
1909: break;
1910: }
1911: }
1912: nwb[i].type = IS_SPARSE_WEIGHT;
1913: nwb[i].body.sparse_weight =
1914: (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
1915: nwb[i].body.sparse_weight[0].pos = onv;
1916: nwb[i].body.sparse_weight[0].value = 1;
1917: new->id = 3;
1918: new->nv = nnv;
1919: new->ord.composite.length = nlen;
1920: new->ord.composite.w_or_b = nwb;
1921: print_composite_order_spec(new);
1.7 noro 1922: break;
1923: default:
1924: error("homogenize_order : invalid input");
1.5 noro 1925: }
1.7 noro 1926: }
1927:
1.20 noro 1928: void qltozl(Q *w,int n,Q *dvr)
1.7 noro 1929: {
1930: N nm,dn;
1931: N g,l1,l2,l3;
1932: Q c,d;
1933: int i;
1934: struct oVECT v;
1.5 noro 1935:
1936: for ( i = 0; i < n; i++ )
1.7 noro 1937: if ( w[i] && !INT(w[i]) )
1938: break;
1939: if ( i == n ) {
1940: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
1941: igcdv(&v,dvr); return;
1942: }
1943: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
1944: for ( i = 1; i < n; i++ ) {
1945: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
1946: gcdn(nm,NM(c),&g); nm = g;
1947: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 1948: }
1.7 noro 1949: if ( UNIN(dn) )
1950: NTOQ(nm,1,d);
1951: else
1952: NDTOQ(nm,dn,1,d);
1953: *dvr = d;
1954: }
1.5 noro 1955:
1.20 noro 1956: int comp_nm(Q *a,Q *b)
1.7 noro 1957: {
1958: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
1959: }
1960:
1.20 noro 1961: void sortbynm(Q *w,int n)
1.7 noro 1962: {
1963: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
1964: }
1.5 noro 1965:
1966:
1.7 noro 1967: /*
1968: * simple operations
1969: *
1970: */
1.5 noro 1971:
1.20 noro 1972: int dp_redble(DP p1,DP p2)
1.7 noro 1973: {
1974: int i,n;
1975: DL d1,d2;
1.5 noro 1976:
1.7 noro 1977: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1978: if ( d1->td < d2->td )
1979: return 0;
1980: else {
1981: for ( i = 0, n = p1->nv; i < n; i++ )
1982: if ( d1->d[i] < d2->d[i] )
1983: return 0;
1984: return 1;
1.5 noro 1985: }
1986: }
1987:
1.20 noro 1988: void dp_subd(DP p1,DP p2,DP *rp)
1.5 noro 1989: {
1.7 noro 1990: int i,n;
1.5 noro 1991: DL d1,d2,d;
1992: MP m;
1.7 noro 1993: DP s;
1.5 noro 1994:
1995: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 noro 1996: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 1997: for ( i = 0; i < n; i++ )
1.7 noro 1998: d->d[i] = d1->d[i]-d2->d[i];
1999: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
2000: MKDP(n,m,s); s->sugar = d->td;
2001: *rp = s;
2002: }
2003:
1.20 noro 2004: void dltod(DL d,int n,DP *rp)
1.7 noro 2005: {
2006: MP m;
2007: DP s;
2008:
2009: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
2010: MKDP(n,m,s); s->sugar = d->td;
2011: *rp = s;
1.5 noro 2012: }
2013:
1.20 noro 2014: void dp_hm(DP p,DP *rp)
1.5 noro 2015: {
2016: MP m,mr;
2017:
2018: if ( !p )
2019: *rp = 0;
2020: else {
2021: m = BDY(p);
2022: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
2023: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
2024: }
2025: }
2026:
1.35 noro 2027: void dp_ht(DP p,DP *rp)
2028: {
2029: MP m,mr;
2030:
2031: if ( !p )
2032: *rp = 0;
2033: else {
2034: m = BDY(p);
2035: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
2036: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
2037: }
2038: }
2039:
1.20 noro 2040: void dp_rest(DP p,DP *rp)
1.5 noro 2041: {
2042: MP m;
2043:
2044: m = BDY(p);
2045: if ( !NEXT(m) )
2046: *rp = 0;
2047: else {
2048: MKDP(p->nv,NEXT(m),*rp);
2049: if ( *rp )
2050: (*rp)->sugar = p->sugar;
2051: }
2052: }
2053:
1.20 noro 2054: DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl)
1.5 noro 2055: {
1.21 noro 2056: register int i, *d1, *d2, *d, td;
1.5 noro 2057:
2058: if ( !dl ) NEWDL(dl,nv);
2059: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1.21 noro 2060: for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) {
2061: *d = *d1 > *d2 ? *d1 : *d2;
2062: td += MUL_WEIGHT(*d,i);
2063: }
1.5 noro 2064: dl->td = td;
2065: return dl;
2066: }
2067:
1.20 noro 2068: int dl_equal(int nv,DL dl1,DL dl2)
1.5 noro 2069: {
2070: register int *d1, *d2, n;
2071:
2072: if ( dl1->td != dl2->td ) return 0;
2073: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
2074: if ( *d1 != *d2 ) return 0;
2075: return 1;
2076: }
2077:
1.20 noro 2078: int dp_nt(DP p)
1.5 noro 2079: {
2080: int i;
2081: MP m;
2082:
2083: if ( !p )
2084: return 0;
2085: else {
2086: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
2087: return i;
2088: }
2089: }
2090:
1.20 noro 2091: int dp_homogeneous(DP p)
1.15 noro 2092: {
2093: MP m;
2094: int d;
2095:
2096: if ( !p )
2097: return 1;
2098: else {
2099: m = BDY(p);
2100: d = m->dl->td;
2101: m = NEXT(m);
2102: for ( ; m; m = NEXT(m) ) {
2103: if ( m->dl->td != d )
2104: return 0;
2105: }
2106: return 1;
2107: }
1.16 noro 2108: }
2109:
1.20 noro 2110: void _print_mp(int nv,MP m)
1.16 noro 2111: {
2112: int i;
2113:
1.17 noro 2114: if ( !m )
1.16 noro 2115: return;
2116: for ( ; m; m = NEXT(m) ) {
2117: fprintf(stderr,"%d<",ITOS(C(m)));
2118: for ( i = 0; i < nv; i++ ) {
2119: fprintf(stderr,"%d",m->dl->d[i]);
2120: if ( i != nv-1 )
2121: fprintf(stderr," ");
2122: }
2123: fprintf(stderr,">",C(m));
2124: }
2125: fprintf(stderr,"\n");
1.15 noro 2126: }
1.26 noro 2127:
2128: static int cmp_mp_nvar;
2129:
2130: int comp_mp(MP *a,MP *b)
2131: {
2132: return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl);
2133: }
2134:
2135: void dp_sort(DP p,DP *rp)
2136: {
2137: MP t,mp,mp0;
2138: int i,n;
2139: DP r;
2140: MP *w;
2141:
2142: if ( !p ) {
2143: *rp = 0;
2144: return;
2145: }
2146: for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ );
2147: w = (MP *)ALLOCA(n*sizeof(MP));
2148: for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ )
2149: w[i] = t;
2150: cmp_mp_nvar = NV(p);
2151: qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp);
2152: mp0 = 0;
2153: for ( i = n-1; i >= 0; i-- ) {
2154: NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]);
2155: NEXT(mp) = mp0; mp0 = mp;
2156: }
2157: MKDP(p->nv,mp0,r);
2158: r->sugar = p->sugar;
2159: *rp = r;
2160: }
2161:
1.32 noro 2162: DP extract_initial_term_from_dp(DP p,int *weight,int n);
2163: LIST extract_initial_term(LIST f,int *weight,int n);
2164:
2165: DP extract_initial_term_from_dp(DP p,int *weight,int n)
2166: {
1.34 noro 2167: int w,t,i,top;
1.32 noro 2168: MP m,r0,r;
2169: DP dp;
2170:
2171: if ( !p ) return 0;
1.34 noro 2172: top = 1;
1.32 noro 2173: for ( m = BDY(p); m; m = NEXT(m) ) {
2174: for ( i = 0, t = 0; i < n; i++ )
2175: t += weight[i]*m->dl->d[i];
1.34 noro 2176: if ( top || t > w ) {
1.32 noro 2177: r0 = 0;
2178: w = t;
1.34 noro 2179: top = 0;
1.32 noro 2180: }
2181: if ( t == w ) {
2182: NEXTMP(r0,r);
2183: r->dl = m->dl;
2184: r->c = m->c;
2185: }
2186: }
2187: NEXT(r) = 0;
2188: MKDP(p->nv,r0,dp);
2189: return dp;
2190: }
2191:
2192: LIST extract_initial_term(LIST f,int *weight,int n)
2193: {
2194: NODE nd,r0,r;
2195: Obj p;
2196: LIST l;
2197:
2198: nd = BDY(f);
2199: for ( r0 = 0; nd; nd = NEXT(nd) ) {
2200: NEXTNODE(r0,r);
2201: p = (Obj)BDY(nd);
2202: BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n);
2203: }
2204: if ( r0 ) NEXT(r) = 0;
2205: MKLIST(l,r0);
2206: return l;
2207: }
2208:
2209: LIST dp_initial_term(LIST f,struct order_spec *ord)
2210: {
2211: int n,l,i;
2212: struct weight_or_block *worb;
2213: int *weight;
2214:
2215: switch ( ord->id ) {
2216: case 2: /* matrix order */
2217: /* extract the first row */
2218: n = ord->nv;
2219: weight = ord->ord.matrix.matrix[0];
2220: return extract_initial_term(f,weight,n);
2221: case 3: /* composite order */
2222: /* the first w_or_b */
2223: worb = ord->ord.composite.w_or_b;
2224: switch ( worb->type ) {
2225: case IS_DENSE_WEIGHT:
2226: n = worb->length;
2227: weight = worb->body.dense_weight;
2228: return extract_initial_term(f,weight,n);
2229: case IS_SPARSE_WEIGHT:
2230: n = ord->nv;
2231: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 2232: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 2233: l = worb->length;
2234: for ( i = 0; i < l; i++ )
2235: weight[worb->body.sparse_weight[i].pos]
2236: = worb->body.sparse_weight[i].value;
2237: return extract_initial_term(f,weight,n);
2238: default:
2239: error("dp_initial_term : unsupported order");
2240: }
2241: default:
2242: error("dp_initial_term : unsupported order");
2243: }
2244: }
2245:
2246: int highest_order_dp(DP p,int *weight,int n);
2247: LIST highest_order(LIST f,int *weight,int n);
2248:
2249: int highest_order_dp(DP p,int *weight,int n)
2250: {
1.34 noro 2251: int w,t,i,top;
1.32 noro 2252: MP m;
2253:
2254: if ( !p ) return -1;
1.34 noro 2255: top = 1;
1.32 noro 2256: for ( m = BDY(p); m; m = NEXT(m) ) {
2257: for ( i = 0, t = 0; i < n; i++ )
2258: t += weight[i]*m->dl->d[i];
1.34 noro 2259: if ( top || t > w ) {
1.32 noro 2260: w = t;
1.34 noro 2261: top = 0;
2262: }
1.32 noro 2263: }
2264: return w;
2265: }
2266:
2267: LIST highest_order(LIST f,int *weight,int n)
2268: {
2269: int h;
2270: NODE nd,r0,r;
2271: Obj p;
2272: LIST l;
2273: Q q;
2274:
2275: nd = BDY(f);
2276: for ( r0 = 0; nd; nd = NEXT(nd) ) {
2277: NEXTNODE(r0,r);
2278: p = (Obj)BDY(nd);
2279: h = highest_order_dp((DP)p,weight,n);
2280: STOQ(h,q);
2281: BDY(r) = (pointer)q;
2282: }
2283: if ( r0 ) NEXT(r) = 0;
2284: MKLIST(l,r0);
2285: return l;
2286: }
2287:
2288: LIST dp_order(LIST f,struct order_spec *ord)
2289: {
2290: int n,l,i;
2291: struct weight_or_block *worb;
2292: int *weight;
2293:
2294: switch ( ord->id ) {
2295: case 2: /* matrix order */
2296: /* extract the first row */
2297: n = ord->nv;
2298: weight = ord->ord.matrix.matrix[0];
2299: return highest_order(f,weight,n);
2300: case 3: /* composite order */
2301: /* the first w_or_b */
2302: worb = ord->ord.composite.w_or_b;
2303: switch ( worb->type ) {
2304: case IS_DENSE_WEIGHT:
2305: n = worb->length;
2306: weight = worb->body.dense_weight;
2307: return highest_order(f,weight,n);
2308: case IS_SPARSE_WEIGHT:
2309: n = ord->nv;
2310: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 2311: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 2312: l = worb->length;
2313: for ( i = 0; i < l; i++ )
2314: weight[worb->body.sparse_weight[i].pos]
2315: = worb->body.sparse_weight[i].value;
2316: return highest_order(f,weight,n);
2317: default:
2318: error("dp_initial_term : unsupported order");
2319: }
2320: default:
2321: error("dp_initial_term : unsupported order");
1.35 noro 2322: }
2323: }
2324:
2325: int dpv_ht(DPV p,DP *h)
2326: {
2327: int len,max,maxi,i,t;
2328: DP *e;
2329: MP m,mr;
2330:
2331: len = p->len;
2332: e = p->body;
2333: max = -1;
2334: maxi = -1;
2335: for ( i = 0; i < len; i++ )
2336: if ( e[i] && (t = BDY(e[i])->dl->td) > max ) {
2337: max = t;
2338: maxi = i;
2339: }
2340: if ( max < 0 ) {
2341: *h = 0;
2342: return -1;
2343: } else {
2344: m = BDY(e[maxi]);
2345: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
2346: MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td; /* XXX */
2347: return maxi;
1.32 noro 2348: }
2349: }