![[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.17 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.16 2001/09/17 01:18:34 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: }
644: NEWDL(d,n); d->td = td - d1->td;
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);
649: NEWDL(d,n); d->td = td - d2->td;
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 ) {
847: _weyl_mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
848: } else {
849: _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
850: }
851: #endif
1.10 noro 852: /* get_eg(&t0); */
1.7 noro 853: _addmd_destructive(mod,p1,t,rp);
1.10 noro 854: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 855: }
856:
857: /*
858: * normal form computation
859: *
860: */
1.5 noro 861:
862: void dp_true_nf(b,g,ps,full,rp,dnp)
863: NODE b;
864: DP g;
865: DP *ps;
866: int full;
867: DP *rp;
868: P *dnp;
869: {
870: DP u,p,d,s,t,dmy;
871: NODE l;
872: MP m,mr;
873: int i,n;
874: int *wb;
875: int sugar,psugar;
876: P dn,tdn,tdn1;
877:
878: dn = (P)ONE;
879: if ( !g ) {
880: *rp = 0; *dnp = dn; return;
881: }
882: for ( n = 0, l = b; l; l = NEXT(l), n++ );
883: wb = (int *)ALLOCA(n*sizeof(int));
884: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
885: wb[i] = QTOS((Q)BDY(l));
886: sugar = g->sugar;
887: for ( d = 0; g; ) {
888: for ( u = 0, i = 0; i < n; i++ ) {
889: if ( dp_redble(g,p = ps[wb[i]]) ) {
890: dp_red(d,g,p,&t,&u,&tdn,&dmy);
891: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
892: sugar = MAX(sugar,psugar);
893: if ( !u ) {
894: if ( d )
895: d->sugar = sugar;
896: *rp = d; *dnp = dn; return;
897: } else {
898: d = t;
899: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
900: }
901: break;
902: }
903: }
904: if ( u )
905: g = u;
906: else if ( !full ) {
907: if ( g ) {
908: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
909: }
910: *rp = g; *dnp = dn; return;
911: } else {
912: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
913: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
914: addd(CO,d,t,&s); d = s;
915: dp_rest(g,&t); g = t;
916: }
917: }
918: if ( d )
919: d->sugar = sugar;
920: *rp = d; *dnp = dn;
921: }
922:
1.13 noro 923: /* nf computation over Z */
924:
925: void dp_nf_z(b,g,ps,full,multiple,rp)
1.5 noro 926: NODE b;
927: DP g;
928: DP *ps;
929: int full,multiple;
930: DP *rp;
931: {
932: DP u,p,d,s,t,dmy1;
933: P dmy;
934: NODE l;
935: MP m,mr;
936: int i,n;
937: int *wb;
938: int hmag;
939: int sugar,psugar;
940:
941: if ( !g ) {
942: *rp = 0; return;
943: }
944: for ( n = 0, l = b; l; l = NEXT(l), n++ );
945: wb = (int *)ALLOCA(n*sizeof(int));
946: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
947: wb[i] = QTOS((Q)BDY(l));
1.12 noro 948:
1.13 noro 949: hmag = multiple*HMAG(g);
1.5 noro 950: sugar = g->sugar;
1.12 noro 951:
1.5 noro 952: for ( d = 0; g; ) {
953: for ( u = 0, i = 0; i < n; i++ ) {
954: if ( dp_redble(g,p = ps[wb[i]]) ) {
955: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
956: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
957: sugar = MAX(sugar,psugar);
958: if ( !u ) {
959: if ( d )
960: d->sugar = sugar;
961: *rp = d; return;
962: }
963: d = t;
964: break;
965: }
966: }
967: if ( u ) {
968: g = u;
969: if ( d ) {
1.13 noro 970: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 971: dp_ptozp2(d,g,&t,&u); d = t; g = u;
972: hmag = multiple*HMAG(d);
973: }
974: } else {
1.13 noro 975: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 976: dp_ptozp(g,&t); g = t;
977: hmag = multiple*HMAG(g);
978: }
979: }
980: }
981: else if ( !full ) {
982: if ( g ) {
983: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
984: }
985: *rp = g; return;
986: } else {
987: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
988: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
989: addd(CO,d,t,&s); d = s;
990: dp_rest(g,&t); g = t;
991:
992: }
993: }
994: if ( d )
995: d->sugar = sugar;
996: *rp = d;
997: }
998:
1.13 noro 999: /* nf computation over a field */
1000:
1001: void dp_nf_f(b,g,ps,full,rp)
1002: NODE b;
1003: DP g;
1004: DP *ps;
1005: int full;
1006: DP *rp;
1007: {
1008: DP u,p,d,s,t;
1009: P dmy;
1010: NODE l;
1011: MP m,mr;
1012: int i,n;
1013: int *wb;
1014: int sugar,psugar;
1015:
1016: if ( !g ) {
1017: *rp = 0; return;
1018: }
1019: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1020: wb = (int *)ALLOCA(n*sizeof(int));
1021: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1022: wb[i] = QTOS((Q)BDY(l));
1023:
1024: sugar = g->sugar;
1025: for ( d = 0; g; ) {
1026: for ( u = 0, i = 0; i < n; i++ ) {
1027: if ( dp_redble(g,p = ps[wb[i]]) ) {
1028: dp_red_f(g,p,&u);
1029: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1030: sugar = MAX(sugar,psugar);
1031: if ( !u ) {
1032: if ( d )
1033: d->sugar = sugar;
1034: *rp = d; return;
1035: }
1036: break;
1037: }
1038: }
1039: if ( u )
1040: g = u;
1041: else if ( !full ) {
1042: if ( g ) {
1043: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1044: }
1045: *rp = g; return;
1046: } else {
1047: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1048: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1049: addd(CO,d,t,&s); d = s;
1050: dp_rest(g,&t); g = t;
1051: }
1052: }
1053: if ( d )
1054: d->sugar = sugar;
1055: *rp = d;
1056: }
1057:
1058: /* nf computation over GF(mod) (only for internal use) */
1059:
1.5 noro 1060: void dp_nf_mod(b,g,ps,mod,full,rp)
1061: NODE b;
1062: DP g;
1063: DP *ps;
1064: int mod,full;
1065: DP *rp;
1066: {
1067: DP u,p,d,s,t;
1068: P dmy;
1069: NODE l;
1070: MP m,mr;
1071: int sugar,psugar;
1072:
1073: if ( !g ) {
1074: *rp = 0; return;
1075: }
1076: sugar = g->sugar;
1077: for ( d = 0; g; ) {
1078: for ( u = 0, l = b; l; l = NEXT(l) ) {
1079: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1080: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
1081: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1082: sugar = MAX(sugar,psugar);
1083: if ( !u ) {
1084: if ( d )
1085: d->sugar = sugar;
1086: *rp = d; return;
1087: }
1088: d = t;
1089: break;
1090: }
1091: }
1092: if ( u )
1093: g = u;
1094: else if ( !full ) {
1095: if ( g ) {
1096: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1097: }
1098: *rp = g; return;
1099: } else {
1100: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1101: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1102: addmd(CO,mod,d,t,&s); d = s;
1103: dp_rest(g,&t); g = t;
1104: }
1105: }
1106: if ( d )
1107: d->sugar = sugar;
1108: *rp = d;
1109: }
1110:
1111: void dp_true_nf_mod(b,g,ps,mod,full,rp,dnp)
1112: NODE b;
1113: DP g;
1114: DP *ps;
1115: int mod,full;
1116: DP *rp;
1117: P *dnp;
1118: {
1119: DP u,p,d,s,t;
1120: NODE l;
1121: MP m,mr;
1122: int i,n;
1123: int *wb;
1124: int sugar,psugar;
1125: P dn,tdn,tdn1;
1126:
1127: dn = (P)ONEM;
1128: if ( !g ) {
1129: *rp = 0; *dnp = dn; return;
1130: }
1131: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1132: wb = (int *)ALLOCA(n*sizeof(int));
1133: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1134: wb[i] = QTOS((Q)BDY(l));
1135: sugar = g->sugar;
1136: for ( d = 0; g; ) {
1137: for ( u = 0, i = 0; i < n; i++ ) {
1138: if ( dp_redble(g,p = ps[wb[i]]) ) {
1139: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
1140: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1141: sugar = MAX(sugar,psugar);
1142: if ( !u ) {
1143: if ( d )
1144: d->sugar = sugar;
1145: *rp = d; *dnp = dn; return;
1146: } else {
1147: d = t;
1148: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1149: }
1150: break;
1151: }
1152: }
1153: if ( u )
1154: g = u;
1155: else if ( !full ) {
1156: if ( g ) {
1157: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1158: }
1159: *rp = g; *dnp = dn; return;
1160: } else {
1161: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1162: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1163: addmd(CO,mod,d,t,&s); d = s;
1164: dp_rest(g,&t); g = t;
1165: }
1166: }
1167: if ( d )
1168: d->sugar = sugar;
1169: *rp = d; *dnp = dn;
1170: }
1171:
1.7 noro 1172: void _dp_nf_mod_destructive(b,g,ps,mod,full,rp)
1173: NODE b;
1174: DP g;
1175: DP *ps;
1176: int mod,full;
1177: DP *rp;
1.5 noro 1178: {
1.7 noro 1179: DP u,p,d,s,t;
1180: NODE l;
1181: MP m,mr,mrd;
1182: int sugar,psugar,n,h_reducible,i;
1.5 noro 1183:
1.7 noro 1184: if ( !g ) {
1185: *rp = 0; return;
1.5 noro 1186: }
1.7 noro 1187: sugar = g->sugar;
1188: n = g->nv;
1189: for ( d = 0; g; ) {
1190: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
1191: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1192: h_reducible = 1;
1193: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1194: _dp_red_mod_destructive(g,p,mod,&u); g = u;
1195: sugar = MAX(sugar,psugar);
1196: if ( !g ) {
1197: if ( d )
1198: d->sugar = sugar;
1199: _dptodp(d,rp); _free_dp(d); return;
1200: }
1201: break;
1202: }
1203: }
1204: if ( !h_reducible ) {
1205: /* head term is not reducible */
1206: if ( !full ) {
1207: if ( g )
1208: g->sugar = sugar;
1209: _dptodp(g,rp); _free_dp(g); return;
1210: } else {
1211: m = BDY(g);
1212: if ( NEXT(m) ) {
1213: BDY(g) = NEXT(m); NEXT(m) = 0;
1214: } else {
1215: _FREEDP(g); g = 0;
1216: }
1217: if ( d ) {
1218: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1219: NEXT(mrd) = m;
1220: } else {
1221: _MKDP(n,m,d);
1222: }
1223: }
1224: }
1.5 noro 1225: }
1.7 noro 1226: if ( d )
1227: d->sugar = sugar;
1228: _dptodp(d,rp); _free_dp(d);
1.5 noro 1229: }
1.13 noro 1230:
1231: /* reduction by linear base over a field */
1232:
1233: void dp_lnf_f(p1,p2,g,r1p,r2p)
1234: DP p1,p2;
1235: NODE g;
1236: DP *r1p,*r2p;
1237: {
1238: DP r1,r2,b1,b2,t,s;
1239: Obj c,c1,c2;
1240: NODE l,b;
1241: int n;
1242:
1243: if ( !p1 ) {
1244: *r1p = p1; *r2p = p2; return;
1245: }
1246: n = p1->nv;
1247: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1248: if ( !r1 ) {
1249: *r1p = r1; *r2p = r2; return;
1250: }
1251: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1252: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1253: b2 = (DP)BDY(NEXT(b));
1254: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
1255: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
1256: muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
1257: muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
1258: }
1259: }
1260: *r1p = r1; *r2p = r2;
1261: }
1262:
1263: /* reduction by linear base over GF(mod) */
1.5 noro 1264:
1.7 noro 1265: void dp_lnf_mod(p1,p2,g,mod,r1p,r2p)
1266: DP p1,p2;
1267: NODE g;
1268: int mod;
1269: DP *r1p,*r2p;
1.5 noro 1270: {
1.7 noro 1271: DP r1,r2,b1,b2,t,s;
1272: P c;
1273: MQ c1,c2;
1274: NODE l,b;
1275: int n;
1276:
1277: if ( !p1 ) {
1278: *r1p = p1; *r2p = p2; return;
1279: }
1280: n = p1->nv;
1281: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1282: if ( !r1 ) {
1283: *r1p = r1; *r2p = r2; return;
1284: }
1285: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1286: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1287: b2 = (DP)BDY(NEXT(b));
1288: invmq(mod,(MQ)BDY(b1)->c,&c1);
1289: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
1290: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
1291: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
1292: }
1293: }
1294: *r1p = r1; *r2p = r2;
1.5 noro 1295: }
1296:
1.7 noro 1297: void dp_nf_tab_mod(p,tab,mod,rp)
1298: DP p;
1299: LIST *tab;
1300: int mod;
1301: DP *rp;
1.5 noro 1302: {
1.7 noro 1303: DP s,t,u;
1304: MP m;
1305: DL h;
1306: int i,n;
1307:
1308: if ( !p ) {
1309: *rp = p; return;
1310: }
1311: n = p->nv;
1312: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1313: h = m->dl;
1314: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1315: i++;
1316: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1317: addmd(CO,mod,s,t,&u); s = u;
1318: }
1319: *rp = s;
1.5 noro 1320: }
1321:
1.7 noro 1322: /*
1323: * setting flags
1324: *
1325: */
1326:
1327: int create_order_spec(obj,spec)
1328: Obj obj;
1329: struct order_spec *spec;
1.5 noro 1330: {
1.7 noro 1331: int i,j,n,s,row,col;
1332: struct order_pair *l;
1333: NODE node,t,tn;
1334: MAT m;
1335: pointer **b;
1336: int **w;
1.5 noro 1337:
1.7 noro 1338: if ( !obj || NUM(obj) ) {
1339: spec->id = 0; spec->obj = obj;
1340: spec->ord.simple = QTOS((Q)obj);
1341: return 1;
1342: } else if ( OID(obj) == O_LIST ) {
1343: node = BDY((LIST)obj);
1344: for ( n = 0, t = node; t; t = NEXT(t), n++ );
1345: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1346: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
1347: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
1348: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
1349: s += l[i].length;
1350: }
1351: spec->id = 1; spec->obj = obj;
1352: spec->ord.block.order_pair = l;
1353: spec->ord.block.length = n; spec->nv = s;
1354: return 1;
1355: } else if ( OID(obj) == O_MAT ) {
1356: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
1357: w = almat(row,col);
1358: for ( i = 0; i < row; i++ )
1359: for ( j = 0; j < col; j++ )
1360: w[i][j] = QTOS((Q)b[i][j]);
1361: spec->id = 2; spec->obj = obj;
1362: spec->nv = col; spec->ord.matrix.row = row;
1363: spec->ord.matrix.matrix = w;
1364: return 1;
1365: } else
1.5 noro 1366: return 0;
1367: }
1368:
1.7 noro 1369: /*
1370: * converters
1371: *
1372: */
1373:
1374: void dp_homo(p,rp)
1375: DP p;
1376: DP *rp;
1.5 noro 1377: {
1.7 noro 1378: MP m,mr,mr0;
1379: int i,n,nv,td;
1380: DL dl,dlh;
1.5 noro 1381:
1.7 noro 1382: if ( !p )
1383: *rp = 0;
1384: else {
1385: n = p->nv; nv = n + 1;
1386: m = BDY(p); td = sugard(m);
1387: for ( mr0 = 0; m; m = NEXT(m) ) {
1388: NEXTMP(mr0,mr); mr->c = m->c;
1389: dl = m->dl;
1390: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1391: dlh->td = td;
1392: for ( i = 0; i < n; i++ )
1393: dlh->d[i] = dl->d[i];
1394: dlh->d[n] = td - dl->td;
1395: }
1396: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 1397: }
1398: }
1399:
1.7 noro 1400: void dp_dehomo(p,rp)
1401: DP p;
1.5 noro 1402: DP *rp;
1403: {
1.7 noro 1404: MP m,mr,mr0;
1405: int i,n,nv;
1406: DL dl,dlh;
1.5 noro 1407:
1.7 noro 1408: if ( !p )
1409: *rp = 0;
1410: else {
1411: n = p->nv; nv = n - 1;
1412: m = BDY(p);
1413: for ( mr0 = 0; m; m = NEXT(m) ) {
1414: NEXTMP(mr0,mr); mr->c = m->c;
1415: dlh = m->dl;
1416: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1417: dl->td = dlh->td - dlh->d[nv];
1418: for ( i = 0; i < nv; i++ )
1419: dl->d[i] = dlh->d[i];
1420: }
1421: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1422: }
1.5 noro 1423: }
1424:
1.7 noro 1425: void dp_mod(p,mod,subst,rp)
1426: DP p;
1427: int mod;
1428: NODE subst;
1.5 noro 1429: DP *rp;
1430: {
1.7 noro 1431: MP m,mr,mr0;
1432: P t,s,s1;
1433: V v;
1434: NODE tn;
1.5 noro 1435:
1.7 noro 1436: if ( !p )
1437: *rp = 0;
1438: else {
1439: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1440: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
1441: v = VR((P)BDY(tn)); tn = NEXT(tn);
1442: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
1443: }
1444: ptomp(mod,s,&t);
1445: if ( t ) {
1446: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
1447: }
1448: }
1449: if ( mr0 ) {
1450: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1451: } else
1452: *rp = 0;
1453: }
1.5 noro 1454: }
1455:
1.7 noro 1456: void dp_rat(p,rp)
1457: DP p;
1458: DP *rp;
1.5 noro 1459: {
1.7 noro 1460: MP m,mr,mr0;
1.5 noro 1461:
1.7 noro 1462: if ( !p )
1463: *rp = 0;
1464: else {
1465: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1466: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5 noro 1467: }
1.7 noro 1468: if ( mr0 ) {
1469: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1470: } else
1471: *rp = 0;
1.5 noro 1472: }
1473: }
1474:
1475:
1.7 noro 1476: void homogenize_order(old,n,new)
1477: struct order_spec *old,*new;
1478: int n;
1.5 noro 1479: {
1.7 noro 1480: struct order_pair *l;
1481: int length,nv,row,i,j;
1482: int **newm,**oldm;
1.5 noro 1483:
1.7 noro 1484: switch ( old->id ) {
1485: case 0:
1486: switch ( old->ord.simple ) {
1487: case 0:
1488: new->id = 0; new->ord.simple = 0; break;
1489: case 1:
1490: l = (struct order_pair *)
1491: MALLOC_ATOMIC(2*sizeof(struct order_pair));
1492: l[0].length = n; l[0].order = 1;
1493: l[1].length = 1; l[1].order = 2;
1494: new->id = 1;
1495: new->ord.block.order_pair = l;
1496: new->ord.block.length = 2; new->nv = n+1;
1497: break;
1498: case 2:
1499: new->id = 0; new->ord.simple = 1; break;
1500: case 3: case 4: case 5:
1501: new->id = 0; new->ord.simple = old->ord.simple+3;
1502: dp_nelim = n-1; break;
1503: case 6: case 7: case 8: case 9:
1504: new->id = 0; new->ord.simple = old->ord.simple; break;
1505: default:
1506: error("homogenize_order : invalid input");
1507: }
1508: break;
1509: case 1:
1510: length = old->ord.block.length;
1511: l = (struct order_pair *)
1512: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
1513: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
1514: l[length].order = 2; l[length].length = 1;
1515: new->id = 1; new->nv = n+1;
1516: new->ord.block.order_pair = l;
1517: new->ord.block.length = length+1;
1518: break;
1519: case 2:
1520: nv = old->nv; row = old->ord.matrix.row;
1521: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
1522: for ( i = 0; i <= nv; i++ )
1523: newm[0][i] = 1;
1524: for ( i = 0; i < row; i++ ) {
1525: for ( j = 0; j < nv; j++ )
1526: newm[i+1][j] = oldm[i][j];
1527: newm[i+1][j] = 0;
1528: }
1529: new->id = 2; new->nv = nv+1;
1530: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1531: break;
1532: default:
1533: error("homogenize_order : invalid input");
1.5 noro 1534: }
1.7 noro 1535: }
1536:
1537: void qltozl(w,n,dvr)
1538: Q *w,*dvr;
1539: int n;
1540: {
1541: N nm,dn;
1542: N g,l1,l2,l3;
1543: Q c,d;
1544: int i;
1545: struct oVECT v;
1.5 noro 1546:
1547: for ( i = 0; i < n; i++ )
1.7 noro 1548: if ( w[i] && !INT(w[i]) )
1549: break;
1550: if ( i == n ) {
1551: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
1552: igcdv(&v,dvr); return;
1553: }
1554: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
1555: for ( i = 1; i < n; i++ ) {
1556: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
1557: gcdn(nm,NM(c),&g); nm = g;
1558: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 1559: }
1.7 noro 1560: if ( UNIN(dn) )
1561: NTOQ(nm,1,d);
1562: else
1563: NDTOQ(nm,dn,1,d);
1564: *dvr = d;
1565: }
1.5 noro 1566:
1.7 noro 1567: int comp_nm(a,b)
1568: Q *a,*b;
1569: {
1570: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
1571: }
1572:
1573: void sortbynm(w,n)
1574: Q *w;
1575: int n;
1576: {
1577: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
1578: }
1.5 noro 1579:
1580:
1.7 noro 1581: /*
1582: * simple operations
1583: *
1584: */
1.5 noro 1585:
1.7 noro 1586: int dp_redble(p1,p2)
1587: DP p1,p2;
1588: {
1589: int i,n;
1590: DL d1,d2;
1.5 noro 1591:
1.7 noro 1592: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1593: if ( d1->td < d2->td )
1594: return 0;
1595: else {
1596: for ( i = 0, n = p1->nv; i < n; i++ )
1597: if ( d1->d[i] < d2->d[i] )
1598: return 0;
1599: return 1;
1.5 noro 1600: }
1601: }
1602:
1.7 noro 1603: void dp_subd(p1,p2,rp)
1.5 noro 1604: DP p1,p2;
1605: DP *rp;
1606: {
1.7 noro 1607: int i,n;
1.5 noro 1608: DL d1,d2,d;
1609: MP m;
1.7 noro 1610: DP s;
1.5 noro 1611:
1612: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 noro 1613: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 1614: for ( i = 0; i < n; i++ )
1.7 noro 1615: d->d[i] = d1->d[i]-d2->d[i];
1616: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1617: MKDP(n,m,s); s->sugar = d->td;
1618: *rp = s;
1619: }
1620:
1621: void dltod(d,n,rp)
1622: DL d;
1623: int n;
1624: DP *rp;
1625: {
1626: MP m;
1627: DP s;
1628:
1629: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1630: MKDP(n,m,s); s->sugar = d->td;
1631: *rp = s;
1.5 noro 1632: }
1633:
1634: void dp_hm(p,rp)
1635: DP p;
1636: DP *rp;
1637: {
1638: MP m,mr;
1639:
1640: if ( !p )
1641: *rp = 0;
1642: else {
1643: m = BDY(p);
1644: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
1645: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1646: }
1647: }
1648:
1649: void dp_rest(p,rp)
1650: DP p,*rp;
1651: {
1652: MP m;
1653:
1654: m = BDY(p);
1655: if ( !NEXT(m) )
1656: *rp = 0;
1657: else {
1658: MKDP(p->nv,NEXT(m),*rp);
1659: if ( *rp )
1660: (*rp)->sugar = p->sugar;
1661: }
1662: }
1663:
1664: DL lcm_of_DL(nv,dl1,dl2,dl)
1665: int nv;
1666: DL dl1,dl2;
1667: register DL dl;
1668: {
1669: register int n, *d1, *d2, *d, td;
1670:
1671: if ( !dl ) NEWDL(dl,nv);
1672: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1673: for ( td = 0, n = nv; --n >= 0; d1++, d2++, d++ )
1674: td += (*d = *d1 > *d2 ? *d1 : *d2 );
1675: dl->td = td;
1676: return dl;
1677: }
1678:
1679: int dl_equal(nv,dl1,dl2)
1680: int nv;
1681: DL dl1, dl2;
1682: {
1683: register int *d1, *d2, n;
1684:
1685: if ( dl1->td != dl2->td ) return 0;
1686: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
1687: if ( *d1 != *d2 ) return 0;
1688: return 1;
1689: }
1690:
1691: int dp_nt(p)
1692: DP p;
1693: {
1694: int i;
1695: MP m;
1696:
1697: if ( !p )
1698: return 0;
1699: else {
1700: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
1701: return i;
1702: }
1703: }
1704:
1.15 noro 1705: int dp_homogeneous(p)
1706: DP p;
1707: {
1708: MP m;
1709: int d;
1710:
1711: if ( !p )
1712: return 1;
1713: else {
1714: m = BDY(p);
1715: d = m->dl->td;
1716: m = NEXT(m);
1717: for ( ; m; m = NEXT(m) ) {
1718: if ( m->dl->td != d )
1719: return 0;
1720: }
1721: return 1;
1722: }
1.16 noro 1723: }
1724:
1725: _print_mp(nv,m)
1726: int nv;
1727: MP m;
1728: {
1729: int i;
1730:
1.17 ! noro 1731: if ( !m )
1.16 noro 1732: return;
1733: for ( ; m; m = NEXT(m) ) {
1734: fprintf(stderr,"%d<",ITOS(C(m)));
1735: for ( i = 0; i < nv; i++ ) {
1736: fprintf(stderr,"%d",m->dl->d[i]);
1737: if ( i != nv-1 )
1738: fprintf(stderr," ");
1739: }
1740: fprintf(stderr,">",C(m));
1741: }
1742: fprintf(stderr,"\n");
1.15 noro 1743: }