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