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