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