![[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.19 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.18 2001/09/17 10:32: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:
65: /*
66: * content reduction
67: *
68: */
69:
70: void dp_ptozp(p,rp)
71: DP p,*rp;
72: {
73: MP m,mr,mr0;
74: int i,n;
75: Q *w;
76: Q dvr;
77: P t;
78:
79: if ( !p )
80: *rp = 0;
81: else {
82: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
83: w = (Q *)ALLOCA(n*sizeof(Q));
84: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
85: if ( NUM(m->c) )
86: w[i] = (Q)m->c;
87: else
88: ptozp(m->c,1,&w[i],&t);
89: sortbynm(w,n);
90: qltozl(w,n,&dvr);
91: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
92: NEXTMP(mr0,mr); divsp(CO,m->c,(P)dvr,&mr->c); mr->dl = m->dl;
93: }
94: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
95: }
96: }
97:
98: void dp_ptozp2(p0,p1,hp,rp)
99: DP p0,p1;
100: DP *hp,*rp;
101: {
102: DP t,s,h,r;
103: MP m,mr,mr0,m0;
104:
105: addd(CO,p0,p1,&t); dp_ptozp(t,&s);
106: if ( !p0 ) {
107: h = 0; r = s;
108: } else if ( !p1 ) {
109: h = s; r = 0;
110: } else {
111: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
112: m = NEXT(m), m0 = NEXT(m0) ) {
113: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
114: }
115: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
116: }
117: if ( h )
118: h->sugar = p0->sugar;
119: if ( r )
120: r->sugar = p1->sugar;
121: *hp = h; *rp = r;
122: }
1.1 noro 123:
124: void dp_idiv(p,c,rp)
125: DP p;
126: Q c;
127: DP *rp;
128: {
129: Q t;
130: N nm,q;
131: int sgn,s;
132: MP mr0,m,mr;
133:
134: if ( !p )
135: *rp = 0;
136: else if ( MUNIQ((Q)c) )
137: *rp = p;
138: else if ( MUNIQ((Q)c) )
139: chsgnd(p,rp);
140: else {
141: nm = NM(c); sgn = SGN(c);
142: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
143: NEXTMP(mr0,mr);
144:
145: divsn(NM((Q)(m->c)),nm,&q);
146: s = sgn*SGN((Q)(m->c));
147: NTOQ(q,s,t);
148: mr->c = (P)t;
149: mr->dl = m->dl;
150: }
151: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
152: if ( *rp )
153: (*rp)->sugar = p->sugar;
154: }
155: }
156:
157: void dp_mbase(hlist,mbase)
158: NODE hlist;
159: NODE *mbase;
160: {
161: DL *dl;
162: DL d;
163: int i,j,n,nvar,td;
164:
165: n = length(hlist); nvar = ((DP)BDY(hlist))->nv;
166: dl = (DL *)MALLOC(n*sizeof(DL));
167: for ( i = 0; i < n; i++, hlist = NEXT(hlist) )
168: dl[i] = BDY((DP)BDY(hlist))->dl;
169: NEWDL(d,nvar); *mbase = 0;
170: while ( 1 ) {
171: insert_to_node(d,mbase,nvar);
172: for ( i = nvar-1; i >= 0; ) {
173: d->d[i]++; d->td++;
174: for ( j = 0; j < n; j++ ) {
175: if ( _dl_redble(dl[j],d,nvar) )
176: break;
177: }
178: if ( j < n ) {
179: for ( j = nvar-1; j >= i; j-- )
180: d->d[j] = 0;
181: for ( j = 0, td = 0; j < i; j++ )
182: td += d->d[j];
183: d->td = td;
184: i--;
185: } else
186: break;
187: }
188: if ( i < 0 )
189: break;
190: }
191: }
192:
193: int _dl_redble(d1,d2,nvar)
194: DL d1,d2;
195: int nvar;
196: {
197: int i;
198:
199: if ( d1->td > d2->td )
200: return 0;
201: for ( i = 0; i < nvar; i++ )
202: if ( d1->d[i] > d2->d[i] )
203: break;
204: if ( i < nvar )
205: return 0;
206: else
207: return 1;
208: }
209:
210: void insert_to_node(d,n,nvar)
211: DL d;
212: NODE *n;
213: int nvar;
214: {
215: DL d1;
216: MP m;
217: DP dp;
218: NODE n0,n1,n2;
219:
220: NEWDL(d1,nvar); d1->td = d->td;
221: bcopy((char *)d->d,(char *)d1->d,nvar*sizeof(int));
222: NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0;
223: MKDP(nvar,m,dp); dp->sugar = d->td;
224: if ( !(*n) ) {
225: MKNODE(n1,dp,0); *n = n1;
226: } else {
227: for ( n1 = *n, n0 = 0; n1; n0 = n1, n1 = NEXT(n1) )
228: if ( (*cmpdl)(nvar,d,BDY((DP)BDY(n1))->dl) > 0 ) {
229: MKNODE(n2,dp,n1);
230: if ( !n0 )
231: *n = n2;
232: else
233: NEXT(n0) = n2;
234: break;
235: }
236: if ( !n1 ) {
237: MKNODE(n2,dp,0); NEXT(n0) = n2;
238: }
239: }
240: }
241:
242: void dp_vtod(c,p,rp)
243: Q *c;
244: DP p;
245: DP *rp;
246: {
247: MP mr0,m,mr;
248: int i;
249:
250: if ( !p )
251: *rp = 0;
252: else {
253: for ( mr0 = 0, m = BDY(p), i = 0; m; m = NEXT(m), i++ ) {
254: NEXTMP(mr0,mr); mr->c = (P)c[i]; mr->dl = m->dl;
255: }
256: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
257: (*rp)->sugar = p->sugar;
258: }
259: }
260:
1.8 noro 261: extern int mpi_mag;
262: extern int PCoeffs;
263:
264: void dp_ptozp_d(p,rp)
1.1 noro 265: DP p,*rp;
266: {
267: int i,j,k,l,n,nsep;
268: MP m;
269: NODE tn,n0,n1,n2,n3;
270: struct oVECT v;
271: VECT c,cs;
272: VECT qi,ri;
273: LIST *qr;
274: int s,id;
275: Obj dmy;
276: Q d0,d1,gcd,a,u,u1;
277: Q *q,*r;
278: STRING iqr_v;
279: pointer *b;
280: N qn,gn;
281: double get_rtime();
282: int blen;
1.8 noro 283: NODE dist;
284: int ndist;
1.1 noro 285: double t0;
286: double t_e,t_d,t_d1,t_c;
1.8 noro 287: extern int DP_NFStat;
288: extern LIST Dist;
1.1 noro 289:
290: if ( !p )
291: *rp = 0;
292: else {
1.8 noro 293: if ( PCoeffs ) {
294: dp_ptozp(p,rp); return;
295: }
1.9 noro 296: if ( !Dist || p_mag(BDY(p)->c) <= mpi_mag ) {
1.8 noro 297: dist = 0; ndist = 0;
298: if ( DP_NFStat ) fprintf(asir_out,"L");
299: } else {
300: dist = BDY(Dist); ndist = length(dist);
301: if ( DP_NFStat ) fprintf(asir_out,"D");
302: }
1.1 noro 303: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
304: nsep = ndist + 1;
305: if ( n <= nsep ) {
306: dp_ptozp(p,rp); return;
307: }
308: t0 = get_rtime();
309: dp_dtov(p,&c);
310: igcdv_estimate(c,&d0);
311: t_e = get_rtime()-t0;
312: t0 = get_rtime();
313: dp_dtov(p,&c);
314: sepvect(c,nsep,&cs);
315: MKSTR(iqr_v,"iqr");
316: qr = (LIST *)CALLOC(nsep,sizeof(LIST));
317: q = (Q *)CALLOC(n,sizeof(Q));
318: r = (Q *)CALLOC(n,sizeof(Q));
319: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
320: MKNODE(n3,d0,0); MKNODE(n2,b[i],n3);
321: MKNODE(n1,iqr_v,n2); MKNODE(n0,BDY(tn),n1);
322: Pox_rpc(n0,&dmy);
323: }
324: iqrv(b[i],d0,&qr[i]);
325: dp_dtov(p,&c);
326: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
327: Pox_pop_local(tn,&qr[i]);
328: if ( OID(qr[i]) == O_ERR ) {
329: printexpr(CO,(Obj)qr[i]);
330: error("dp_ptozp_d : aborted");
331: }
332: }
333: t_d = get_rtime()-t0;
334: t_d1 = t_d/n;
335: t0 = get_rtime();
336: for ( i = j = 0; i < nsep; i++ ) {
337: tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
338: for ( k = 0, l = qi->len; k < l; k++, j++ ) {
339: q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
340: }
341: }
342: v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
343: if ( d1 ) {
344: gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
345: divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
346: for ( i = 0; i < n; i++ ) {
347: mulq(a,q[i],&u);
348: if ( r[i] ) {
349: divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
350: addq(u,u1,&q[i]);
351: } else
352: q[i] = u;
353: }
354: } else
355: gcd = d0;
356: dp_vtod(q,p,rp);
357: t_c = get_rtime()-t0;
358: blen=p_mag((P)gcd);
359: pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
360: if ( 0 )
361: fprintf(stderr,"(%d,%d)",p_mag((P)d0)-blen,blen);
362: }
363: }
364:
1.8 noro 365: void dp_ptozp2_d(p0,p1,hp,rp)
1.1 noro 366: DP p0,p1;
367: DP *hp,*rp;
368: {
369: DP t,s,h,r;
370: MP m,mr,mr0,m0;
371:
1.8 noro 372: addd(CO,p0,p1,&t); dp_ptozp_d(t,&s);
1.1 noro 373: if ( !p0 ) {
374: h = 0; r = s;
375: } else if ( !p1 ) {
376: h = s; r = 0;
377: } else {
378: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
379: m = NEXT(m), m0 = NEXT(m0) ) {
380: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
381: }
382: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
383: }
384: if ( h )
385: h->sugar = p0->sugar;
386: if ( r )
387: r->sugar = p1->sugar;
388: *hp = h; *rp = r;
1.5 noro 389: }
390:
1.7 noro 391: void dp_prim(p,rp)
392: DP p,*rp;
1.5 noro 393: {
1.7 noro 394: P t,g;
395: DP p1;
396: MP m,mr,mr0;
397: int i,n;
398: P *w;
399: Q *c;
400: Q dvr;
1.5 noro 401:
1.7 noro 402: if ( !p )
403: *rp = 0;
404: else if ( dp_fcoeffs )
405: *rp = p;
406: else if ( NoGCD )
407: dp_ptozp(p,rp);
408: else {
409: dp_ptozp(p,&p1); p = p1;
410: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
411: if ( n == 1 ) {
412: m = BDY(p);
413: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
414: MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
415: return;
416: }
417: w = (P *)ALLOCA(n*sizeof(P));
418: c = (Q *)ALLOCA(n*sizeof(Q));
419: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
420: if ( NUM(m->c) ) {
421: c[i] = (Q)m->c; w[i] = (P)ONE;
422: } else
423: ptozp(m->c,1,&c[i],&w[i]);
424: qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
425: if ( NUM(g) )
426: *rp = p;
427: else {
428: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
429: NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
430: }
431: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 432: }
1.7 noro 433: }
1.5 noro 434: }
435:
436: void heu_nezgcdnpz(vl,pl,m,pr)
437: VL vl;
438: P *pl,*pr;
439: int m;
440: {
441: int i,r;
442: P gcd,t,s1,s2,u;
443: Q rq;
444:
445: while ( 1 ) {
446: for ( i = 0, s1 = 0; i < m; i++ ) {
447: r = random(); UTOQ(r,rq);
448: mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
449: }
450: for ( i = 0, s2 = 0; i < m; i++ ) {
451: r = random(); UTOQ(r,rq);
452: mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
453: }
454: ezgcdp(vl,s1,s2,&gcd);
455: for ( i = 0; i < m; i++ ) {
456: if ( !divtpz(vl,pl[i],gcd,&t) )
457: break;
458: }
459: if ( i == m )
460: break;
461: }
462: *pr = gcd;
463: }
464:
465: void dp_prim_mod(p,mod,rp)
466: int mod;
467: DP p,*rp;
468: {
469: P t,g;
470: MP m,mr,mr0;
471:
472: if ( !p )
473: *rp = 0;
474: else if ( NoGCD )
475: *rp = p;
476: else {
477: for ( m = BDY(p), g = m->c, m = NEXT(m); m; m = NEXT(m) ) {
478: gcdprsmp(CO,mod,g,m->c,&t); g = t;
479: }
480: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
481: NEXTMP(mr0,mr); divsmp(CO,mod,m->c,g,&mr->c); mr->dl = m->dl;
482: }
483: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
484: }
485: }
486:
1.7 noro 487: void dp_cont(p,rp)
1.5 noro 488: DP p;
1.7 noro 489: Q *rp;
1.5 noro 490: {
1.7 noro 491: VECT v;
1.5 noro 492:
1.7 noro 493: dp_dtov(p,&v); igcdv(v,rp);
1.5 noro 494: }
495:
1.7 noro 496: void dp_dtov(dp,rp)
497: DP dp;
498: VECT *rp;
1.5 noro 499: {
1.7 noro 500: MP m,t;
501: int i,n;
502: VECT v;
503: pointer *p;
1.5 noro 504:
1.7 noro 505: m = BDY(dp);
506: for ( t = m, n = 0; t; t = NEXT(t), n++ );
507: MKVECT(v,n);
508: for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
509: p[i] = (pointer)(t->c);
510: *rp = v;
1.5 noro 511: }
512:
1.7 noro 513: /*
514: * s-poly computation
515: *
516: */
1.5 noro 517:
1.7 noro 518: void dp_sp(p1,p2,rp)
519: DP p1,p2;
1.5 noro 520: DP *rp;
521: {
1.7 noro 522: int i,n,td;
523: int *w;
524: DL d1,d2,d;
525: MP m;
526: DP t,s1,s2,u;
527: Q c,c1,c2;
528: N gn,tn;
1.5 noro 529:
1.7 noro 530: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
531: w = (int *)ALLOCA(n*sizeof(int));
532: for ( i = 0, td = 0; i < n; i++ ) {
533: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
1.5 noro 534: }
1.7 noro 535:
536: NEWDL(d,n); d->td = td - d1->td;
537: for ( i = 0; i < n; i++ )
538: d->d[i] = w[i] - d1->d[i];
539: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
540: if ( INT(c1) && INT(c2) ) {
541: gcdn(NM(c1),NM(c2),&gn);
542: if ( !UNIN(gn) ) {
543: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
544: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1.5 noro 545: }
546: }
1.7 noro 547:
548: NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
549: MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
550:
551: NEWDL(d,n); d->td = td - d2->td;
552: for ( i = 0; i < n; i++ )
553: d->d[i] = w[i] - d2->d[i];
554: NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0;
555: MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
556:
557: subd(CO,t,u,rp);
1.14 noro 558: if ( GenTrace ) {
559: LIST hist;
560: NODE node;
561:
562: node = mknode(4,ONE,0,s1,ONE);
563: MKLIST(hist,node);
564: MKNODE(TraceList,hist,0);
565:
566: node = mknode(4,ONE,0,0,ONE);
567: chsgnd(s2,(DP *)&ARG2(node));
568: MKLIST(hist,node);
569: MKNODE(node,hist,TraceList); TraceList = node;
570: }
571: }
572:
573: void _dp_sp_dup(p1,p2,rp)
574: DP p1,p2;
575: DP *rp;
576: {
577: int i,n,td;
578: int *w;
579: DL d1,d2,d;
580: MP m;
581: DP t,s1,s2,u;
582: Q c,c1,c2;
583: N gn,tn;
584:
585: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
586: w = (int *)ALLOCA(n*sizeof(int));
587: for ( i = 0, td = 0; i < n; i++ ) {
588: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
589: }
590:
591: _NEWDL(d,n); d->td = td - d1->td;
592: for ( i = 0; i < n; i++ )
593: d->d[i] = w[i] - d1->d[i];
594: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
595: if ( INT(c1) && INT(c2) ) {
596: gcdn(NM(c1),NM(c2),&gn);
597: if ( !UNIN(gn) ) {
598: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
599: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
600: }
601: }
602:
603: _NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
604: _MKDP(n,m,s1); s1->sugar = d->td; _muld_dup(CO,s1,p1,&t); _free_dp(s1);
605:
606: _NEWDL(d,n); d->td = td - d2->td;
607: for ( i = 0; i < n; i++ )
608: d->d[i] = w[i] - d2->d[i];
609: _NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0;
610: _MKDP(n,m,s2); s2->sugar = d->td; _muld_dup(CO,s2,p2,&u); _free_dp(s2);
611:
612: _addd_destructive(CO,t,u,rp);
1.7 noro 613: if ( GenTrace ) {
614: LIST hist;
615: NODE node;
616:
617: node = mknode(4,ONE,0,s1,ONE);
618: MKLIST(hist,node);
619: MKNODE(TraceList,hist,0);
620:
621: node = mknode(4,ONE,0,0,ONE);
622: chsgnd(s2,(DP *)&ARG2(node));
623: MKLIST(hist,node);
624: MKNODE(node,hist,TraceList); TraceList = node;
625: }
626: }
627:
628: void dp_sp_mod(p1,p2,mod,rp)
629: DP p1,p2;
630: int mod;
631: DP *rp;
632: {
633: int i,n,td;
634: int *w;
635: DL d1,d2,d;
636: MP m;
637: DP t,s,u;
638:
639: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
640: w = (int *)ALLOCA(n*sizeof(int));
641: for ( i = 0, td = 0; i < n; i++ ) {
642: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
643: }
1.18 noro 644: NEWDL_NOINIT(d,n); d->td = td - d1->td;
1.7 noro 645: for ( i = 0; i < n; i++ )
646: d->d[i] = w[i] - d1->d[i];
647: NEWMP(m); m->dl = d; m->c = (P)BDY(p2)->c; NEXT(m) = 0;
648: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
1.18 noro 649: NEWDL_NOINIT(d,n); d->td = td - d2->td;
1.7 noro 650: for ( i = 0; i < n; i++ )
651: d->d[i] = w[i] - d2->d[i];
652: NEWMP(m); m->dl = d; m->c = (P)BDY(p1)->c; NEXT(m) = 0;
653: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
654: submd(CO,mod,t,u,rp);
655: }
656:
657: void _dp_sp_mod_dup(p1,p2,mod,rp)
658: DP p1,p2;
659: int mod;
660: DP *rp;
661: {
662: int i,n,td;
663: int *w;
664: DL d1,d2,d;
665: MP m;
666: DP t,s,u;
667:
668: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
669: w = (int *)ALLOCA(n*sizeof(int));
670: for ( i = 0, td = 0; i < n; i++ ) {
671: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
672: }
673: _NEWDL(d,n); d->td = td - d1->td;
674: for ( i = 0; i < n; i++ )
675: d->d[i] = w[i] - d1->d[i];
676: _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
677: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
678: _NEWDL(d,n); d->td = td - d2->td;
679: for ( i = 0; i < n; i++ )
680: d->d[i] = w[i] - d2->d[i];
681: _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
682: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
683: _addmd_destructive(mod,t,u,rp);
684: }
685:
686: void _dp_sp_mod(p1,p2,mod,rp)
687: DP p1,p2;
688: int mod;
689: DP *rp;
690: {
691: int i,n,td;
692: int *w;
693: DL d1,d2,d;
694: MP m;
695: DP t,s,u;
696:
697: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
698: w = (int *)ALLOCA(n*sizeof(int));
699: for ( i = 0, td = 0; i < n; i++ ) {
700: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
701: }
702: NEWDL(d,n); d->td = td - d1->td;
703: for ( i = 0; i < n; i++ )
704: d->d[i] = w[i] - d1->d[i];
705: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
706: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
707: NEWDL(d,n); d->td = td - d2->td;
708: for ( i = 0; i < n; i++ )
709: d->d[i] = w[i] - d2->d[i];
710: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
711: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
712: addmd_destructive(mod,t,u,rp);
713: }
714:
715: /*
716: * m-reduction
1.13 noro 717: * do content reduction over Z or Q(x,...)
718: * do nothing over finite fields
1.7 noro 719: *
720: */
721:
722: void dp_red(p0,p1,p2,head,rest,dnp,multp)
723: DP p0,p1,p2;
724: DP *head,*rest;
725: P *dnp;
726: DP *multp;
727: {
728: int i,n;
729: DL d1,d2,d;
730: MP m;
731: DP t,s,r,h;
732: Q c,c1,c2;
733: N gn,tn;
734: P g,a;
735:
736: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
737: NEWDL(d,n); d->td = d1->td - d2->td;
738: for ( i = 0; i < n; i++ )
739: d->d[i] = d1->d[i]-d2->d[i];
740: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
741: if ( dp_fcoeffs ) {
742: /* do nothing */
743: } else if ( INT(c1) && INT(c2) ) {
744: gcdn(NM(c1),NM(c2),&gn);
745: if ( !UNIN(gn) ) {
746: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
747: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
748: }
749: } else {
750: ezgcdpz(CO,(P)c1,(P)c2,&g);
751: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
752: }
753: NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
754: *multp = s;
755: muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
756: muldc(CO,p0,(P)c2,&h);
757: *head = h; *rest = r; *dnp = (P)c2;
758: }
759:
1.13 noro 760: /* m-reduction over a field */
761:
762: void dp_red_f(p1,p2,rest)
763: DP p1,p2;
764: DP *rest;
765: {
766: int i,n;
767: DL d1,d2,d;
768: MP m;
769: DP t,s,r,h;
770: Obj a,b;
771:
772: n = p1->nv;
773: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
774:
775: NEWDL(d,n); d->td = d1->td - d2->td;
776: for ( i = 0; i < n; i++ )
777: d->d[i] = d1->d[i]-d2->d[i];
778:
779: NEWMP(m); m->dl = d;
780: divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
781: C(m) = (P)b;
782: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
783:
784: muld(CO,s,p2,&t); addd(CO,p1,t,rest);
785: }
786:
1.7 noro 787: void dp_red_mod(p0,p1,p2,mod,head,rest,dnp)
788: DP p0,p1,p2;
789: int mod;
790: DP *head,*rest;
791: P *dnp;
792: {
793: int i,n;
794: DL d1,d2,d;
795: MP m;
796: DP t,s,r,h;
797: P c1,c2,g,u;
798:
799: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
800: NEWDL(d,n); d->td = d1->td - d2->td;
801: for ( i = 0; i < n; i++ )
802: d->d[i] = d1->d[i]-d2->d[i];
803: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
804: gcdprsmp(CO,mod,c1,c2,&g);
805: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
806: if ( NUM(c2) ) {
807: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
808: }
809: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1.11 noro 810: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7 noro 811: if ( NUM(c2) ) {
812: addmd(CO,mod,p1,t,&r); h = p0;
813: } else {
814: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
815: }
816: *head = h; *rest = r; *dnp = c2;
817: }
818:
1.10 noro 819: struct oEGT eg_red_mod;
820:
1.7 noro 821: void _dp_red_mod_destructive(p1,p2,mod,rp)
822: DP p1,p2;
823: int mod;
824: DP *rp;
825: {
826: int i,n;
827: DL d1,d2,d;
828: MP m;
829: DP t,s;
1.16 noro 830: int c,c1,c2;
831: struct oEGT t0,t1;
832: extern int do_weyl;
1.7 noro 833:
834: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
835: _NEWDL(d,n); d->td = d1->td - d2->td;
836: for ( i = 0; i < n; i++ )
837: d->d[i] = d1->d[i]-d2->d[i];
1.16 noro 838: c = invm(ITOS(BDY(p2)->c),mod);
839: c2 = ITOS(BDY(p1)->c);
840: DMAR(c,c2,0,mod,c1);
1.7 noro 841: _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
1.16 noro 842: #if 0
1.7 noro 843: _MKDP(n,m,s); s->sugar = d->td;
844: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 845: #else
846: if ( do_weyl ) {
1.19 ! noro 847: _MKDP(n,m,s); s->sugar = d->td;
! 848: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 849: } else {
850: _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
851: }
852: #endif
1.10 noro 853: /* get_eg(&t0); */
1.7 noro 854: _addmd_destructive(mod,p1,t,rp);
1.10 noro 855: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 856: }
857:
858: /*
859: * normal form computation
860: *
861: */
1.5 noro 862:
863: void dp_true_nf(b,g,ps,full,rp,dnp)
864: NODE b;
865: DP g;
866: DP *ps;
867: int full;
868: DP *rp;
869: P *dnp;
870: {
871: DP u,p,d,s,t,dmy;
872: NODE l;
873: MP m,mr;
874: int i,n;
875: int *wb;
876: int sugar,psugar;
877: P dn,tdn,tdn1;
878:
879: dn = (P)ONE;
880: if ( !g ) {
881: *rp = 0; *dnp = dn; return;
882: }
883: for ( n = 0, l = b; l; l = NEXT(l), n++ );
884: wb = (int *)ALLOCA(n*sizeof(int));
885: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
886: wb[i] = QTOS((Q)BDY(l));
887: sugar = g->sugar;
888: for ( d = 0; g; ) {
889: for ( u = 0, i = 0; i < n; i++ ) {
890: if ( dp_redble(g,p = ps[wb[i]]) ) {
891: dp_red(d,g,p,&t,&u,&tdn,&dmy);
892: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
893: sugar = MAX(sugar,psugar);
894: if ( !u ) {
895: if ( d )
896: d->sugar = sugar;
897: *rp = d; *dnp = dn; return;
898: } else {
899: d = t;
900: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
901: }
902: break;
903: }
904: }
905: if ( u )
906: g = u;
907: else if ( !full ) {
908: if ( g ) {
909: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
910: }
911: *rp = g; *dnp = dn; return;
912: } else {
913: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
914: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
915: addd(CO,d,t,&s); d = s;
916: dp_rest(g,&t); g = t;
917: }
918: }
919: if ( d )
920: d->sugar = sugar;
921: *rp = d; *dnp = dn;
922: }
923:
1.13 noro 924: /* nf computation over Z */
925:
926: void dp_nf_z(b,g,ps,full,multiple,rp)
1.5 noro 927: NODE b;
928: DP g;
929: DP *ps;
930: int full,multiple;
931: DP *rp;
932: {
933: DP u,p,d,s,t,dmy1;
934: P dmy;
935: NODE l;
936: MP m,mr;
937: int i,n;
938: int *wb;
939: int hmag;
940: int sugar,psugar;
941:
942: if ( !g ) {
943: *rp = 0; return;
944: }
945: for ( n = 0, l = b; l; l = NEXT(l), n++ );
946: wb = (int *)ALLOCA(n*sizeof(int));
947: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
948: wb[i] = QTOS((Q)BDY(l));
1.12 noro 949:
1.13 noro 950: hmag = multiple*HMAG(g);
1.5 noro 951: sugar = g->sugar;
1.12 noro 952:
1.5 noro 953: for ( d = 0; g; ) {
954: for ( u = 0, i = 0; i < n; i++ ) {
955: if ( dp_redble(g,p = ps[wb[i]]) ) {
956: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
957: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
958: sugar = MAX(sugar,psugar);
959: if ( !u ) {
960: if ( d )
961: d->sugar = sugar;
962: *rp = d; return;
963: }
964: d = t;
965: break;
966: }
967: }
968: if ( u ) {
969: g = u;
970: if ( d ) {
1.13 noro 971: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 972: dp_ptozp2(d,g,&t,&u); d = t; g = u;
973: hmag = multiple*HMAG(d);
974: }
975: } else {
1.13 noro 976: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 977: dp_ptozp(g,&t); g = t;
978: hmag = multiple*HMAG(g);
979: }
980: }
981: }
982: else if ( !full ) {
983: if ( g ) {
984: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
985: }
986: *rp = g; return;
987: } else {
988: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
989: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
990: addd(CO,d,t,&s); d = s;
991: dp_rest(g,&t); g = t;
992:
993: }
994: }
995: if ( d )
996: d->sugar = sugar;
997: *rp = d;
998: }
999:
1.13 noro 1000: /* nf computation over a field */
1001:
1002: void dp_nf_f(b,g,ps,full,rp)
1003: NODE b;
1004: DP g;
1005: DP *ps;
1006: int full;
1007: DP *rp;
1008: {
1009: DP u,p,d,s,t;
1010: P dmy;
1011: NODE l;
1012: MP m,mr;
1013: int i,n;
1014: int *wb;
1015: int sugar,psugar;
1016:
1017: if ( !g ) {
1018: *rp = 0; return;
1019: }
1020: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1021: wb = (int *)ALLOCA(n*sizeof(int));
1022: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1023: wb[i] = QTOS((Q)BDY(l));
1024:
1025: sugar = g->sugar;
1026: for ( d = 0; g; ) {
1027: for ( u = 0, i = 0; i < n; i++ ) {
1028: if ( dp_redble(g,p = ps[wb[i]]) ) {
1029: dp_red_f(g,p,&u);
1030: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1031: sugar = MAX(sugar,psugar);
1032: if ( !u ) {
1033: if ( d )
1034: d->sugar = sugar;
1035: *rp = d; return;
1036: }
1037: break;
1038: }
1039: }
1040: if ( u )
1041: g = u;
1042: else if ( !full ) {
1043: if ( g ) {
1044: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1045: }
1046: *rp = g; return;
1047: } else {
1048: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1049: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1050: addd(CO,d,t,&s); d = s;
1051: dp_rest(g,&t); g = t;
1052: }
1053: }
1054: if ( d )
1055: d->sugar = sugar;
1056: *rp = d;
1057: }
1058:
1059: /* nf computation over GF(mod) (only for internal use) */
1060:
1.5 noro 1061: void dp_nf_mod(b,g,ps,mod,full,rp)
1062: NODE b;
1063: DP g;
1064: DP *ps;
1065: int mod,full;
1066: DP *rp;
1067: {
1068: DP u,p,d,s,t;
1069: P dmy;
1070: NODE l;
1071: MP m,mr;
1072: int sugar,psugar;
1073:
1074: if ( !g ) {
1075: *rp = 0; return;
1076: }
1077: sugar = g->sugar;
1078: for ( d = 0; g; ) {
1079: for ( u = 0, l = b; l; l = NEXT(l) ) {
1080: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1081: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
1082: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1083: sugar = MAX(sugar,psugar);
1084: if ( !u ) {
1085: if ( d )
1086: d->sugar = sugar;
1087: *rp = d; return;
1088: }
1089: d = t;
1090: break;
1091: }
1092: }
1093: if ( u )
1094: g = u;
1095: else if ( !full ) {
1096: if ( g ) {
1097: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1098: }
1099: *rp = g; return;
1100: } else {
1101: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1102: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1103: addmd(CO,mod,d,t,&s); d = s;
1104: dp_rest(g,&t); g = t;
1105: }
1106: }
1107: if ( d )
1108: d->sugar = sugar;
1109: *rp = d;
1110: }
1111:
1112: void dp_true_nf_mod(b,g,ps,mod,full,rp,dnp)
1113: NODE b;
1114: DP g;
1115: DP *ps;
1116: int mod,full;
1117: DP *rp;
1118: P *dnp;
1119: {
1120: DP u,p,d,s,t;
1121: NODE l;
1122: MP m,mr;
1123: int i,n;
1124: int *wb;
1125: int sugar,psugar;
1126: P dn,tdn,tdn1;
1127:
1128: dn = (P)ONEM;
1129: if ( !g ) {
1130: *rp = 0; *dnp = dn; return;
1131: }
1132: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1133: wb = (int *)ALLOCA(n*sizeof(int));
1134: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1135: wb[i] = QTOS((Q)BDY(l));
1136: sugar = g->sugar;
1137: for ( d = 0; g; ) {
1138: for ( u = 0, i = 0; i < n; i++ ) {
1139: if ( dp_redble(g,p = ps[wb[i]]) ) {
1140: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
1141: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1142: sugar = MAX(sugar,psugar);
1143: if ( !u ) {
1144: if ( d )
1145: d->sugar = sugar;
1146: *rp = d; *dnp = dn; return;
1147: } else {
1148: d = t;
1149: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1150: }
1151: break;
1152: }
1153: }
1154: if ( u )
1155: g = u;
1156: else if ( !full ) {
1157: if ( g ) {
1158: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1159: }
1160: *rp = g; *dnp = dn; return;
1161: } else {
1162: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1163: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1164: addmd(CO,mod,d,t,&s); d = s;
1165: dp_rest(g,&t); g = t;
1166: }
1167: }
1168: if ( d )
1169: d->sugar = sugar;
1170: *rp = d; *dnp = dn;
1171: }
1172:
1.7 noro 1173: void _dp_nf_mod_destructive(b,g,ps,mod,full,rp)
1174: NODE b;
1175: DP g;
1176: DP *ps;
1177: int mod,full;
1178: DP *rp;
1.5 noro 1179: {
1.7 noro 1180: DP u,p,d,s,t;
1181: NODE l;
1182: MP m,mr,mrd;
1183: int sugar,psugar,n,h_reducible,i;
1.5 noro 1184:
1.7 noro 1185: if ( !g ) {
1186: *rp = 0; return;
1.5 noro 1187: }
1.7 noro 1188: sugar = g->sugar;
1189: n = g->nv;
1190: for ( d = 0; g; ) {
1191: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
1192: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1193: h_reducible = 1;
1194: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1195: _dp_red_mod_destructive(g,p,mod,&u); g = u;
1196: sugar = MAX(sugar,psugar);
1197: if ( !g ) {
1198: if ( d )
1199: d->sugar = sugar;
1200: _dptodp(d,rp); _free_dp(d); return;
1201: }
1202: break;
1203: }
1204: }
1205: if ( !h_reducible ) {
1206: /* head term is not reducible */
1207: if ( !full ) {
1208: if ( g )
1209: g->sugar = sugar;
1210: _dptodp(g,rp); _free_dp(g); return;
1211: } else {
1212: m = BDY(g);
1213: if ( NEXT(m) ) {
1214: BDY(g) = NEXT(m); NEXT(m) = 0;
1215: } else {
1216: _FREEDP(g); g = 0;
1217: }
1218: if ( d ) {
1219: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1220: NEXT(mrd) = m;
1221: } else {
1222: _MKDP(n,m,d);
1223: }
1224: }
1225: }
1.5 noro 1226: }
1.7 noro 1227: if ( d )
1228: d->sugar = sugar;
1229: _dptodp(d,rp); _free_dp(d);
1.5 noro 1230: }
1.13 noro 1231:
1232: /* reduction by linear base over a field */
1233:
1234: void dp_lnf_f(p1,p2,g,r1p,r2p)
1235: DP p1,p2;
1236: NODE g;
1237: DP *r1p,*r2p;
1238: {
1239: DP r1,r2,b1,b2,t,s;
1240: Obj c,c1,c2;
1241: NODE l,b;
1242: int n;
1243:
1244: if ( !p1 ) {
1245: *r1p = p1; *r2p = p2; return;
1246: }
1247: n = p1->nv;
1248: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1249: if ( !r1 ) {
1250: *r1p = r1; *r2p = r2; return;
1251: }
1252: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1253: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1254: b2 = (DP)BDY(NEXT(b));
1255: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
1256: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
1257: muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
1258: muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
1259: }
1260: }
1261: *r1p = r1; *r2p = r2;
1262: }
1263:
1264: /* reduction by linear base over GF(mod) */
1.5 noro 1265:
1.7 noro 1266: void dp_lnf_mod(p1,p2,g,mod,r1p,r2p)
1267: DP p1,p2;
1268: NODE g;
1269: int mod;
1270: DP *r1p,*r2p;
1.5 noro 1271: {
1.7 noro 1272: DP r1,r2,b1,b2,t,s;
1273: P c;
1274: MQ c1,c2;
1275: NODE l,b;
1276: int n;
1277:
1278: if ( !p1 ) {
1279: *r1p = p1; *r2p = p2; return;
1280: }
1281: n = p1->nv;
1282: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1283: if ( !r1 ) {
1284: *r1p = r1; *r2p = r2; return;
1285: }
1286: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1287: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1288: b2 = (DP)BDY(NEXT(b));
1289: invmq(mod,(MQ)BDY(b1)->c,&c1);
1290: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
1291: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
1292: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
1293: }
1294: }
1295: *r1p = r1; *r2p = r2;
1.5 noro 1296: }
1297:
1.7 noro 1298: void dp_nf_tab_mod(p,tab,mod,rp)
1299: DP p;
1300: LIST *tab;
1301: int mod;
1302: DP *rp;
1.5 noro 1303: {
1.7 noro 1304: DP s,t,u;
1305: MP m;
1306: DL h;
1307: int i,n;
1308:
1309: if ( !p ) {
1310: *rp = p; return;
1311: }
1312: n = p->nv;
1313: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1314: h = m->dl;
1315: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1316: i++;
1317: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1318: addmd(CO,mod,s,t,&u); s = u;
1319: }
1320: *rp = s;
1.5 noro 1321: }
1322:
1.7 noro 1323: /*
1324: * setting flags
1325: *
1326: */
1327:
1328: int create_order_spec(obj,spec)
1329: Obj obj;
1330: struct order_spec *spec;
1.5 noro 1331: {
1.7 noro 1332: int i,j,n,s,row,col;
1333: struct order_pair *l;
1334: NODE node,t,tn;
1335: MAT m;
1336: pointer **b;
1337: int **w;
1.5 noro 1338:
1.7 noro 1339: if ( !obj || NUM(obj) ) {
1340: spec->id = 0; spec->obj = obj;
1341: spec->ord.simple = QTOS((Q)obj);
1342: return 1;
1343: } else if ( OID(obj) == O_LIST ) {
1344: node = BDY((LIST)obj);
1345: for ( n = 0, t = node; t; t = NEXT(t), n++ );
1346: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1347: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
1348: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
1349: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
1350: s += l[i].length;
1351: }
1352: spec->id = 1; spec->obj = obj;
1353: spec->ord.block.order_pair = l;
1354: spec->ord.block.length = n; spec->nv = s;
1355: return 1;
1356: } else if ( OID(obj) == O_MAT ) {
1357: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
1358: w = almat(row,col);
1359: for ( i = 0; i < row; i++ )
1360: for ( j = 0; j < col; j++ )
1361: w[i][j] = QTOS((Q)b[i][j]);
1362: spec->id = 2; spec->obj = obj;
1363: spec->nv = col; spec->ord.matrix.row = row;
1364: spec->ord.matrix.matrix = w;
1365: return 1;
1366: } else
1.5 noro 1367: return 0;
1368: }
1369:
1.7 noro 1370: /*
1371: * converters
1372: *
1373: */
1374:
1375: void dp_homo(p,rp)
1376: DP p;
1377: DP *rp;
1.5 noro 1378: {
1.7 noro 1379: MP m,mr,mr0;
1380: int i,n,nv,td;
1381: DL dl,dlh;
1.5 noro 1382:
1.7 noro 1383: if ( !p )
1384: *rp = 0;
1385: else {
1386: n = p->nv; nv = n + 1;
1387: m = BDY(p); td = sugard(m);
1388: for ( mr0 = 0; m; m = NEXT(m) ) {
1389: NEXTMP(mr0,mr); mr->c = m->c;
1390: dl = m->dl;
1391: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1392: dlh->td = td;
1393: for ( i = 0; i < n; i++ )
1394: dlh->d[i] = dl->d[i];
1395: dlh->d[n] = td - dl->td;
1396: }
1397: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 1398: }
1399: }
1400:
1.7 noro 1401: void dp_dehomo(p,rp)
1402: DP p;
1.5 noro 1403: DP *rp;
1404: {
1.7 noro 1405: MP m,mr,mr0;
1406: int i,n,nv;
1407: DL dl,dlh;
1.5 noro 1408:
1.7 noro 1409: if ( !p )
1410: *rp = 0;
1411: else {
1412: n = p->nv; nv = n - 1;
1413: m = BDY(p);
1414: for ( mr0 = 0; m; m = NEXT(m) ) {
1415: NEXTMP(mr0,mr); mr->c = m->c;
1416: dlh = m->dl;
1417: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1418: dl->td = dlh->td - dlh->d[nv];
1419: for ( i = 0; i < nv; i++ )
1420: dl->d[i] = dlh->d[i];
1421: }
1422: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1423: }
1.5 noro 1424: }
1425:
1.7 noro 1426: void dp_mod(p,mod,subst,rp)
1427: DP p;
1428: int mod;
1429: NODE subst;
1.5 noro 1430: DP *rp;
1431: {
1.7 noro 1432: MP m,mr,mr0;
1433: P t,s,s1;
1434: V v;
1435: NODE tn;
1.5 noro 1436:
1.7 noro 1437: if ( !p )
1438: *rp = 0;
1439: else {
1440: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1441: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
1442: v = VR((P)BDY(tn)); tn = NEXT(tn);
1443: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
1444: }
1445: ptomp(mod,s,&t);
1446: if ( t ) {
1447: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
1448: }
1449: }
1450: if ( mr0 ) {
1451: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1452: } else
1453: *rp = 0;
1454: }
1.5 noro 1455: }
1456:
1.7 noro 1457: void dp_rat(p,rp)
1458: DP p;
1459: DP *rp;
1.5 noro 1460: {
1.7 noro 1461: MP m,mr,mr0;
1.5 noro 1462:
1.7 noro 1463: if ( !p )
1464: *rp = 0;
1465: else {
1466: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1467: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5 noro 1468: }
1.7 noro 1469: if ( mr0 ) {
1470: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1471: } else
1472: *rp = 0;
1.5 noro 1473: }
1474: }
1475:
1476:
1.7 noro 1477: void homogenize_order(old,n,new)
1478: struct order_spec *old,*new;
1479: int n;
1.5 noro 1480: {
1.7 noro 1481: struct order_pair *l;
1482: int length,nv,row,i,j;
1483: int **newm,**oldm;
1.5 noro 1484:
1.7 noro 1485: switch ( old->id ) {
1486: case 0:
1487: switch ( old->ord.simple ) {
1488: case 0:
1489: new->id = 0; new->ord.simple = 0; break;
1490: case 1:
1491: l = (struct order_pair *)
1492: MALLOC_ATOMIC(2*sizeof(struct order_pair));
1493: l[0].length = n; l[0].order = 1;
1494: l[1].length = 1; l[1].order = 2;
1495: new->id = 1;
1496: new->ord.block.order_pair = l;
1497: new->ord.block.length = 2; new->nv = n+1;
1498: break;
1499: case 2:
1500: new->id = 0; new->ord.simple = 1; break;
1501: case 3: case 4: case 5:
1502: new->id = 0; new->ord.simple = old->ord.simple+3;
1503: dp_nelim = n-1; break;
1504: case 6: case 7: case 8: case 9:
1505: new->id = 0; new->ord.simple = old->ord.simple; break;
1506: default:
1507: error("homogenize_order : invalid input");
1508: }
1509: break;
1510: case 1:
1511: length = old->ord.block.length;
1512: l = (struct order_pair *)
1513: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
1514: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
1515: l[length].order = 2; l[length].length = 1;
1516: new->id = 1; new->nv = n+1;
1517: new->ord.block.order_pair = l;
1518: new->ord.block.length = length+1;
1519: break;
1520: case 2:
1521: nv = old->nv; row = old->ord.matrix.row;
1522: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
1523: for ( i = 0; i <= nv; i++ )
1524: newm[0][i] = 1;
1525: for ( i = 0; i < row; i++ ) {
1526: for ( j = 0; j < nv; j++ )
1527: newm[i+1][j] = oldm[i][j];
1528: newm[i+1][j] = 0;
1529: }
1530: new->id = 2; new->nv = nv+1;
1531: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1532: break;
1533: default:
1534: error("homogenize_order : invalid input");
1.5 noro 1535: }
1.7 noro 1536: }
1537:
1538: void qltozl(w,n,dvr)
1539: Q *w,*dvr;
1540: int n;
1541: {
1542: N nm,dn;
1543: N g,l1,l2,l3;
1544: Q c,d;
1545: int i;
1546: struct oVECT v;
1.5 noro 1547:
1548: for ( i = 0; i < n; i++ )
1.7 noro 1549: if ( w[i] && !INT(w[i]) )
1550: break;
1551: if ( i == n ) {
1552: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
1553: igcdv(&v,dvr); return;
1554: }
1555: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
1556: for ( i = 1; i < n; i++ ) {
1557: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
1558: gcdn(nm,NM(c),&g); nm = g;
1559: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 1560: }
1.7 noro 1561: if ( UNIN(dn) )
1562: NTOQ(nm,1,d);
1563: else
1564: NDTOQ(nm,dn,1,d);
1565: *dvr = d;
1566: }
1.5 noro 1567:
1.7 noro 1568: int comp_nm(a,b)
1569: Q *a,*b;
1570: {
1571: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
1572: }
1573:
1574: void sortbynm(w,n)
1575: Q *w;
1576: int n;
1577: {
1578: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
1579: }
1.5 noro 1580:
1581:
1.7 noro 1582: /*
1583: * simple operations
1584: *
1585: */
1.5 noro 1586:
1.7 noro 1587: int dp_redble(p1,p2)
1588: DP p1,p2;
1589: {
1590: int i,n;
1591: DL d1,d2;
1.5 noro 1592:
1.7 noro 1593: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1594: if ( d1->td < d2->td )
1595: return 0;
1596: else {
1597: for ( i = 0, n = p1->nv; i < n; i++ )
1598: if ( d1->d[i] < d2->d[i] )
1599: return 0;
1600: return 1;
1.5 noro 1601: }
1602: }
1603:
1.7 noro 1604: void dp_subd(p1,p2,rp)
1.5 noro 1605: DP p1,p2;
1606: DP *rp;
1607: {
1.7 noro 1608: int i,n;
1.5 noro 1609: DL d1,d2,d;
1610: MP m;
1.7 noro 1611: DP s;
1.5 noro 1612:
1613: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 noro 1614: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 1615: for ( i = 0; i < n; i++ )
1.7 noro 1616: d->d[i] = d1->d[i]-d2->d[i];
1617: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1618: MKDP(n,m,s); s->sugar = d->td;
1619: *rp = s;
1620: }
1621:
1622: void dltod(d,n,rp)
1623: DL d;
1624: int n;
1625: DP *rp;
1626: {
1627: MP m;
1628: DP s;
1629:
1630: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1631: MKDP(n,m,s); s->sugar = d->td;
1632: *rp = s;
1.5 noro 1633: }
1634:
1635: void dp_hm(p,rp)
1636: DP p;
1637: DP *rp;
1638: {
1639: MP m,mr;
1640:
1641: if ( !p )
1642: *rp = 0;
1643: else {
1644: m = BDY(p);
1645: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
1646: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1647: }
1648: }
1649:
1650: void dp_rest(p,rp)
1651: DP p,*rp;
1652: {
1653: MP m;
1654:
1655: m = BDY(p);
1656: if ( !NEXT(m) )
1657: *rp = 0;
1658: else {
1659: MKDP(p->nv,NEXT(m),*rp);
1660: if ( *rp )
1661: (*rp)->sugar = p->sugar;
1662: }
1663: }
1664:
1665: DL lcm_of_DL(nv,dl1,dl2,dl)
1666: int nv;
1667: DL dl1,dl2;
1668: register DL dl;
1669: {
1670: register int n, *d1, *d2, *d, td;
1671:
1672: if ( !dl ) NEWDL(dl,nv);
1673: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1674: for ( td = 0, n = nv; --n >= 0; d1++, d2++, d++ )
1675: td += (*d = *d1 > *d2 ? *d1 : *d2 );
1676: dl->td = td;
1677: return dl;
1678: }
1679:
1680: int dl_equal(nv,dl1,dl2)
1681: int nv;
1682: DL dl1, dl2;
1683: {
1684: register int *d1, *d2, n;
1685:
1686: if ( dl1->td != dl2->td ) return 0;
1687: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
1688: if ( *d1 != *d2 ) return 0;
1689: return 1;
1690: }
1691:
1692: int dp_nt(p)
1693: DP p;
1694: {
1695: int i;
1696: MP m;
1697:
1698: if ( !p )
1699: return 0;
1700: else {
1701: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
1702: return i;
1703: }
1704: }
1705:
1.15 noro 1706: int dp_homogeneous(p)
1707: DP p;
1708: {
1709: MP m;
1710: int d;
1711:
1712: if ( !p )
1713: return 1;
1714: else {
1715: m = BDY(p);
1716: d = m->dl->td;
1717: m = NEXT(m);
1718: for ( ; m; m = NEXT(m) ) {
1719: if ( m->dl->td != d )
1720: return 0;
1721: }
1722: return 1;
1723: }
1.16 noro 1724: }
1725:
1726: _print_mp(nv,m)
1727: int nv;
1728: MP m;
1729: {
1730: int i;
1731:
1.17 noro 1732: if ( !m )
1.16 noro 1733: return;
1734: for ( ; m; m = NEXT(m) ) {
1735: fprintf(stderr,"%d<",ITOS(C(m)));
1736: for ( i = 0; i < nv; i++ ) {
1737: fprintf(stderr,"%d",m->dl->d[i]);
1738: if ( i != nv-1 )
1739: fprintf(stderr," ");
1740: }
1741: fprintf(stderr,">",C(m));
1742: }
1743: fprintf(stderr,"\n");
1.15 noro 1744: }