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