![[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.45 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.44 2007/09/15 10:17:08 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.44 noro 863: void dp_red_marked_mod(DP p0,DP p1,DP p2,DP hp2,int mod,DP *head,DP *rest,P *dnp)
864: {
865: int i,n;
866: DL d1,d2,d;
867: MP m;
868: DP t,s,r,h;
869: P c1,c2,g,u;
870:
871: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl;
872: NEWDL(d,n); d->td = d1->td - d2->td;
873: for ( i = 0; i < n; i++ )
874: d->d[i] = d1->d[i]-d2->d[i];
875: c1 = (P)BDY(p1)->c; c2 = (P)BDY(hp2)->c;
876: gcdprsmp(CO,mod,c1,c2,&g);
877: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
878: if ( NUM(c2) ) {
879: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
880: }
881: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
882: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
883: if ( NUM(c2) ) {
884: addmd(CO,mod,p1,t,&r); h = p0;
885: } else {
886: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
887: }
888: *head = h; *rest = r; *dnp = c2;
889: }
890:
1.13 noro 891: /* m-reduction over a field */
892:
1.20 noro 893: void dp_red_f(DP p1,DP p2,DP *rest)
1.13 noro 894: {
895: int i,n;
896: DL d1,d2,d;
897: MP m;
1.20 noro 898: DP t,s;
1.13 noro 899: Obj a,b;
900:
901: n = p1->nv;
902: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
903:
904: NEWDL(d,n); d->td = d1->td - d2->td;
905: for ( i = 0; i < n; i++ )
906: d->d[i] = d1->d[i]-d2->d[i];
907:
908: NEWMP(m); m->dl = d;
909: divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
910: C(m) = (P)b;
911: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
912:
913: muld(CO,s,p2,&t); addd(CO,p1,t,rest);
914: }
915:
1.20 noro 916: void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp)
1.7 noro 917: {
918: int i,n;
919: DL d1,d2,d;
920: MP m;
921: DP t,s,r,h;
922: P c1,c2,g,u;
923:
924: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
925: NEWDL(d,n); d->td = d1->td - d2->td;
926: for ( i = 0; i < n; i++ )
927: d->d[i] = d1->d[i]-d2->d[i];
928: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
929: gcdprsmp(CO,mod,c1,c2,&g);
930: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
931: if ( NUM(c2) ) {
932: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
933: }
934: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1.11 noro 935: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7 noro 936: if ( NUM(c2) ) {
937: addmd(CO,mod,p1,t,&r); h = p0;
938: } else {
939: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
940: }
941: *head = h; *rest = r; *dnp = c2;
942: }
943:
1.10 noro 944: struct oEGT eg_red_mod;
945:
1.20 noro 946: void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp)
1.7 noro 947: {
948: int i,n;
949: DL d1,d2,d;
950: MP m;
951: DP t,s;
1.16 noro 952: int c,c1,c2;
953: extern int do_weyl;
1.7 noro 954:
955: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
956: _NEWDL(d,n); d->td = d1->td - d2->td;
957: for ( i = 0; i < n; i++ )
958: d->d[i] = d1->d[i]-d2->d[i];
1.16 noro 959: c = invm(ITOS(BDY(p2)->c),mod);
960: c2 = ITOS(BDY(p1)->c);
961: DMAR(c,c2,0,mod,c1);
1.7 noro 962: _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
1.16 noro 963: #if 0
1.7 noro 964: _MKDP(n,m,s); s->sugar = d->td;
965: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 966: #else
967: if ( do_weyl ) {
1.19 noro 968: _MKDP(n,m,s); s->sugar = d->td;
969: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 970: } else {
971: _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
972: }
973: #endif
1.10 noro 974: /* get_eg(&t0); */
1.7 noro 975: _addmd_destructive(mod,p1,t,rp);
1.10 noro 976: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 977: }
978:
979: /*
980: * normal form computation
981: *
982: */
1.5 noro 983:
1.20 noro 984: void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp)
1.5 noro 985: {
986: DP u,p,d,s,t,dmy;
987: NODE l;
988: MP m,mr;
989: int i,n;
990: int *wb;
991: int sugar,psugar;
992: P dn,tdn,tdn1;
993:
994: dn = (P)ONE;
995: if ( !g ) {
996: *rp = 0; *dnp = dn; return;
997: }
998: for ( n = 0, l = b; l; l = NEXT(l), n++ );
999: wb = (int *)ALLOCA(n*sizeof(int));
1000: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1001: wb[i] = QTOS((Q)BDY(l));
1002: sugar = g->sugar;
1003: for ( d = 0; g; ) {
1004: for ( u = 0, i = 0; i < n; i++ ) {
1005: if ( dp_redble(g,p = ps[wb[i]]) ) {
1006: dp_red(d,g,p,&t,&u,&tdn,&dmy);
1007: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1008: sugar = MAX(sugar,psugar);
1009: if ( !u ) {
1010: if ( d )
1011: d->sugar = sugar;
1012: *rp = d; *dnp = dn; return;
1013: } else {
1014: d = t;
1015: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
1016: }
1017: break;
1018: }
1019: }
1020: if ( u )
1021: g = u;
1022: else if ( !full ) {
1023: if ( g ) {
1024: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1025: }
1026: *rp = g; *dnp = dn; return;
1027: } else {
1028: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1029: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1030: addd(CO,d,t,&s); d = s;
1031: dp_rest(g,&t); g = t;
1032: }
1033: }
1034: if ( d )
1035: d->sugar = sugar;
1036: *rp = d; *dnp = dn;
1037: }
1038:
1.43 noro 1039: void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Q *contp)
1040: {
1041: struct oVECT v;
1042: int i,n1,n2,n;
1043: MP m,m0,t;
1044: Q *w;
1045: Q h;
1046:
1047: if ( p1 ) {
1048: for ( i = 0, m = BDY(p1); m; m = NEXT(m), i++ );
1049: n1 = i;
1050: } else
1051: n1 = 0;
1052: if ( p2 ) {
1053: for ( i = 0, m = BDY(p2); m; m = NEXT(m), i++ );
1054: n2 = i;
1055: } else
1056: n2 = 0;
1057: n = n1+n2;
1058: if ( !n ) {
1059: *r1p = 0; *r2p = 0; *contp = ONE; return;
1060: }
1061: w = (Q *)ALLOCA(n*sizeof(Q));
1062: v.len = n;
1063: v.body = (pointer *)w;
1064: i = 0;
1065: if ( p1 )
1066: for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) w[i] = (Q)m->c;
1067: if ( p2 )
1068: for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = (Q)m->c;
1069: h = w[0]; removecont_array((P *)w,n,1); divq(h,w[0],contp);
1070: i = 0;
1071: if ( p1 ) {
1072: for ( m0 = 0, t = BDY(p1); i < n1; i++, t = NEXT(t) ) {
1073: NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl;
1074: }
1075: NEXT(m) = 0;
1076: MKDP(p1->nv,m0,*r1p); (*r1p)->sugar = p1->sugar;
1077: } else
1078: *r1p = 0;
1079: if ( p2 ) {
1080: for ( m0 = 0, t = BDY(p2); i < n; i++, t = NEXT(t) ) {
1081: NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl;
1082: }
1083: NEXT(m) = 0;
1084: MKDP(p2->nv,m0,*r2p); (*r2p)->sugar = p2->sugar;
1085: } else
1086: *r2p = 0;
1087: }
1088:
1.41 noro 1089: /* true nf by a marked GB */
1090:
1.43 noro 1091: void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP *rp,P *nmp,P *dnp)
1.41 noro 1092: {
1093: DP u,p,d,s,t,dmy,hp;
1094: NODE l;
1095: MP m,mr;
1.43 noro 1096: int i,n,hmag;
1.41 noro 1097: int *wb;
1.43 noro 1098: int sugar,psugar,multiple;
1099: P nm,tnm1,dn,tdn,tdn1;
1100: Q cont;
1.41 noro 1101:
1.43 noro 1102: multiple = 0;
1103: hmag = multiple*HMAG(g);
1104: nm = (P)ONE;
1.41 noro 1105: dn = (P)ONE;
1106: if ( !g ) {
1107: *rp = 0; *dnp = dn; return;
1108: }
1109: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1110: wb = (int *)ALLOCA(n*sizeof(int));
1111: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1112: wb[i] = QTOS((Q)BDY(l));
1113: sugar = g->sugar;
1114: for ( d = 0; g; ) {
1115: for ( u = 0, i = 0; i < n; i++ ) {
1116: if ( dp_redble(g,hp = hps[wb[i]]) ) {
1117: p = ps[wb[i]];
1118: dp_red_marked(d,g,p,hp,&t,&u,&tdn,&dmy);
1119: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1120: sugar = MAX(sugar,psugar);
1121: if ( !u ) {
1.43 noro 1122: goto last;
1.41 noro 1123: } else {
1124: d = t;
1125: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
1126: }
1127: break;
1128: }
1129: }
1.43 noro 1130: if ( u ) {
1.41 noro 1131: g = u;
1.43 noro 1132: if ( multiple && ((d && HMAG(d)>hmag) || (HMAG(g)>hmag)) ) {
1133: dp_removecont2(d,g,&t,&u,&cont); d = t; g = u;
1134: mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
1135: if ( d )
1136: hmag = multiple*HMAG(d);
1137: else
1138: hmag = multiple*HMAG(g);
1139: }
1140: } else {
1.41 noro 1141: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1142: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1143: addd(CO,d,t,&s); d = s;
1144: dp_rest(g,&t); g = t;
1145: }
1146: }
1.43 noro 1147: last:
1148: if ( d ) {
1149: dp_removecont2(d,0,&t,&u,&cont); d = t;
1150: mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
1.41 noro 1151: d->sugar = sugar;
1.43 noro 1152: }
1153: *rp = d; *nmp = nm; *dnp = dn;
1.41 noro 1154: }
1155:
1.44 noro 1156: void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
1157: {
1158: DP hp,u,p,d,s,t;
1159: NODE l;
1160: MP m,mr;
1161: int i,n;
1162: int *wb;
1163: int sugar,psugar;
1164: P dn,tdn,tdn1;
1165:
1166: dn = (P)ONEM;
1167: if ( !g ) {
1168: *rp = 0; *dnp = dn; return;
1169: }
1170: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1171: wb = (int *)ALLOCA(n*sizeof(int));
1172: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1173: wb[i] = QTOS((Q)BDY(l));
1174: sugar = g->sugar;
1175: for ( d = 0; g; ) {
1176: for ( u = 0, i = 0; i < n; i++ ) {
1177: if ( dp_redble(g,hp = hps[wb[i]]) ) {
1178: p = ps[wb[i]];
1179: dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn);
1180: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1181: sugar = MAX(sugar,psugar);
1182: if ( !u ) {
1183: if ( d )
1184: d->sugar = sugar;
1185: *rp = d; *dnp = dn; return;
1186: } else {
1187: d = t;
1188: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1189: }
1190: break;
1191: }
1192: }
1193: if ( u )
1194: g = u;
1195: else {
1196: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1197: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1198: addmd(CO,mod,d,t,&s); d = s;
1199: dp_rest(g,&t); g = t;
1200: }
1201: }
1202: if ( d )
1203: d->sugar = sugar;
1204: *rp = d; *dnp = dn;
1205: }
1206:
1.13 noro 1207: /* nf computation over Z */
1208:
1.20 noro 1209: void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
1.5 noro 1210: {
1211: DP u,p,d,s,t,dmy1;
1212: P dmy;
1213: NODE l;
1214: MP m,mr;
1215: int i,n;
1216: int *wb;
1217: int hmag;
1218: int sugar,psugar;
1219:
1220: if ( !g ) {
1221: *rp = 0; return;
1222: }
1223: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1224: wb = (int *)ALLOCA(n*sizeof(int));
1225: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1226: wb[i] = QTOS((Q)BDY(l));
1.12 noro 1227:
1.13 noro 1228: hmag = multiple*HMAG(g);
1.5 noro 1229: sugar = g->sugar;
1.12 noro 1230:
1.5 noro 1231: for ( d = 0; g; ) {
1232: for ( u = 0, i = 0; i < n; i++ ) {
1233: if ( dp_redble(g,p = ps[wb[i]]) ) {
1234: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
1235: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1236: sugar = MAX(sugar,psugar);
1237: if ( !u ) {
1238: if ( d )
1239: d->sugar = sugar;
1240: *rp = d; return;
1241: }
1242: d = t;
1243: break;
1244: }
1245: }
1246: if ( u ) {
1247: g = u;
1248: if ( d ) {
1.13 noro 1249: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 1250: dp_ptozp2(d,g,&t,&u); d = t; g = u;
1251: hmag = multiple*HMAG(d);
1252: }
1253: } else {
1.13 noro 1254: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 1255: dp_ptozp(g,&t); g = t;
1256: hmag = multiple*HMAG(g);
1257: }
1258: }
1259: }
1260: else if ( !full ) {
1261: if ( g ) {
1262: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1263: }
1264: *rp = g; return;
1265: } else {
1266: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1267: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1268: addd(CO,d,t,&s); d = s;
1269: dp_rest(g,&t); g = t;
1270:
1271: }
1272: }
1273: if ( d )
1274: d->sugar = sugar;
1275: *rp = d;
1276: }
1277:
1.13 noro 1278: /* nf computation over a field */
1279:
1.20 noro 1280: void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp)
1.13 noro 1281: {
1282: DP u,p,d,s,t;
1283: NODE l;
1284: MP m,mr;
1285: int i,n;
1286: int *wb;
1287: int sugar,psugar;
1288:
1289: if ( !g ) {
1290: *rp = 0; return;
1291: }
1292: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1293: wb = (int *)ALLOCA(n*sizeof(int));
1294: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1295: wb[i] = QTOS((Q)BDY(l));
1296:
1297: sugar = g->sugar;
1298: for ( d = 0; g; ) {
1299: for ( u = 0, i = 0; i < n; i++ ) {
1300: if ( dp_redble(g,p = ps[wb[i]]) ) {
1301: dp_red_f(g,p,&u);
1302: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1303: sugar = MAX(sugar,psugar);
1304: if ( !u ) {
1305: if ( d )
1306: d->sugar = sugar;
1307: *rp = d; return;
1308: }
1309: break;
1310: }
1311: }
1312: if ( u )
1313: g = u;
1314: else if ( !full ) {
1315: if ( g ) {
1316: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1317: }
1318: *rp = g; return;
1319: } else {
1320: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1321: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1322: addd(CO,d,t,&s); d = s;
1323: dp_rest(g,&t); g = t;
1324: }
1325: }
1326: if ( d )
1327: d->sugar = sugar;
1328: *rp = d;
1329: }
1330:
1331: /* nf computation over GF(mod) (only for internal use) */
1332:
1.20 noro 1333: void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1334: {
1335: DP u,p,d,s,t;
1336: P dmy;
1337: NODE l;
1338: MP m,mr;
1339: int sugar,psugar;
1340:
1341: if ( !g ) {
1342: *rp = 0; return;
1343: }
1344: sugar = g->sugar;
1345: for ( d = 0; g; ) {
1346: for ( u = 0, l = b; l; l = NEXT(l) ) {
1347: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1348: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
1349: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1350: sugar = MAX(sugar,psugar);
1351: if ( !u ) {
1352: if ( d )
1353: d->sugar = sugar;
1354: *rp = d; return;
1355: }
1356: d = t;
1357: break;
1358: }
1359: }
1360: if ( u )
1361: g = u;
1362: else if ( !full ) {
1363: if ( g ) {
1364: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1365: }
1366: *rp = g; return;
1367: } else {
1368: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1369: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1370: addmd(CO,mod,d,t,&s); d = s;
1371: dp_rest(g,&t); g = t;
1372: }
1373: }
1374: if ( d )
1375: d->sugar = sugar;
1376: *rp = d;
1377: }
1378:
1.20 noro 1379: void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp)
1.5 noro 1380: {
1381: DP u,p,d,s,t;
1382: NODE l;
1383: MP m,mr;
1384: int i,n;
1385: int *wb;
1386: int sugar,psugar;
1387: P dn,tdn,tdn1;
1388:
1389: dn = (P)ONEM;
1390: if ( !g ) {
1391: *rp = 0; *dnp = dn; return;
1392: }
1393: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1394: wb = (int *)ALLOCA(n*sizeof(int));
1395: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1396: wb[i] = QTOS((Q)BDY(l));
1397: sugar = g->sugar;
1398: for ( d = 0; g; ) {
1399: for ( u = 0, i = 0; i < n; i++ ) {
1400: if ( dp_redble(g,p = ps[wb[i]]) ) {
1401: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
1402: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1403: sugar = MAX(sugar,psugar);
1404: if ( !u ) {
1405: if ( d )
1406: d->sugar = sugar;
1407: *rp = d; *dnp = dn; return;
1408: } else {
1409: d = t;
1410: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1411: }
1412: break;
1413: }
1414: }
1415: if ( u )
1416: g = u;
1417: else if ( !full ) {
1418: if ( g ) {
1419: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1420: }
1421: *rp = g; *dnp = dn; return;
1422: } else {
1423: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1424: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1425: addmd(CO,mod,d,t,&s); d = s;
1426: dp_rest(g,&t); g = t;
1427: }
1428: }
1429: if ( d )
1430: d->sugar = sugar;
1431: *rp = d; *dnp = dn;
1432: }
1433:
1.20 noro 1434: void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1435: {
1.20 noro 1436: DP u,p,d;
1.7 noro 1437: NODE l;
1.20 noro 1438: MP m,mrd;
1439: int sugar,psugar,n,h_reducible;
1.5 noro 1440:
1.7 noro 1441: if ( !g ) {
1442: *rp = 0; return;
1.5 noro 1443: }
1.7 noro 1444: sugar = g->sugar;
1445: n = g->nv;
1446: for ( d = 0; g; ) {
1447: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
1448: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1449: h_reducible = 1;
1450: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1451: _dp_red_mod_destructive(g,p,mod,&u); g = u;
1452: sugar = MAX(sugar,psugar);
1453: if ( !g ) {
1454: if ( d )
1455: d->sugar = sugar;
1456: _dptodp(d,rp); _free_dp(d); return;
1457: }
1458: break;
1459: }
1460: }
1461: if ( !h_reducible ) {
1462: /* head term is not reducible */
1463: if ( !full ) {
1464: if ( g )
1465: g->sugar = sugar;
1466: _dptodp(g,rp); _free_dp(g); return;
1467: } else {
1468: m = BDY(g);
1469: if ( NEXT(m) ) {
1470: BDY(g) = NEXT(m); NEXT(m) = 0;
1471: } else {
1472: _FREEDP(g); g = 0;
1473: }
1474: if ( d ) {
1475: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1476: NEXT(mrd) = m;
1477: } else {
1478: _MKDP(n,m,d);
1479: }
1480: }
1481: }
1.5 noro 1482: }
1.7 noro 1483: if ( d )
1484: d->sugar = sugar;
1485: _dptodp(d,rp); _free_dp(d);
1.5 noro 1486: }
1.13 noro 1487:
1488: /* reduction by linear base over a field */
1489:
1.20 noro 1490: void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p)
1.13 noro 1491: {
1492: DP r1,r2,b1,b2,t,s;
1493: Obj c,c1,c2;
1494: NODE l,b;
1495: int n;
1496:
1497: if ( !p1 ) {
1498: *r1p = p1; *r2p = p2; return;
1499: }
1500: n = p1->nv;
1501: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1502: if ( !r1 ) {
1503: *r1p = r1; *r2p = r2; return;
1504: }
1505: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1506: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1507: b2 = (DP)BDY(NEXT(b));
1508: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
1509: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
1510: muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
1511: muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
1512: }
1513: }
1514: *r1p = r1; *r2p = r2;
1515: }
1516:
1517: /* reduction by linear base over GF(mod) */
1.5 noro 1518:
1.20 noro 1519: void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p)
1.5 noro 1520: {
1.7 noro 1521: DP r1,r2,b1,b2,t,s;
1522: P c;
1523: MQ c1,c2;
1524: NODE l,b;
1525: int n;
1526:
1527: if ( !p1 ) {
1528: *r1p = p1; *r2p = p2; return;
1529: }
1530: n = p1->nv;
1531: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1532: if ( !r1 ) {
1533: *r1p = r1; *r2p = r2; return;
1534: }
1535: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1536: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1537: b2 = (DP)BDY(NEXT(b));
1538: invmq(mod,(MQ)BDY(b1)->c,&c1);
1539: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
1540: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
1541: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
1542: }
1543: }
1544: *r1p = r1; *r2p = r2;
1.5 noro 1545: }
1546:
1.20 noro 1547: void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp)
1.5 noro 1548: {
1.7 noro 1549: DP s,t,u;
1550: MP m;
1551: DL h;
1552: int i,n;
1553:
1554: if ( !p ) {
1555: *rp = p; return;
1556: }
1557: n = p->nv;
1558: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1559: h = m->dl;
1560: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1561: i++;
1562: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1563: addmd(CO,mod,s,t,&u); s = u;
1.24 noro 1564: }
1565: *rp = s;
1566: }
1567:
1568: void dp_nf_tab_f(DP p,LIST *tab,DP *rp)
1569: {
1570: DP s,t,u;
1571: MP m;
1572: DL h;
1573: int i,n;
1574:
1575: if ( !p ) {
1576: *rp = p; return;
1577: }
1578: n = p->nv;
1579: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1580: h = m->dl;
1581: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1582: i++;
1583: muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1584: addd(CO,s,t,&u); s = u;
1.7 noro 1585: }
1586: *rp = s;
1.5 noro 1587: }
1588:
1.7 noro 1589: /*
1590: * setting flags
1.30 noro 1591: * call create_order_spec with vl=0 to set old type order.
1.7 noro 1592: *
1593: */
1594:
1.27 noro 1595: int create_order_spec(VL vl,Obj obj,struct order_spec **specp)
1.5 noro 1596: {
1.37 noro 1597: int i,j,n,s,row,col,ret;
1.27 noro 1598: struct order_spec *spec;
1.7 noro 1599: struct order_pair *l;
1600: NODE node,t,tn;
1601: MAT m;
1602: pointer **b;
1603: int **w;
1.5 noro 1604:
1.37 noro 1605: if ( vl && obj && OID(obj) == O_LIST ) {
1606: ret = create_composite_order_spec(vl,(LIST)obj,specp);
1607: if ( show_orderspec )
1608: print_composite_order_spec(*specp);
1609: return ret;
1610: }
1.27 noro 1611:
1612: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 1613: if ( !obj || NUM(obj) ) {
1614: spec->id = 0; spec->obj = obj;
1615: spec->ord.simple = QTOS((Q)obj);
1616: return 1;
1617: } else if ( OID(obj) == O_LIST ) {
1618: node = BDY((LIST)obj);
1619: for ( n = 0, t = node; t; t = NEXT(t), n++ );
1620: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1621: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
1622: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
1623: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
1624: s += l[i].length;
1625: }
1626: spec->id = 1; spec->obj = obj;
1627: spec->ord.block.order_pair = l;
1628: spec->ord.block.length = n; spec->nv = s;
1629: return 1;
1630: } else if ( OID(obj) == O_MAT ) {
1631: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
1632: w = almat(row,col);
1633: for ( i = 0; i < row; i++ )
1634: for ( j = 0; j < col; j++ )
1635: w[i][j] = QTOS((Q)b[i][j]);
1636: spec->id = 2; spec->obj = obj;
1637: spec->nv = col; spec->ord.matrix.row = row;
1638: spec->ord.matrix.matrix = w;
1639: return 1;
1640: } else
1.5 noro 1641: return 0;
1642: }
1643:
1.28 noro 1644: void print_composite_order_spec(struct order_spec *spec)
1645: {
1646: int nv,n,len,i,j,k,start;
1647: struct weight_or_block *worb;
1648:
1649: nv = spec->nv;
1650: n = spec->ord.composite.length;
1651: worb = spec->ord.composite.w_or_b;
1652: for ( i = 0; i < n; i++, worb++ ) {
1653: len = worb->length;
1654: printf("[ ");
1655: switch ( worb->type ) {
1656: case IS_DENSE_WEIGHT:
1657: for ( j = 0; j < len; j++ )
1658: printf("%d ",worb->body.dense_weight[j]);
1659: for ( ; j < nv; j++ )
1660: printf("0 ");
1661: break;
1662: case IS_SPARSE_WEIGHT:
1663: for ( j = 0, k = 0; j < nv; j++ )
1664: if ( j == worb->body.sparse_weight[k].pos )
1665: printf("%d ",worb->body.sparse_weight[k++].value);
1666: else
1667: printf("0 ");
1668: break;
1669: case IS_BLOCK:
1670: start = worb->body.block.start;
1671: for ( j = 0; j < start; j++ ) printf("0 ");
1672: switch ( worb->body.block.order ) {
1673: case 0:
1674: for ( k = 0; k < len; k++, j++ ) printf("R ");
1675: break;
1676: case 1:
1677: for ( k = 0; k < len; k++, j++ ) printf("G ");
1678: break;
1679: case 2:
1680: for ( k = 0; k < len; k++, j++ ) printf("L ");
1681: break;
1682: }
1683: for ( ; j < nv; j++ ) printf("0 ");
1684: break;
1685: }
1686: printf("]\n");
1687: }
1.38 noro 1688: }
1689:
1690: struct order_spec *append_block(struct order_spec *spec,
1691: int nv,int nalg,int ord)
1692: {
1693: MAT m,mat;
1694: int i,j,row,col,n;
1695: Q **b,**wp;
1696: int **w;
1697: NODE t,s,s0;
1698: struct order_pair *l,*l0;
1699: int n0,nv0;
1700: LIST list0,list1,list;
1701: Q oq,nq;
1702: struct order_spec *r;
1703:
1704: r = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1705: switch ( spec->id ) {
1706: case 0:
1707: STOQ(spec->ord.simple,oq); STOQ(nv,nq);
1708: t = mknode(2,oq,nq); MKLIST(list0,t);
1709: STOQ(ord,oq); STOQ(nalg,nq);
1710: t = mknode(2,oq,nq); MKLIST(list1,t);
1711: t = mknode(2,list0,list1); MKLIST(list,t);
1712: l = (struct order_pair *)MALLOC_ATOMIC(2*sizeof(struct order_pair));
1713: l[0].order = spec->ord.simple; l[0].length = nv;
1714: l[1].order = ord; l[1].length = nalg;
1715: r->id = 1; r->obj = (Obj)list;
1716: r->ord.block.order_pair = l;
1717: r->ord.block.length = 2;
1718: r->nv = nv+nalg;
1719: break;
1720: case 1:
1721: if ( spec->nv != nv )
1722: error("append_block : number of variables mismatch");
1723: l0 = spec->ord.block.order_pair;
1724: n0 = spec->ord.block.length;
1725: nv0 = spec->nv;
1726: list0 = (LIST)spec->obj;
1727: n = n0+1;
1728: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1729: for ( i = 0; i < n0; i++ )
1730: l[i] = l0[i];
1731: l[i].order = ord; l[i].length = nalg;
1732: for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) {
1733: NEXTNODE(s0,s); BDY(s) = BDY(t);
1734: }
1735: STOQ(ord,oq); STOQ(nalg,nq);
1736: t = mknode(2,oq,nq); MKLIST(list,t);
1737: NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0;
1738: MKLIST(list,s0);
1739: r->id = 1; r->obj = (Obj)list;
1740: r->ord.block.order_pair = l;
1741: r->ord.block.length = n;
1742: r->nv = nv+nalg;
1743: break;
1744: case 2:
1745: if ( spec->nv != nv )
1746: error("append_block : number of variables mismatch");
1747: m = (MAT)spec->obj;
1748: row = m->row; col = m->col; b = (Q **)BDY(m);
1749: w = almat(row+nalg,col+nalg);
1750: MKMAT(mat,row+nalg,col+nalg); wp = (Q **)BDY(mat);
1751: for ( i = 0; i < row; i++ )
1752: for ( j = 0; j < col; j++ ) {
1753: w[i][j] = QTOS(b[i][j]);
1754: wp[i][j] = b[i][j];
1755: }
1756: for ( i = 0; i < nalg; i++ ) {
1757: w[i+row][i+col] = 1;
1758: wp[i+row][i+col] = ONE;
1759: }
1760: r->id = 2; r->obj = (Obj)mat;
1761: r->nv = col+nalg; r->ord.matrix.row = row+nalg;
1762: r->ord.matrix.matrix = w;
1763: break;
1764: case 3:
1765: default:
1766: /* XXX */
1767: error("append_block : not implemented yet");
1768: }
1769: return r;
1.28 noro 1770: }
1771:
1.37 noro 1772: int comp_sw(struct sparse_weight *a, struct sparse_weight *b)
1773: {
1774: if ( a->pos > b->pos ) return 1;
1775: else if ( a->pos < b->pos ) return -1;
1776: else return 0;
1777: }
1778:
1.27 noro 1779: /* order = [w_or_b, w_or_b, ... ] */
1780: /* w_or_b = w or b */
1781: /* w = [1,2,...] or [x,1,y,2,...] */
1782: /* b = [@lex,x,y,...,z] etc */
1783:
1784: int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp)
1785: {
1786: NODE wb,t,p;
1787: struct order_spec *spec;
1788: VL tvl;
1.29 noro 1789: int n,i,j,k,l,start,end,len,w;
1.27 noro 1790: int *dw;
1791: struct sparse_weight *sw;
1792: struct weight_or_block *w_or_b;
1793: Obj a0;
1794: NODE a;
1.29 noro 1795: V v,sv,ev;
1796: SYMBOL sym;
1797: int *top;
1.27 noro 1798:
1799: /* l = number of vars in vl */
1800: for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ );
1801: /* n = number of primitives in order */
1802: wb = BDY(order);
1803: n = length(wb);
1804: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1805: spec->id = 3;
1806: spec->obj = (Obj)order;
1807: spec->nv = l;
1808: spec->ord.composite.length = n;
1.28 noro 1809: w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *)
1.29 noro 1810: MALLOC(sizeof(struct weight_or_block)*(n+1));
1811:
1812: /* top : register the top variable in each w_or_b specification */
1813: top = (int *)ALLOCA(l*sizeof(int));
1814: for ( i = 0; i < l; i++ ) top[i] = 0;
1815:
1.28 noro 1816: for ( t = wb, i = 0; t; t = NEXT(t), i++ ) {
1.30 noro 1817: if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST )
1818: error("a list of lists must be specified for the key \"order\"");
1.28 noro 1819: a = BDY((LIST)BDY(t));
1.27 noro 1820: len = length(a);
1821: a0 = (Obj)BDY(a);
1822: if ( !a0 || OID(a0) == O_N ) {
1.28 noro 1823: /* a is a dense weight vector */
1.27 noro 1824: dw = (int *)MALLOC(sizeof(int)*len);
1.30 noro 1825: for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) {
1826: if ( !INT((Q)BDY(p)) )
1827: error("a dense weight vector must be specified as a list of integers");
1.27 noro 1828: dw[j] = QTOS((Q)BDY(p));
1.30 noro 1829: }
1.27 noro 1830: w_or_b[i].type = IS_DENSE_WEIGHT;
1831: w_or_b[i].length = len;
1832: w_or_b[i].body.dense_weight = dw;
1.29 noro 1833:
1834: /* find the top */
1835: for ( k = 0; k < len && !dw[k]; k++ );
1836: if ( k < len ) top[k] = 1;
1837:
1.27 noro 1838: } else if ( OID(a0) == O_P ) {
1.28 noro 1839: /* a is a sparse weight vector */
1840: len >>= 1;
1.27 noro 1841: sw = (struct sparse_weight *)
1842: MALLOC(sizeof(struct sparse_weight)*len);
1843: for ( j = 0, p = a; j < len; j++ ) {
1.30 noro 1844: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
1845: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 1846: v = VR((P)BDY(p)); p = NEXT(p);
1.27 noro 1847: for ( tvl = vl, k = 0; tvl && tvl->v != v;
1848: k++, tvl = NEXT(tvl) );
1849: if ( !tvl )
1.30 noro 1850: error("invalid variable name in a sparse weight vector");
1.27 noro 1851: sw[j].pos = k;
1.30 noro 1852: if ( !INT((Q)BDY(p)) )
1853: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 1854: sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p);
1.27 noro 1855: }
1.37 noro 1856: qsort(sw,len,sizeof(struct sparse_weight),
1857: (int (*)(const void *,const void *))comp_sw);
1.27 noro 1858: w_or_b[i].type = IS_SPARSE_WEIGHT;
1859: w_or_b[i].length = len;
1860: w_or_b[i].body.sparse_weight = sw;
1.29 noro 1861:
1862: /* find the top */
1863: for ( k = 0; k < len && !sw[k].value; k++ );
1864: if ( k < len ) top[sw[k].pos] = 1;
1865: } else if ( OID(a0) == O_RANGE ) {
1866: /* [range(v1,v2),w] */
1867: sv = VR((P)(((RANGE)a0)->start));
1868: ev = VR((P)(((RANGE)a0)->end));
1869: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
1870: if ( !tvl )
1871: error("invalid range");
1872: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
1873: if ( !tvl )
1874: error("invalid range");
1875: len = end-start+1;
1876: sw = (struct sparse_weight *)
1877: MALLOC(sizeof(struct sparse_weight)*len);
1878: w = QTOS((Q)BDY(NEXT(a)));
1879: for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) );
1880: for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) {
1881: sw[j].pos = k;
1882: sw[j].value = w;
1883: }
1884: w_or_b[i].type = IS_SPARSE_WEIGHT;
1885: w_or_b[i].length = len;
1886: w_or_b[i].body.sparse_weight = sw;
1887:
1888: /* register the top */
1889: if ( w ) top[start] = 1;
1.28 noro 1890: } else if ( OID(a0) == O_SYMBOL ) {
1891: /* a is a block */
1.29 noro 1892: sym = (SYMBOL)a0; a = NEXT(a); len--;
1893: if ( OID((Obj)BDY(a)) == O_RANGE ) {
1894: sv = VR((P)(((RANGE)BDY(a))->start));
1895: ev = VR((P)(((RANGE)BDY(a))->end));
1896: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
1897: if ( !tvl )
1898: error("invalid range");
1899: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
1900: if ( !tvl )
1901: error("invalid range");
1902: len = end-start+1;
1903: } else {
1904: for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a));
1.28 noro 1905: tvl = NEXT(tvl), start++ );
1.29 noro 1906: for ( p = NEXT(a), tvl = NEXT(tvl); p;
1.30 noro 1907: p = NEXT(p), tvl = NEXT(tvl) ) {
1908: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
1909: error("a block must be specified as [ordsymbol,var1,var2,...]");
1.29 noro 1910: if ( tvl->v != VR((P)BDY(p)) ) break;
1.30 noro 1911: }
1.29 noro 1912: if ( p )
1.30 noro 1913: error("a block must be contiguous in the variable list");
1.29 noro 1914: }
1.28 noro 1915: w_or_b[i].type = IS_BLOCK;
1916: w_or_b[i].length = len;
1917: w_or_b[i].body.block.start = start;
1918: if ( !strcmp(sym->name,"@grlex") )
1919: w_or_b[i].body.block.order = 0;
1920: else if ( !strcmp(sym->name,"@glex") )
1921: w_or_b[i].body.block.order = 1;
1922: else if ( !strcmp(sym->name,"@lex") )
1923: w_or_b[i].body.block.order = 2;
1924: else
1.29 noro 1925: error("invalid ordername");
1926: /* register the tops */
1927: for ( j = 0, k = start; j < len; j++, k++ )
1928: top[k] = 1;
1.28 noro 1929: }
1.29 noro 1930: }
1931: for ( k = 0; k < l && top[k]; k++ );
1932: if ( k < l ) {
1933: /* incomplete order specification; add @grlex */
1934: w_or_b[n].type = IS_BLOCK;
1935: w_or_b[n].length = l;
1936: w_or_b[n].body.block.start = 0;
1937: w_or_b[n].body.block.order = 0;
1938: spec->ord.composite.length = n+1;
1.27 noro 1939: }
1940: }
1941:
1.35 noro 1942: /* module order spec */
1943:
1944: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s)
1945: {
1946: struct modorder_spec *spec;
1947: NODE n,t;
1948: LIST list;
1949: int *ds;
1950: int i,l;
1951: Q q;
1952:
1953: *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec));
1954: spec->id = id;
1955: if ( shift ) {
1956: n = BDY(shift);
1957: spec->len = l = length(n);
1958: spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int));
1959: for ( t = n, i = 0; t; t = NEXT(t), i++ )
1960: ds[i] = QTOS((Q)BDY(t));
1961: } else {
1962: spec->len = 0;
1963: spec->degree_shift = 0;
1964: }
1965: STOQ(id,q);
1966: n = mknode(2,q,shift);
1967: MKLIST(list,n);
1968: spec->obj = (Obj)list;
1969: }
1970:
1.7 noro 1971: /*
1972: * converters
1973: *
1974: */
1975:
1.20 noro 1976: void dp_homo(DP p,DP *rp)
1.5 noro 1977: {
1.7 noro 1978: MP m,mr,mr0;
1979: int i,n,nv,td;
1980: DL dl,dlh;
1.5 noro 1981:
1.7 noro 1982: if ( !p )
1983: *rp = 0;
1984: else {
1985: n = p->nv; nv = n + 1;
1986: m = BDY(p); td = sugard(m);
1987: for ( mr0 = 0; m; m = NEXT(m) ) {
1988: NEXTMP(mr0,mr); mr->c = m->c;
1989: dl = m->dl;
1990: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1991: dlh->td = td;
1992: for ( i = 0; i < n; i++ )
1993: dlh->d[i] = dl->d[i];
1994: dlh->d[n] = td - dl->td;
1995: }
1996: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 1997: }
1998: }
1999:
1.20 noro 2000: void dp_dehomo(DP p,DP *rp)
1.5 noro 2001: {
1.7 noro 2002: MP m,mr,mr0;
2003: int i,n,nv;
2004: DL dl,dlh;
1.5 noro 2005:
1.7 noro 2006: if ( !p )
2007: *rp = 0;
2008: else {
2009: n = p->nv; nv = n - 1;
2010: m = BDY(p);
2011: for ( mr0 = 0; m; m = NEXT(m) ) {
2012: NEXTMP(mr0,mr); mr->c = m->c;
2013: dlh = m->dl;
2014: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
2015: dl->td = dlh->td - dlh->d[nv];
2016: for ( i = 0; i < nv; i++ )
2017: dl->d[i] = dlh->d[i];
2018: }
2019: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
2020: }
1.5 noro 2021: }
2022:
1.20 noro 2023: void dp_mod(DP p,int mod,NODE subst,DP *rp)
1.5 noro 2024: {
1.7 noro 2025: MP m,mr,mr0;
2026: P t,s,s1;
2027: V v;
2028: NODE tn;
1.5 noro 2029:
1.7 noro 2030: if ( !p )
2031: *rp = 0;
2032: else {
2033: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
2034: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
2035: v = VR((P)BDY(tn)); tn = NEXT(tn);
2036: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
2037: }
2038: ptomp(mod,s,&t);
2039: if ( t ) {
2040: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
2041: }
2042: }
2043: if ( mr0 ) {
2044: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
2045: } else
2046: *rp = 0;
2047: }
1.5 noro 2048: }
2049:
1.20 noro 2050: void dp_rat(DP p,DP *rp)
1.5 noro 2051: {
1.7 noro 2052: MP m,mr,mr0;
1.5 noro 2053:
1.7 noro 2054: if ( !p )
2055: *rp = 0;
2056: else {
2057: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
2058: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5 noro 2059: }
1.7 noro 2060: if ( mr0 ) {
2061: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
2062: } else
2063: *rp = 0;
1.5 noro 2064: }
2065: }
2066:
2067:
1.27 noro 2068: void homogenize_order(struct order_spec *old,int n,struct order_spec **newp)
1.5 noro 2069: {
1.7 noro 2070: struct order_pair *l;
2071: int length,nv,row,i,j;
2072: int **newm,**oldm;
1.27 noro 2073: struct order_spec *new;
1.31 noro 2074: int onv,nnv,nlen,olen,owlen;
2075: struct weight_or_block *owb,*nwb;
1.5 noro 2076:
1.27 noro 2077: *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 2078: switch ( old->id ) {
2079: case 0:
2080: switch ( old->ord.simple ) {
2081: case 0:
2082: new->id = 0; new->ord.simple = 0; break;
2083: case 1:
2084: l = (struct order_pair *)
2085: MALLOC_ATOMIC(2*sizeof(struct order_pair));
2086: l[0].length = n; l[0].order = 1;
2087: l[1].length = 1; l[1].order = 2;
2088: new->id = 1;
2089: new->ord.block.order_pair = l;
2090: new->ord.block.length = 2; new->nv = n+1;
2091: break;
2092: case 2:
2093: new->id = 0; new->ord.simple = 1; break;
2094: case 3: case 4: case 5:
2095: new->id = 0; new->ord.simple = old->ord.simple+3;
2096: dp_nelim = n-1; break;
2097: case 6: case 7: case 8: case 9:
2098: new->id = 0; new->ord.simple = old->ord.simple; break;
2099: default:
2100: error("homogenize_order : invalid input");
2101: }
2102: break;
2103: case 1:
2104: length = old->ord.block.length;
2105: l = (struct order_pair *)
2106: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
2107: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
2108: l[length].order = 2; l[length].length = 1;
2109: new->id = 1; new->nv = n+1;
2110: new->ord.block.order_pair = l;
2111: new->ord.block.length = length+1;
2112: break;
2113: case 2:
2114: nv = old->nv; row = old->ord.matrix.row;
2115: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
2116: for ( i = 0; i <= nv; i++ )
2117: newm[0][i] = 1;
2118: for ( i = 0; i < row; i++ ) {
2119: for ( j = 0; j < nv; j++ )
2120: newm[i+1][j] = oldm[i][j];
2121: newm[i+1][j] = 0;
2122: }
2123: new->id = 2; new->nv = nv+1;
2124: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1.31 noro 2125: break;
2126: case 3:
2127: onv = old->nv;
2128: nnv = onv+1;
2129: olen = old->ord.composite.length;
2130: nlen = olen+1;
2131: owb = old->ord.composite.w_or_b;
2132: nwb = (struct weight_or_block *)
2133: MALLOC(nlen*sizeof(struct weight_or_block));
2134: for ( i = 0; i < olen; i++ ) {
2135: nwb[i].type = owb[i].type;
2136: switch ( owb[i].type ) {
2137: case IS_DENSE_WEIGHT:
2138: owlen = owb[i].length;
2139: nwb[i].length = owlen+1;
2140: nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int));
2141: for ( j = 0; j < owlen; j++ )
2142: nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j];
2143: nwb[i].body.dense_weight[owlen] = 0;
2144: break;
2145: case IS_SPARSE_WEIGHT:
2146: nwb[i].length = owb[i].length;
2147: nwb[i].body.sparse_weight = owb[i].body.sparse_weight;
2148: break;
2149: case IS_BLOCK:
2150: nwb[i].length = owb[i].length;
2151: nwb[i].body.block = owb[i].body.block;
2152: break;
2153: }
2154: }
2155: nwb[i].type = IS_SPARSE_WEIGHT;
2156: nwb[i].body.sparse_weight =
2157: (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
2158: nwb[i].body.sparse_weight[0].pos = onv;
2159: nwb[i].body.sparse_weight[0].value = 1;
2160: new->id = 3;
2161: new->nv = nnv;
2162: new->ord.composite.length = nlen;
2163: new->ord.composite.w_or_b = nwb;
2164: print_composite_order_spec(new);
1.7 noro 2165: break;
2166: default:
2167: error("homogenize_order : invalid input");
1.5 noro 2168: }
1.7 noro 2169: }
2170:
1.20 noro 2171: void qltozl(Q *w,int n,Q *dvr)
1.7 noro 2172: {
2173: N nm,dn;
2174: N g,l1,l2,l3;
2175: Q c,d;
2176: int i;
2177: struct oVECT v;
1.5 noro 2178:
2179: for ( i = 0; i < n; i++ )
1.7 noro 2180: if ( w[i] && !INT(w[i]) )
2181: break;
2182: if ( i == n ) {
2183: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
2184: igcdv(&v,dvr); return;
2185: }
2186: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
2187: for ( i = 1; i < n; i++ ) {
2188: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
2189: gcdn(nm,NM(c),&g); nm = g;
2190: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 2191: }
1.7 noro 2192: if ( UNIN(dn) )
2193: NTOQ(nm,1,d);
2194: else
2195: NDTOQ(nm,dn,1,d);
2196: *dvr = d;
2197: }
1.5 noro 2198:
1.20 noro 2199: int comp_nm(Q *a,Q *b)
1.7 noro 2200: {
2201: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
2202: }
2203:
1.20 noro 2204: void sortbynm(Q *w,int n)
1.7 noro 2205: {
2206: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
2207: }
1.5 noro 2208:
2209:
1.7 noro 2210: /*
2211: * simple operations
2212: *
2213: */
1.5 noro 2214:
1.20 noro 2215: int dp_redble(DP p1,DP p2)
1.7 noro 2216: {
2217: int i,n;
2218: DL d1,d2;
1.5 noro 2219:
1.7 noro 2220: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
2221: if ( d1->td < d2->td )
2222: return 0;
2223: else {
2224: for ( i = 0, n = p1->nv; i < n; i++ )
2225: if ( d1->d[i] < d2->d[i] )
2226: return 0;
2227: return 1;
1.5 noro 2228: }
2229: }
2230:
1.20 noro 2231: void dp_subd(DP p1,DP p2,DP *rp)
1.5 noro 2232: {
1.7 noro 2233: int i,n;
1.5 noro 2234: DL d1,d2,d;
2235: MP m;
1.7 noro 2236: DP s;
1.5 noro 2237:
2238: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 noro 2239: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 2240: for ( i = 0; i < n; i++ )
1.7 noro 2241: d->d[i] = d1->d[i]-d2->d[i];
2242: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
2243: MKDP(n,m,s); s->sugar = d->td;
2244: *rp = s;
2245: }
2246:
1.20 noro 2247: void dltod(DL d,int n,DP *rp)
1.7 noro 2248: {
2249: MP m;
2250: DP s;
2251:
2252: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
2253: MKDP(n,m,s); s->sugar = d->td;
2254: *rp = s;
1.5 noro 2255: }
2256:
1.20 noro 2257: void dp_hm(DP p,DP *rp)
1.5 noro 2258: {
2259: MP m,mr;
2260:
2261: if ( !p )
2262: *rp = 0;
2263: else {
2264: m = BDY(p);
2265: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
2266: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
2267: }
2268: }
2269:
1.35 noro 2270: void dp_ht(DP p,DP *rp)
2271: {
2272: MP m,mr;
2273:
2274: if ( !p )
2275: *rp = 0;
2276: else {
2277: m = BDY(p);
2278: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
2279: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
2280: }
2281: }
2282:
1.20 noro 2283: void dp_rest(DP p,DP *rp)
1.5 noro 2284: {
2285: MP m;
2286:
2287: m = BDY(p);
2288: if ( !NEXT(m) )
2289: *rp = 0;
2290: else {
2291: MKDP(p->nv,NEXT(m),*rp);
2292: if ( *rp )
2293: (*rp)->sugar = p->sugar;
2294: }
2295: }
2296:
1.20 noro 2297: DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl)
1.5 noro 2298: {
1.21 noro 2299: register int i, *d1, *d2, *d, td;
1.5 noro 2300:
2301: if ( !dl ) NEWDL(dl,nv);
2302: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1.21 noro 2303: for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) {
2304: *d = *d1 > *d2 ? *d1 : *d2;
2305: td += MUL_WEIGHT(*d,i);
2306: }
1.5 noro 2307: dl->td = td;
2308: return dl;
2309: }
2310:
1.20 noro 2311: int dl_equal(int nv,DL dl1,DL dl2)
1.5 noro 2312: {
2313: register int *d1, *d2, n;
2314:
2315: if ( dl1->td != dl2->td ) return 0;
2316: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
2317: if ( *d1 != *d2 ) return 0;
2318: return 1;
2319: }
2320:
1.20 noro 2321: int dp_nt(DP p)
1.5 noro 2322: {
2323: int i;
2324: MP m;
2325:
2326: if ( !p )
2327: return 0;
2328: else {
2329: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
2330: return i;
2331: }
2332: }
2333:
1.20 noro 2334: int dp_homogeneous(DP p)
1.15 noro 2335: {
2336: MP m;
2337: int d;
2338:
2339: if ( !p )
2340: return 1;
2341: else {
2342: m = BDY(p);
2343: d = m->dl->td;
2344: m = NEXT(m);
2345: for ( ; m; m = NEXT(m) ) {
2346: if ( m->dl->td != d )
2347: return 0;
2348: }
2349: return 1;
2350: }
1.16 noro 2351: }
2352:
1.20 noro 2353: void _print_mp(int nv,MP m)
1.16 noro 2354: {
2355: int i;
2356:
1.17 noro 2357: if ( !m )
1.16 noro 2358: return;
2359: for ( ; m; m = NEXT(m) ) {
2360: fprintf(stderr,"%d<",ITOS(C(m)));
2361: for ( i = 0; i < nv; i++ ) {
2362: fprintf(stderr,"%d",m->dl->d[i]);
2363: if ( i != nv-1 )
2364: fprintf(stderr," ");
2365: }
2366: fprintf(stderr,">",C(m));
2367: }
2368: fprintf(stderr,"\n");
1.15 noro 2369: }
1.26 noro 2370:
2371: static int cmp_mp_nvar;
2372:
2373: int comp_mp(MP *a,MP *b)
2374: {
2375: return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl);
2376: }
2377:
2378: void dp_sort(DP p,DP *rp)
2379: {
2380: MP t,mp,mp0;
2381: int i,n;
2382: DP r;
2383: MP *w;
2384:
2385: if ( !p ) {
2386: *rp = 0;
2387: return;
2388: }
2389: for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ );
2390: w = (MP *)ALLOCA(n*sizeof(MP));
2391: for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ )
2392: w[i] = t;
2393: cmp_mp_nvar = NV(p);
2394: qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp);
2395: mp0 = 0;
2396: for ( i = n-1; i >= 0; i-- ) {
2397: NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]);
2398: NEXT(mp) = mp0; mp0 = mp;
2399: }
2400: MKDP(p->nv,mp0,r);
2401: r->sugar = p->sugar;
2402: *rp = r;
2403: }
2404:
1.32 noro 2405: DP extract_initial_term_from_dp(DP p,int *weight,int n);
2406: LIST extract_initial_term(LIST f,int *weight,int n);
2407:
2408: DP extract_initial_term_from_dp(DP p,int *weight,int n)
2409: {
1.34 noro 2410: int w,t,i,top;
1.32 noro 2411: MP m,r0,r;
2412: DP dp;
2413:
2414: if ( !p ) return 0;
1.34 noro 2415: top = 1;
1.32 noro 2416: for ( m = BDY(p); m; m = NEXT(m) ) {
2417: for ( i = 0, t = 0; i < n; i++ )
2418: t += weight[i]*m->dl->d[i];
1.34 noro 2419: if ( top || t > w ) {
1.32 noro 2420: r0 = 0;
2421: w = t;
1.34 noro 2422: top = 0;
1.32 noro 2423: }
2424: if ( t == w ) {
2425: NEXTMP(r0,r);
2426: r->dl = m->dl;
2427: r->c = m->c;
2428: }
2429: }
2430: NEXT(r) = 0;
2431: MKDP(p->nv,r0,dp);
2432: return dp;
2433: }
2434:
2435: LIST extract_initial_term(LIST f,int *weight,int n)
2436: {
2437: NODE nd,r0,r;
2438: Obj p;
2439: LIST l;
2440:
2441: nd = BDY(f);
2442: for ( r0 = 0; nd; nd = NEXT(nd) ) {
2443: NEXTNODE(r0,r);
2444: p = (Obj)BDY(nd);
2445: BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n);
2446: }
2447: if ( r0 ) NEXT(r) = 0;
2448: MKLIST(l,r0);
2449: return l;
2450: }
2451:
2452: LIST dp_initial_term(LIST f,struct order_spec *ord)
2453: {
2454: int n,l,i;
2455: struct weight_or_block *worb;
2456: int *weight;
2457:
2458: switch ( ord->id ) {
2459: case 2: /* matrix order */
2460: /* extract the first row */
2461: n = ord->nv;
2462: weight = ord->ord.matrix.matrix[0];
2463: return extract_initial_term(f,weight,n);
2464: case 3: /* composite order */
2465: /* the first w_or_b */
2466: worb = ord->ord.composite.w_or_b;
2467: switch ( worb->type ) {
2468: case IS_DENSE_WEIGHT:
2469: n = worb->length;
2470: weight = worb->body.dense_weight;
2471: return extract_initial_term(f,weight,n);
2472: case IS_SPARSE_WEIGHT:
2473: n = ord->nv;
2474: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 2475: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 2476: l = worb->length;
2477: for ( i = 0; i < l; i++ )
2478: weight[worb->body.sparse_weight[i].pos]
2479: = worb->body.sparse_weight[i].value;
2480: return extract_initial_term(f,weight,n);
2481: default:
2482: error("dp_initial_term : unsupported order");
2483: }
2484: default:
2485: error("dp_initial_term : unsupported order");
2486: }
2487: }
2488:
2489: int highest_order_dp(DP p,int *weight,int n);
2490: LIST highest_order(LIST f,int *weight,int n);
2491:
2492: int highest_order_dp(DP p,int *weight,int n)
2493: {
1.34 noro 2494: int w,t,i,top;
1.32 noro 2495: MP m;
2496:
2497: if ( !p ) return -1;
1.34 noro 2498: top = 1;
1.32 noro 2499: for ( m = BDY(p); m; m = NEXT(m) ) {
2500: for ( i = 0, t = 0; i < n; i++ )
2501: t += weight[i]*m->dl->d[i];
1.34 noro 2502: if ( top || t > w ) {
1.32 noro 2503: w = t;
1.34 noro 2504: top = 0;
2505: }
1.32 noro 2506: }
2507: return w;
2508: }
2509:
2510: LIST highest_order(LIST f,int *weight,int n)
2511: {
2512: int h;
2513: NODE nd,r0,r;
2514: Obj p;
2515: LIST l;
2516: Q q;
2517:
2518: nd = BDY(f);
2519: for ( r0 = 0; nd; nd = NEXT(nd) ) {
2520: NEXTNODE(r0,r);
2521: p = (Obj)BDY(nd);
2522: h = highest_order_dp((DP)p,weight,n);
2523: STOQ(h,q);
2524: BDY(r) = (pointer)q;
2525: }
2526: if ( r0 ) NEXT(r) = 0;
2527: MKLIST(l,r0);
2528: return l;
2529: }
2530:
2531: LIST dp_order(LIST f,struct order_spec *ord)
2532: {
2533: int n,l,i;
2534: struct weight_or_block *worb;
2535: int *weight;
2536:
2537: switch ( ord->id ) {
2538: case 2: /* matrix order */
2539: /* extract the first row */
2540: n = ord->nv;
2541: weight = ord->ord.matrix.matrix[0];
2542: return highest_order(f,weight,n);
2543: case 3: /* composite order */
2544: /* the first w_or_b */
2545: worb = ord->ord.composite.w_or_b;
2546: switch ( worb->type ) {
2547: case IS_DENSE_WEIGHT:
2548: n = worb->length;
2549: weight = worb->body.dense_weight;
2550: return highest_order(f,weight,n);
2551: case IS_SPARSE_WEIGHT:
2552: n = ord->nv;
2553: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 2554: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 2555: l = worb->length;
2556: for ( i = 0; i < l; i++ )
2557: weight[worb->body.sparse_weight[i].pos]
2558: = worb->body.sparse_weight[i].value;
2559: return highest_order(f,weight,n);
2560: default:
2561: error("dp_initial_term : unsupported order");
2562: }
2563: default:
2564: error("dp_initial_term : unsupported order");
1.35 noro 2565: }
2566: }
2567:
2568: int dpv_ht(DPV p,DP *h)
2569: {
2570: int len,max,maxi,i,t;
2571: DP *e;
2572: MP m,mr;
2573:
2574: len = p->len;
2575: e = p->body;
2576: max = -1;
2577: maxi = -1;
2578: for ( i = 0; i < len; i++ )
2579: if ( e[i] && (t = BDY(e[i])->dl->td) > max ) {
2580: max = t;
2581: maxi = i;
2582: }
2583: if ( max < 0 ) {
2584: *h = 0;
2585: return -1;
2586: } else {
2587: m = BDY(e[maxi]);
2588: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
2589: MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td; /* XXX */
2590: return maxi;
1.32 noro 2591: }
2592: }
1.42 noro 2593:
2594: /* return 1 if 0 <_w1 v && v <_w2 0 */
2595:
2596: int in_c12(int n,int *v,int row1,int **w1,int row2, int **w2)
2597: {
2598: int t1,t2;
2599:
2600: t1 = compare_zero(n,v,row1,w1);
2601: t2 = compare_zero(n,v,row2,w2);
2602: if ( t1 > 0 && t2 < 0 ) return 1;
2603: else return 0;
2604: }
2605:
2606: /* 0 < u => 1, 0 > u => -1 */
2607:
2608: int compare_zero(int n,int *u,int row,int **w)
2609: {
2610: int i,j,t;
2611: int *wi;
2612:
2613: for ( i = 0; i < row; i++ ) {
2614: wi = w[i];
2615: for ( j = 0, t = 0; j < n; j++ ) t += u[j]*wi[j];
2616: if ( t > 0 ) return 1;
2617: else if ( t < 0 ) return -1;
2618: }
2619: return 0;
2620: }
2621:
2622: /* functions for generic groebner walk */
2623: /* u=0 means u=-infty */
2624:
2625: int compare_facet_preorder(int n,int *u,int *v,
2626: int row1,int **w1,int row2,int **w2)
2627: {
2628: int i,j,s,t,tu,tv;
2629: int *w2i,*uv;
2630:
2631: if ( !u ) return 1;
2632: uv = W_ALLOC(n);
2633: for ( i = 0; i < row2; i++ ) {
2634: w2i = w2[i];
2635: for ( j = 0, tu = tv = 0; j < n; j++ )
2636: if ( s = w2i[j] ) {
2637: tu += s*u[j]; tv += s*v[j];
2638: }
2639: for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu;
2640: t = compare_zero(n,uv,row1,w1);
2641: if ( t > 0 ) return 1;
2642: else if ( t < 0 ) return 0;
2643: }
2644: return 1;
2645: }
2646:
2647: /* return 0 if last_w = infty */
2648:
2649: NODE compute_last_w(NODE g,NODE gh,int n,int **w,
2650: int row1,int **w1,int row2,int **w2)
2651: {
2652: DP d;
2653: MP f,m0,m;
2654: int *wt,*v,*h;
2655: NODE t,s,n0,tn,n1,r0,r;
2656: int i;
2657:
2658: wt = W_ALLOC(n);
2659: n0 = 0;
2660: for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
2661: f = BDY((DP)BDY(t));
2662: h = BDY((DP)BDY(s))->dl->d;
2663: for ( ; f; f = NEXT(f) ) {
2664: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
2665: for ( i = 0; i < n && !wt[i]; i++ );
2666: if ( i == n ) continue;
2667:
2668: if ( in_c12(n,wt,row1,w1,row2,w2) &&
2669: compare_facet_preorder(n,*w,wt,row1,w1,row2,w2) ) {
2670: v = (int *)MALLOC_ATOMIC(n*sizeof(int));
2671: for ( i = 0; i < n; i++ ) v[i] = wt[i];
2672: MKNODE(n1,v,n0); n0 = n1;
2673: }
2674: }
2675: }
2676: if ( !n0 ) return 0;
2677: for ( t = n0; t; t = NEXT(t) ) {
2678: v = (int *)BDY(t);
2679: for ( s = n0; s; s = NEXT(s) )
2680: if ( !compare_facet_preorder(n,v,(int *)BDY(s),row1,w1,row2,w2) )
2681: break;
2682: if ( !s ) {
2683: *w = v;
2684: break;
2685: }
2686: }
2687: if ( !t )
2688: error("compute_last_w : cannot happen");
2689: r0 = 0;
2690: for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
2691: f = BDY((DP)BDY(t));
2692: h = BDY((DP)BDY(s))->dl->d;
2693: for ( m0 = 0; f; f = NEXT(f) ) {
2694: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
2695: for ( i = 0; i < n && !wt[i]; i++ );
2696: if ( i == n ||
2697: (compare_facet_preorder(n,wt,*w,row1,w1,row2,w2)
2698: && compare_facet_preorder(n,*w,wt,row1,w1,row2,w2)) ) {
2699: NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
2700: }
2701: }
1.43 noro 2702: NEXT(m) = 0;
1.42 noro 2703: MKDP(((DP)BDY(t))->nv,m0,d); d->sugar = ((DP)BDY(t))->sugar;
2704: NEXTNODE(r0,r); BDY(r) = (pointer)d;
2705: }
2706: NEXT(r) = 0;
2707: return r0;
2708: }
1.44 noro 2709:
2710: /* compute a sufficient set of d(f)=u-v */
2711:
2712: NODE compute_essential_df(DP *g,DP *gh,int ng)
2713: {
1.45 ! noro 2714: int nv,i,j,k,t,lj;
! 2715: NODE r,r1,ri,rt,r0;
! 2716: MP m;
! 2717: MP *mj;
! 2718: DL di,hj,dl,dlt;
! 2719: int *d,*dt;
! 2720: LIST l;
1.44 noro 2721: Q q;
1.45 ! noro 2722:
! 2723: nv = g[0]->nv;
! 2724: r = 0;
! 2725: for ( j = 0; j < ng; j++ ) {
! 2726: for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ );
! 2727: mj = (MP *)ALLOCA(lj*sizeof(MP));
! 2728: for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ )
! 2729: mj[k] = m;
! 2730: for ( i = 0; i < lj; i++ ) {
! 2731: for ( di = mj[i]->dl, k = i+1; k < lj; k++ )
! 2732: if ( _dl_redble(di,mj[k]->dl,nv) ) break;
! 2733: if ( k < lj ) mj[i] = 0;
! 2734: }
! 2735: hj = BDY(gh[j])->dl;
! 2736: _NEWDL(dl,nv); d = dl->d;
! 2737: r0 = r;
! 2738: for ( i = 0; i < lj; i++ ) {
! 2739: if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) {
! 2740: for ( k = 0, t = 0; k < nv; k++ ) {
! 2741: d[k] = hj->d[k]-di->d[k];
! 2742: t += d[k];
! 2743: }
! 2744: dl->td = t;
! 2745: #if 1
! 2746: for ( rt = r0; rt; rt = NEXT(rt) ) {
! 2747: dlt = (DL)BDY(rt);
! 2748: if ( dlt->td != dl->td ) continue;
! 2749: for ( dt = dlt->d, k = 0; k < nv; k++ )
! 2750: if ( d[k] != dt[k] ) break;
! 2751: if ( k == nv ) break;
! 2752: }
! 2753: #else
! 2754: rt = 0;
! 2755: #endif
! 2756: if ( !rt ) {
! 2757: MKNODE(r1,dl,r); r = r1;
! 2758: _NEWDL(dl,nv); d = dl->d;
! 2759: }
1.44 noro 2760: }
2761: }
2762: }
1.45 ! noro 2763: for ( rt = r; rt; rt = NEXT(rt) ) {
! 2764: dl = (DL)BDY(rt); d = dl->d;
! 2765: ri = 0;
! 2766: for ( k = nv-1; k >= 0; k-- ) {
! 2767: STOQ(d[k],q);
! 2768: MKNODE(r1,q,ri); ri = r1;
1.44 noro 2769: }
1.45 ! noro 2770: MKNODE(r1,0,ri); MKLIST(l,r1);
! 2771: BDY(rt) = (pointer)l;
1.44 noro 2772: }
2773: return r;
2774: }