![[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.13 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.12 2001/02/21 07:10:17 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);
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_mod(p1,p2,mod,rp)
573: DP p1,p2;
574: int mod;
575: DP *rp;
576: {
577: int i,n,td;
578: int *w;
579: DL d1,d2,d;
580: MP m;
581: DP t,s,u;
582:
583: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
584: w = (int *)ALLOCA(n*sizeof(int));
585: for ( i = 0, td = 0; i < n; i++ ) {
586: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
587: }
588: NEWDL(d,n); d->td = td - d1->td;
589: for ( i = 0; i < n; i++ )
590: d->d[i] = w[i] - d1->d[i];
591: NEWMP(m); m->dl = d; m->c = (P)BDY(p2)->c; NEXT(m) = 0;
592: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
593: NEWDL(d,n); d->td = td - d2->td;
594: for ( i = 0; i < n; i++ )
595: d->d[i] = w[i] - d2->d[i];
596: NEWMP(m); m->dl = d; m->c = (P)BDY(p1)->c; NEXT(m) = 0;
597: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
598: submd(CO,mod,t,u,rp);
599: }
600:
601: void _dp_sp_mod_dup(p1,p2,mod,rp)
602: DP p1,p2;
603: int mod;
604: DP *rp;
605: {
606: int i,n,td;
607: int *w;
608: DL d1,d2,d;
609: MP m;
610: DP t,s,u;
611:
612: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
613: w = (int *)ALLOCA(n*sizeof(int));
614: for ( i = 0, td = 0; i < n; i++ ) {
615: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
616: }
617: _NEWDL(d,n); d->td = td - d1->td;
618: for ( i = 0; i < n; i++ )
619: d->d[i] = w[i] - d1->d[i];
620: _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
621: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
622: _NEWDL(d,n); d->td = td - d2->td;
623: for ( i = 0; i < n; i++ )
624: d->d[i] = w[i] - d2->d[i];
625: _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
626: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
627: _addmd_destructive(mod,t,u,rp);
628: }
629:
630: void _dp_sp_mod(p1,p2,mod,rp)
631: DP p1,p2;
632: int mod;
633: DP *rp;
634: {
635: int i,n,td;
636: int *w;
637: DL d1,d2,d;
638: MP m;
639: DP t,s,u;
640:
641: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
642: w = (int *)ALLOCA(n*sizeof(int));
643: for ( i = 0, td = 0; i < n; i++ ) {
644: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
645: }
646: NEWDL(d,n); d->td = td - d1->td;
647: for ( i = 0; i < n; i++ )
648: d->d[i] = w[i] - d1->d[i];
649: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
650: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
651: NEWDL(d,n); d->td = td - d2->td;
652: for ( i = 0; i < n; i++ )
653: d->d[i] = w[i] - d2->d[i];
654: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
655: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
656: addmd_destructive(mod,t,u,rp);
657: }
658:
659: /*
660: * m-reduction
1.13 ! noro 661: * do content reduction over Z or Q(x,...)
! 662: * do nothing over finite fields
1.7 noro 663: *
664: */
665:
666: void dp_red(p0,p1,p2,head,rest,dnp,multp)
667: DP p0,p1,p2;
668: DP *head,*rest;
669: P *dnp;
670: DP *multp;
671: {
672: int i,n;
673: DL d1,d2,d;
674: MP m;
675: DP t,s,r,h;
676: Q c,c1,c2;
677: N gn,tn;
678: P g,a;
679:
680: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
681: NEWDL(d,n); d->td = d1->td - d2->td;
682: for ( i = 0; i < n; i++ )
683: d->d[i] = d1->d[i]-d2->d[i];
684: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
685: if ( dp_fcoeffs ) {
686: /* do nothing */
687: } else if ( INT(c1) && INT(c2) ) {
688: gcdn(NM(c1),NM(c2),&gn);
689: if ( !UNIN(gn) ) {
690: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
691: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
692: }
693: } else {
694: ezgcdpz(CO,(P)c1,(P)c2,&g);
695: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
696: }
697: NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
698: *multp = s;
699: muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
700: muldc(CO,p0,(P)c2,&h);
701: *head = h; *rest = r; *dnp = (P)c2;
702: }
703:
1.13 ! noro 704: /* m-reduction over a field */
! 705:
! 706: void dp_red_f(p1,p2,rest)
! 707: DP p1,p2;
! 708: DP *rest;
! 709: {
! 710: int i,n;
! 711: DL d1,d2,d;
! 712: MP m;
! 713: DP t,s,r,h;
! 714: Obj a,b;
! 715:
! 716: n = p1->nv;
! 717: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 718:
! 719: NEWDL(d,n); d->td = d1->td - d2->td;
! 720: for ( i = 0; i < n; i++ )
! 721: d->d[i] = d1->d[i]-d2->d[i];
! 722:
! 723: NEWMP(m); m->dl = d;
! 724: divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
! 725: C(m) = (P)b;
! 726: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
! 727:
! 728: muld(CO,s,p2,&t); addd(CO,p1,t,rest);
! 729: }
! 730:
1.7 noro 731: void dp_red_mod(p0,p1,p2,mod,head,rest,dnp)
732: DP p0,p1,p2;
733: int mod;
734: DP *head,*rest;
735: P *dnp;
736: {
737: int i,n;
738: DL d1,d2,d;
739: MP m;
740: DP t,s,r,h;
741: P c1,c2,g,u;
742:
743: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
744: NEWDL(d,n); d->td = d1->td - d2->td;
745: for ( i = 0; i < n; i++ )
746: d->d[i] = d1->d[i]-d2->d[i];
747: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
748: gcdprsmp(CO,mod,c1,c2,&g);
749: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
750: if ( NUM(c2) ) {
751: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
752: }
753: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1.11 noro 754: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7 noro 755: if ( NUM(c2) ) {
756: addmd(CO,mod,p1,t,&r); h = p0;
757: } else {
758: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
759: }
760: *head = h; *rest = r; *dnp = c2;
761: }
762:
1.10 noro 763: struct oEGT eg_red_mod;
764:
1.7 noro 765: void _dp_red_mod_destructive(p1,p2,mod,rp)
766: DP p1,p2;
767: int mod;
768: DP *rp;
769: {
770: int i,n;
771: DL d1,d2,d;
772: MP m;
773: DP t,s;
774: int c,c1;
1.10 noro 775: struct oEGT t0,t1;
1.7 noro 776:
777: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
778: _NEWDL(d,n); d->td = d1->td - d2->td;
779: for ( i = 0; i < n; i++ )
780: d->d[i] = d1->d[i]-d2->d[i];
781: c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
782: _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
783: _MKDP(n,m,s); s->sugar = d->td;
784: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.10 noro 785: /* get_eg(&t0); */
1.7 noro 786: _addmd_destructive(mod,p1,t,rp);
1.10 noro 787: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 788: }
789:
790: /*
791: * normal form computation
792: *
793: */
1.5 noro 794:
795: void dp_true_nf(b,g,ps,full,rp,dnp)
796: NODE b;
797: DP g;
798: DP *ps;
799: int full;
800: DP *rp;
801: P *dnp;
802: {
803: DP u,p,d,s,t,dmy;
804: NODE l;
805: MP m,mr;
806: int i,n;
807: int *wb;
808: int sugar,psugar;
809: P dn,tdn,tdn1;
810:
811: dn = (P)ONE;
812: if ( !g ) {
813: *rp = 0; *dnp = dn; return;
814: }
815: for ( n = 0, l = b; l; l = NEXT(l), n++ );
816: wb = (int *)ALLOCA(n*sizeof(int));
817: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
818: wb[i] = QTOS((Q)BDY(l));
819: sugar = g->sugar;
820: for ( d = 0; g; ) {
821: for ( u = 0, i = 0; i < n; i++ ) {
822: if ( dp_redble(g,p = ps[wb[i]]) ) {
823: dp_red(d,g,p,&t,&u,&tdn,&dmy);
824: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
825: sugar = MAX(sugar,psugar);
826: if ( !u ) {
827: if ( d )
828: d->sugar = sugar;
829: *rp = d; *dnp = dn; return;
830: } else {
831: d = t;
832: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
833: }
834: break;
835: }
836: }
837: if ( u )
838: g = u;
839: else if ( !full ) {
840: if ( g ) {
841: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
842: }
843: *rp = g; *dnp = dn; return;
844: } else {
845: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
846: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
847: addd(CO,d,t,&s); d = s;
848: dp_rest(g,&t); g = t;
849: }
850: }
851: if ( d )
852: d->sugar = sugar;
853: *rp = d; *dnp = dn;
854: }
855:
1.13 ! noro 856: /* nf computation over Z */
! 857:
! 858: void dp_nf_z(b,g,ps,full,multiple,rp)
1.5 noro 859: NODE b;
860: DP g;
861: DP *ps;
862: int full,multiple;
863: DP *rp;
864: {
865: DP u,p,d,s,t,dmy1;
866: P dmy;
867: NODE l;
868: MP m,mr;
869: int i,n;
870: int *wb;
871: int hmag;
872: int sugar,psugar;
873:
874: if ( !g ) {
875: *rp = 0; return;
876: }
877: for ( n = 0, l = b; l; l = NEXT(l), n++ );
878: wb = (int *)ALLOCA(n*sizeof(int));
879: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
880: wb[i] = QTOS((Q)BDY(l));
1.12 noro 881:
1.13 ! noro 882: hmag = multiple*HMAG(g);
1.5 noro 883: sugar = g->sugar;
1.12 noro 884:
1.5 noro 885: for ( d = 0; g; ) {
886: for ( u = 0, i = 0; i < n; i++ ) {
887: if ( dp_redble(g,p = ps[wb[i]]) ) {
888: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
889: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
890: sugar = MAX(sugar,psugar);
891: if ( !u ) {
892: if ( d )
893: d->sugar = sugar;
894: *rp = d; return;
895: }
896: d = t;
897: break;
898: }
899: }
900: if ( u ) {
901: g = u;
902: if ( d ) {
1.13 ! noro 903: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 904: dp_ptozp2(d,g,&t,&u); d = t; g = u;
905: hmag = multiple*HMAG(d);
906: }
907: } else {
1.13 ! noro 908: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 909: dp_ptozp(g,&t); g = t;
910: hmag = multiple*HMAG(g);
911: }
912: }
913: }
914: else if ( !full ) {
915: if ( g ) {
916: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
917: }
918: *rp = g; return;
919: } else {
920: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
921: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
922: addd(CO,d,t,&s); d = s;
923: dp_rest(g,&t); g = t;
924:
925: }
926: }
927: if ( d )
928: d->sugar = sugar;
929: *rp = d;
930: }
931:
1.13 ! noro 932: /* nf computation over a field */
! 933:
! 934: void dp_nf_f(b,g,ps,full,rp)
! 935: NODE b;
! 936: DP g;
! 937: DP *ps;
! 938: int full;
! 939: DP *rp;
! 940: {
! 941: DP u,p,d,s,t;
! 942: P dmy;
! 943: NODE l;
! 944: MP m,mr;
! 945: int i,n;
! 946: int *wb;
! 947: int sugar,psugar;
! 948:
! 949: if ( !g ) {
! 950: *rp = 0; return;
! 951: }
! 952: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 953: wb = (int *)ALLOCA(n*sizeof(int));
! 954: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 955: wb[i] = QTOS((Q)BDY(l));
! 956:
! 957: sugar = g->sugar;
! 958: for ( d = 0; g; ) {
! 959: for ( u = 0, i = 0; i < n; i++ ) {
! 960: if ( dp_redble(g,p = ps[wb[i]]) ) {
! 961: dp_red_f(g,p,&u);
! 962: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 963: sugar = MAX(sugar,psugar);
! 964: if ( !u ) {
! 965: if ( d )
! 966: d->sugar = sugar;
! 967: *rp = d; return;
! 968: }
! 969: break;
! 970: }
! 971: }
! 972: if ( u )
! 973: g = u;
! 974: else if ( !full ) {
! 975: if ( g ) {
! 976: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 977: }
! 978: *rp = g; return;
! 979: } else {
! 980: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
! 981: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
! 982: addd(CO,d,t,&s); d = s;
! 983: dp_rest(g,&t); g = t;
! 984: }
! 985: }
! 986: if ( d )
! 987: d->sugar = sugar;
! 988: *rp = d;
! 989: }
! 990:
! 991: /* nf computation over GF(mod) (only for internal use) */
! 992:
1.5 noro 993: void dp_nf_mod(b,g,ps,mod,full,rp)
994: NODE b;
995: DP g;
996: DP *ps;
997: int mod,full;
998: DP *rp;
999: {
1000: DP u,p,d,s,t;
1001: P dmy;
1002: NODE l;
1003: MP m,mr;
1004: int sugar,psugar;
1005:
1006: if ( !g ) {
1007: *rp = 0; return;
1008: }
1009: sugar = g->sugar;
1010: for ( d = 0; g; ) {
1011: for ( u = 0, l = b; l; l = NEXT(l) ) {
1012: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1013: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
1014: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1015: sugar = MAX(sugar,psugar);
1016: if ( !u ) {
1017: if ( d )
1018: d->sugar = sugar;
1019: *rp = d; return;
1020: }
1021: d = t;
1022: break;
1023: }
1024: }
1025: if ( u )
1026: g = u;
1027: else if ( !full ) {
1028: if ( g ) {
1029: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1030: }
1031: *rp = g; return;
1032: } else {
1033: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1034: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1035: addmd(CO,mod,d,t,&s); d = s;
1036: dp_rest(g,&t); g = t;
1037: }
1038: }
1039: if ( d )
1040: d->sugar = sugar;
1041: *rp = d;
1042: }
1043:
1044: void dp_true_nf_mod(b,g,ps,mod,full,rp,dnp)
1045: NODE b;
1046: DP g;
1047: DP *ps;
1048: int mod,full;
1049: DP *rp;
1050: P *dnp;
1051: {
1052: DP u,p,d,s,t;
1053: NODE l;
1054: MP m,mr;
1055: int i,n;
1056: int *wb;
1057: int sugar,psugar;
1058: P dn,tdn,tdn1;
1059:
1060: dn = (P)ONEM;
1061: if ( !g ) {
1062: *rp = 0; *dnp = dn; return;
1063: }
1064: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1065: wb = (int *)ALLOCA(n*sizeof(int));
1066: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1067: wb[i] = QTOS((Q)BDY(l));
1068: sugar = g->sugar;
1069: for ( d = 0; g; ) {
1070: for ( u = 0, i = 0; i < n; i++ ) {
1071: if ( dp_redble(g,p = ps[wb[i]]) ) {
1072: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
1073: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1074: sugar = MAX(sugar,psugar);
1075: if ( !u ) {
1076: if ( d )
1077: d->sugar = sugar;
1078: *rp = d; *dnp = dn; return;
1079: } else {
1080: d = t;
1081: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1082: }
1083: break;
1084: }
1085: }
1086: if ( u )
1087: g = u;
1088: else if ( !full ) {
1089: if ( g ) {
1090: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1091: }
1092: *rp = g; *dnp = dn; return;
1093: } else {
1094: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1095: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1096: addmd(CO,mod,d,t,&s); d = s;
1097: dp_rest(g,&t); g = t;
1098: }
1099: }
1100: if ( d )
1101: d->sugar = sugar;
1102: *rp = d; *dnp = dn;
1103: }
1104:
1.7 noro 1105: void _dp_nf_mod_destructive(b,g,ps,mod,full,rp)
1106: NODE b;
1107: DP g;
1108: DP *ps;
1109: int mod,full;
1110: DP *rp;
1.5 noro 1111: {
1.7 noro 1112: DP u,p,d,s,t;
1113: NODE l;
1114: MP m,mr,mrd;
1115: int sugar,psugar,n,h_reducible,i;
1.5 noro 1116:
1.7 noro 1117: if ( !g ) {
1118: *rp = 0; return;
1.5 noro 1119: }
1.7 noro 1120: sugar = g->sugar;
1121: n = g->nv;
1122: for ( d = 0; g; ) {
1123: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
1124: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1125: h_reducible = 1;
1126: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1127: _dp_red_mod_destructive(g,p,mod,&u); g = u;
1128: sugar = MAX(sugar,psugar);
1129: if ( !g ) {
1130: if ( d )
1131: d->sugar = sugar;
1132: _dptodp(d,rp); _free_dp(d); return;
1133: }
1134: break;
1135: }
1136: }
1137: if ( !h_reducible ) {
1138: /* head term is not reducible */
1139: if ( !full ) {
1140: if ( g )
1141: g->sugar = sugar;
1142: _dptodp(g,rp); _free_dp(g); return;
1143: } else {
1144: m = BDY(g);
1145: if ( NEXT(m) ) {
1146: BDY(g) = NEXT(m); NEXT(m) = 0;
1147: } else {
1148: _FREEDP(g); g = 0;
1149: }
1150: if ( d ) {
1151: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1152: NEXT(mrd) = m;
1153: } else {
1154: _MKDP(n,m,d);
1155: }
1156: }
1157: }
1.5 noro 1158: }
1.7 noro 1159: if ( d )
1160: d->sugar = sugar;
1161: _dptodp(d,rp); _free_dp(d);
1.5 noro 1162: }
1.13 ! noro 1163:
! 1164: /* reduction by linear base over a field */
! 1165:
! 1166: void dp_lnf_f(p1,p2,g,r1p,r2p)
! 1167: DP p1,p2;
! 1168: NODE g;
! 1169: DP *r1p,*r2p;
! 1170: {
! 1171: DP r1,r2,b1,b2,t,s;
! 1172: Obj c,c1,c2;
! 1173: NODE l,b;
! 1174: int n;
! 1175:
! 1176: if ( !p1 ) {
! 1177: *r1p = p1; *r2p = p2; return;
! 1178: }
! 1179: n = p1->nv;
! 1180: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
! 1181: if ( !r1 ) {
! 1182: *r1p = r1; *r2p = r2; return;
! 1183: }
! 1184: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
! 1185: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
! 1186: b2 = (DP)BDY(NEXT(b));
! 1187: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
! 1188: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
! 1189: muldc(CO,b1,(P)c,&t); addd(CO,r1,t,&s); r1 = s;
! 1190: muldc(CO,b2,(P)c,&t); addd(CO,r2,t,&s); r2 = s;
! 1191: }
! 1192: }
! 1193: *r1p = r1; *r2p = r2;
! 1194: }
! 1195:
! 1196: /* reduction by linear base over GF(mod) */
1.5 noro 1197:
1.7 noro 1198: void dp_lnf_mod(p1,p2,g,mod,r1p,r2p)
1199: DP p1,p2;
1200: NODE g;
1201: int mod;
1202: DP *r1p,*r2p;
1.5 noro 1203: {
1.7 noro 1204: DP r1,r2,b1,b2,t,s;
1205: P c;
1206: MQ c1,c2;
1207: NODE l,b;
1208: int n;
1209:
1210: if ( !p1 ) {
1211: *r1p = p1; *r2p = p2; return;
1212: }
1213: n = p1->nv;
1214: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1215: if ( !r1 ) {
1216: *r1p = r1; *r2p = r2; return;
1217: }
1218: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1219: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1220: b2 = (DP)BDY(NEXT(b));
1221: invmq(mod,(MQ)BDY(b1)->c,&c1);
1222: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
1223: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
1224: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
1225: }
1226: }
1227: *r1p = r1; *r2p = r2;
1.5 noro 1228: }
1229:
1.7 noro 1230: void dp_nf_tab_mod(p,tab,mod,rp)
1231: DP p;
1232: LIST *tab;
1233: int mod;
1234: DP *rp;
1.5 noro 1235: {
1.7 noro 1236: DP s,t,u;
1237: MP m;
1238: DL h;
1239: int i,n;
1240:
1241: if ( !p ) {
1242: *rp = p; return;
1243: }
1244: n = p->nv;
1245: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
1246: h = m->dl;
1247: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
1248: i++;
1249: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
1250: addmd(CO,mod,s,t,&u); s = u;
1251: }
1252: *rp = s;
1.5 noro 1253: }
1254:
1.7 noro 1255: /*
1256: * setting flags
1257: *
1258: */
1259:
1260: int create_order_spec(obj,spec)
1261: Obj obj;
1262: struct order_spec *spec;
1.5 noro 1263: {
1.7 noro 1264: int i,j,n,s,row,col;
1265: struct order_pair *l;
1266: NODE node,t,tn;
1267: MAT m;
1268: pointer **b;
1269: int **w;
1.5 noro 1270:
1.7 noro 1271: if ( !obj || NUM(obj) ) {
1272: spec->id = 0; spec->obj = obj;
1273: spec->ord.simple = QTOS((Q)obj);
1274: return 1;
1275: } else if ( OID(obj) == O_LIST ) {
1276: node = BDY((LIST)obj);
1277: for ( n = 0, t = node; t; t = NEXT(t), n++ );
1278: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
1279: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
1280: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
1281: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
1282: s += l[i].length;
1283: }
1284: spec->id = 1; spec->obj = obj;
1285: spec->ord.block.order_pair = l;
1286: spec->ord.block.length = n; spec->nv = s;
1287: return 1;
1288: } else if ( OID(obj) == O_MAT ) {
1289: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
1290: w = almat(row,col);
1291: for ( i = 0; i < row; i++ )
1292: for ( j = 0; j < col; j++ )
1293: w[i][j] = QTOS((Q)b[i][j]);
1294: spec->id = 2; spec->obj = obj;
1295: spec->nv = col; spec->ord.matrix.row = row;
1296: spec->ord.matrix.matrix = w;
1297: return 1;
1298: } else
1.5 noro 1299: return 0;
1300: }
1301:
1.7 noro 1302: /*
1303: * converters
1304: *
1305: */
1306:
1307: void dp_homo(p,rp)
1308: DP p;
1309: DP *rp;
1.5 noro 1310: {
1.7 noro 1311: MP m,mr,mr0;
1312: int i,n,nv,td;
1313: DL dl,dlh;
1.5 noro 1314:
1.7 noro 1315: if ( !p )
1316: *rp = 0;
1317: else {
1318: n = p->nv; nv = n + 1;
1319: m = BDY(p); td = sugard(m);
1320: for ( mr0 = 0; m; m = NEXT(m) ) {
1321: NEXTMP(mr0,mr); mr->c = m->c;
1322: dl = m->dl;
1323: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1324: dlh->td = td;
1325: for ( i = 0; i < n; i++ )
1326: dlh->d[i] = dl->d[i];
1327: dlh->d[n] = td - dl->td;
1328: }
1329: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 1330: }
1331: }
1332:
1.7 noro 1333: void dp_dehomo(p,rp)
1334: DP p;
1.5 noro 1335: DP *rp;
1336: {
1.7 noro 1337: MP m,mr,mr0;
1338: int i,n,nv;
1339: DL dl,dlh;
1.5 noro 1340:
1.7 noro 1341: if ( !p )
1342: *rp = 0;
1343: else {
1344: n = p->nv; nv = n - 1;
1345: m = BDY(p);
1346: for ( mr0 = 0; m; m = NEXT(m) ) {
1347: NEXTMP(mr0,mr); mr->c = m->c;
1348: dlh = m->dl;
1349: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1350: dl->td = dlh->td - dlh->d[nv];
1351: for ( i = 0; i < nv; i++ )
1352: dl->d[i] = dlh->d[i];
1353: }
1354: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1355: }
1.5 noro 1356: }
1357:
1.7 noro 1358: void dp_mod(p,mod,subst,rp)
1359: DP p;
1360: int mod;
1361: NODE subst;
1.5 noro 1362: DP *rp;
1363: {
1.7 noro 1364: MP m,mr,mr0;
1365: P t,s,s1;
1366: V v;
1367: NODE tn;
1.5 noro 1368:
1.7 noro 1369: if ( !p )
1370: *rp = 0;
1371: else {
1372: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1373: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
1374: v = VR((P)BDY(tn)); tn = NEXT(tn);
1375: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
1376: }
1377: ptomp(mod,s,&t);
1378: if ( t ) {
1379: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
1380: }
1381: }
1382: if ( mr0 ) {
1383: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1384: } else
1385: *rp = 0;
1386: }
1.5 noro 1387: }
1388:
1.7 noro 1389: void dp_rat(p,rp)
1390: DP p;
1391: DP *rp;
1.5 noro 1392: {
1.7 noro 1393: MP m,mr,mr0;
1.5 noro 1394:
1.7 noro 1395: if ( !p )
1396: *rp = 0;
1397: else {
1398: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1399: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5 noro 1400: }
1.7 noro 1401: if ( mr0 ) {
1402: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1403: } else
1404: *rp = 0;
1.5 noro 1405: }
1406: }
1407:
1408:
1.7 noro 1409: void homogenize_order(old,n,new)
1410: struct order_spec *old,*new;
1411: int n;
1.5 noro 1412: {
1.7 noro 1413: struct order_pair *l;
1414: int length,nv,row,i,j;
1415: int **newm,**oldm;
1.5 noro 1416:
1.7 noro 1417: switch ( old->id ) {
1418: case 0:
1419: switch ( old->ord.simple ) {
1420: case 0:
1421: new->id = 0; new->ord.simple = 0; break;
1422: case 1:
1423: l = (struct order_pair *)
1424: MALLOC_ATOMIC(2*sizeof(struct order_pair));
1425: l[0].length = n; l[0].order = 1;
1426: l[1].length = 1; l[1].order = 2;
1427: new->id = 1;
1428: new->ord.block.order_pair = l;
1429: new->ord.block.length = 2; new->nv = n+1;
1430: break;
1431: case 2:
1432: new->id = 0; new->ord.simple = 1; break;
1433: case 3: case 4: case 5:
1434: new->id = 0; new->ord.simple = old->ord.simple+3;
1435: dp_nelim = n-1; break;
1436: case 6: case 7: case 8: case 9:
1437: new->id = 0; new->ord.simple = old->ord.simple; break;
1438: default:
1439: error("homogenize_order : invalid input");
1440: }
1441: break;
1442: case 1:
1443: length = old->ord.block.length;
1444: l = (struct order_pair *)
1445: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
1446: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
1447: l[length].order = 2; l[length].length = 1;
1448: new->id = 1; new->nv = n+1;
1449: new->ord.block.order_pair = l;
1450: new->ord.block.length = length+1;
1451: break;
1452: case 2:
1453: nv = old->nv; row = old->ord.matrix.row;
1454: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
1455: for ( i = 0; i <= nv; i++ )
1456: newm[0][i] = 1;
1457: for ( i = 0; i < row; i++ ) {
1458: for ( j = 0; j < nv; j++ )
1459: newm[i+1][j] = oldm[i][j];
1460: newm[i+1][j] = 0;
1461: }
1462: new->id = 2; new->nv = nv+1;
1463: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1464: break;
1465: default:
1466: error("homogenize_order : invalid input");
1.5 noro 1467: }
1.7 noro 1468: }
1469:
1470: void qltozl(w,n,dvr)
1471: Q *w,*dvr;
1472: int n;
1473: {
1474: N nm,dn;
1475: N g,l1,l2,l3;
1476: Q c,d;
1477: int i;
1478: struct oVECT v;
1.5 noro 1479:
1480: for ( i = 0; i < n; i++ )
1.7 noro 1481: if ( w[i] && !INT(w[i]) )
1482: break;
1483: if ( i == n ) {
1484: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
1485: igcdv(&v,dvr); return;
1486: }
1487: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
1488: for ( i = 1; i < n; i++ ) {
1489: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
1490: gcdn(nm,NM(c),&g); nm = g;
1491: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 1492: }
1.7 noro 1493: if ( UNIN(dn) )
1494: NTOQ(nm,1,d);
1495: else
1496: NDTOQ(nm,dn,1,d);
1497: *dvr = d;
1498: }
1.5 noro 1499:
1.7 noro 1500: int comp_nm(a,b)
1501: Q *a,*b;
1502: {
1503: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
1504: }
1505:
1506: void sortbynm(w,n)
1507: Q *w;
1508: int n;
1509: {
1510: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
1511: }
1.5 noro 1512:
1513:
1.7 noro 1514: /*
1515: * simple operations
1516: *
1517: */
1.5 noro 1518:
1.7 noro 1519: int dp_redble(p1,p2)
1520: DP p1,p2;
1521: {
1522: int i,n;
1523: DL d1,d2;
1.5 noro 1524:
1.7 noro 1525: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1526: if ( d1->td < d2->td )
1527: return 0;
1528: else {
1529: for ( i = 0, n = p1->nv; i < n; i++ )
1530: if ( d1->d[i] < d2->d[i] )
1531: return 0;
1532: return 1;
1.5 noro 1533: }
1534: }
1535:
1.7 noro 1536: void dp_subd(p1,p2,rp)
1.5 noro 1537: DP p1,p2;
1538: DP *rp;
1539: {
1.7 noro 1540: int i,n;
1.5 noro 1541: DL d1,d2,d;
1542: MP m;
1.7 noro 1543: DP s;
1.5 noro 1544:
1545: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 noro 1546: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 1547: for ( i = 0; i < n; i++ )
1.7 noro 1548: d->d[i] = d1->d[i]-d2->d[i];
1549: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1550: MKDP(n,m,s); s->sugar = d->td;
1551: *rp = s;
1552: }
1553:
1554: void dltod(d,n,rp)
1555: DL d;
1556: int n;
1557: DP *rp;
1558: {
1559: MP m;
1560: DP s;
1561:
1562: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1563: MKDP(n,m,s); s->sugar = d->td;
1564: *rp = s;
1.5 noro 1565: }
1566:
1567: void dp_hm(p,rp)
1568: DP p;
1569: DP *rp;
1570: {
1571: MP m,mr;
1572:
1573: if ( !p )
1574: *rp = 0;
1575: else {
1576: m = BDY(p);
1577: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
1578: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1579: }
1580: }
1581:
1582: void dp_rest(p,rp)
1583: DP p,*rp;
1584: {
1585: MP m;
1586:
1587: m = BDY(p);
1588: if ( !NEXT(m) )
1589: *rp = 0;
1590: else {
1591: MKDP(p->nv,NEXT(m),*rp);
1592: if ( *rp )
1593: (*rp)->sugar = p->sugar;
1594: }
1595: }
1596:
1597: DL lcm_of_DL(nv,dl1,dl2,dl)
1598: int nv;
1599: DL dl1,dl2;
1600: register DL dl;
1601: {
1602: register int n, *d1, *d2, *d, td;
1603:
1604: if ( !dl ) NEWDL(dl,nv);
1605: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1606: for ( td = 0, n = nv; --n >= 0; d1++, d2++, d++ )
1607: td += (*d = *d1 > *d2 ? *d1 : *d2 );
1608: dl->td = td;
1609: return dl;
1610: }
1611:
1612: int dl_equal(nv,dl1,dl2)
1613: int nv;
1614: DL dl1, dl2;
1615: {
1616: register int *d1, *d2, n;
1617:
1618: if ( dl1->td != dl2->td ) return 0;
1619: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
1620: if ( *d1 != *d2 ) return 0;
1621: return 1;
1622: }
1623:
1624: int dp_nt(p)
1625: DP p;
1626: {
1627: int i;
1628: MP m;
1629:
1630: if ( !p )
1631: return 0;
1632: else {
1633: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
1634: return i;
1635: }
1636: }
1637: