![[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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.6 2000/12/05 08:29:43 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:
260: void dp_ptozp_d(dist,ndist,p,rp)
261: NODE dist;
262: int ndist;
263: DP p,*rp;
264: {
265: int i,j,k,l,n,nsep;
266: MP m;
267: NODE tn,n0,n1,n2,n3;
268: struct oVECT v;
269: VECT c,cs;
270: VECT qi,ri;
271: LIST *qr;
272: int s,id;
273: Obj dmy;
274: Q d0,d1,gcd,a,u,u1;
275: Q *q,*r;
276: STRING iqr_v;
277: pointer *b;
278: N qn,gn;
279: double get_rtime();
280: int blen;
281: double t0;
282: double t_e,t_d,t_d1,t_c;
283:
284: if ( !p )
285: *rp = 0;
286: else {
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:
349: void dp_ptozp2_d(dist,ndist,p0,p1,hp,rp)
350: NODE dist;
351: int ndist;
352: DP p0,p1;
353: DP *hp,*rp;
354: {
355: DP t,s,h,r;
356: MP m,mr,mr0,m0;
357:
358: addd(CO,p0,p1,&t); dp_ptozp_d(dist,ndist,t,&s);
359: if ( !p0 ) {
360: h = 0; r = s;
361: } else if ( !p1 ) {
362: h = s; r = 0;
363: } else {
364: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
365: m = NEXT(m), m0 = NEXT(m0) ) {
366: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
367: }
368: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
369: }
370: if ( h )
371: h->sugar = p0->sugar;
372: if ( r )
373: r->sugar = p1->sugar;
374: *hp = h; *rp = r;
1.5 noro 375: }
376:
1.7 ! noro 377: void dp_prim(p,rp)
! 378: DP p,*rp;
1.5 noro 379: {
1.7 ! noro 380: P t,g;
! 381: DP p1;
! 382: MP m,mr,mr0;
! 383: int i,n;
! 384: P *w;
! 385: Q *c;
! 386: Q dvr;
1.5 noro 387:
1.7 ! noro 388: if ( !p )
! 389: *rp = 0;
! 390: else if ( dp_fcoeffs )
! 391: *rp = p;
! 392: else if ( NoGCD )
! 393: dp_ptozp(p,rp);
! 394: else {
! 395: dp_ptozp(p,&p1); p = p1;
! 396: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
! 397: if ( n == 1 ) {
! 398: m = BDY(p);
! 399: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
! 400: MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
! 401: return;
! 402: }
! 403: w = (P *)ALLOCA(n*sizeof(P));
! 404: c = (Q *)ALLOCA(n*sizeof(Q));
! 405: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
! 406: if ( NUM(m->c) ) {
! 407: c[i] = (Q)m->c; w[i] = (P)ONE;
! 408: } else
! 409: ptozp(m->c,1,&c[i],&w[i]);
! 410: qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
! 411: if ( NUM(g) )
! 412: *rp = p;
! 413: else {
! 414: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 415: NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
! 416: }
! 417: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 418: }
1.7 ! noro 419: }
1.5 noro 420: }
421:
422: void heu_nezgcdnpz(vl,pl,m,pr)
423: VL vl;
424: P *pl,*pr;
425: int m;
426: {
427: int i,r;
428: P gcd,t,s1,s2,u;
429: Q rq;
430:
431: while ( 1 ) {
432: for ( i = 0, s1 = 0; i < m; i++ ) {
433: r = random(); UTOQ(r,rq);
434: mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
435: }
436: for ( i = 0, s2 = 0; i < m; i++ ) {
437: r = random(); UTOQ(r,rq);
438: mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
439: }
440: ezgcdp(vl,s1,s2,&gcd);
441: for ( i = 0; i < m; i++ ) {
442: if ( !divtpz(vl,pl[i],gcd,&t) )
443: break;
444: }
445: if ( i == m )
446: break;
447: }
448: *pr = gcd;
449: }
450:
451: void dp_prim_mod(p,mod,rp)
452: int mod;
453: DP p,*rp;
454: {
455: P t,g;
456: MP m,mr,mr0;
457:
458: if ( !p )
459: *rp = 0;
460: else if ( NoGCD )
461: *rp = p;
462: else {
463: for ( m = BDY(p), g = m->c, m = NEXT(m); m; m = NEXT(m) ) {
464: gcdprsmp(CO,mod,g,m->c,&t); g = t;
465: }
466: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
467: NEXTMP(mr0,mr); divsmp(CO,mod,m->c,g,&mr->c); mr->dl = m->dl;
468: }
469: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
470: }
471: }
472:
1.7 ! noro 473: void dp_cont(p,rp)
1.5 noro 474: DP p;
1.7 ! noro 475: Q *rp;
1.5 noro 476: {
1.7 ! noro 477: VECT v;
1.5 noro 478:
1.7 ! noro 479: dp_dtov(p,&v); igcdv(v,rp);
1.5 noro 480: }
481:
1.7 ! noro 482: void dp_dtov(dp,rp)
! 483: DP dp;
! 484: VECT *rp;
1.5 noro 485: {
1.7 ! noro 486: MP m,t;
! 487: int i,n;
! 488: VECT v;
! 489: pointer *p;
1.5 noro 490:
1.7 ! noro 491: m = BDY(dp);
! 492: for ( t = m, n = 0; t; t = NEXT(t), n++ );
! 493: MKVECT(v,n);
! 494: for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
! 495: p[i] = (pointer)(t->c);
! 496: *rp = v;
1.5 noro 497: }
498:
1.7 ! noro 499: /*
! 500: * s-poly computation
! 501: *
! 502: */
1.5 noro 503:
1.7 ! noro 504: void dp_sp(p1,p2,rp)
! 505: DP p1,p2;
1.5 noro 506: DP *rp;
507: {
1.7 ! noro 508: int i,n,td;
! 509: int *w;
! 510: DL d1,d2,d;
! 511: MP m;
! 512: DP t,s1,s2,u;
! 513: Q c,c1,c2;
! 514: N gn,tn;
1.5 noro 515:
1.7 ! noro 516: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 517: w = (int *)ALLOCA(n*sizeof(int));
! 518: for ( i = 0, td = 0; i < n; i++ ) {
! 519: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
1.5 noro 520: }
1.7 ! noro 521:
! 522: NEWDL(d,n); d->td = td - d1->td;
! 523: for ( i = 0; i < n; i++ )
! 524: d->d[i] = w[i] - d1->d[i];
! 525: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 526: if ( INT(c1) && INT(c2) ) {
! 527: gcdn(NM(c1),NM(c2),&gn);
! 528: if ( !UNIN(gn) ) {
! 529: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 530: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1.5 noro 531: }
532: }
1.7 ! noro 533:
! 534: NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
! 535: MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
! 536:
! 537: NEWDL(d,n); d->td = td - d2->td;
! 538: for ( i = 0; i < n; i++ )
! 539: d->d[i] = w[i] - d2->d[i];
! 540: NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0;
! 541: MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
! 542:
! 543: subd(CO,t,u,rp);
! 544: if ( GenTrace ) {
! 545: LIST hist;
! 546: NODE node;
! 547:
! 548: node = mknode(4,ONE,0,s1,ONE);
! 549: MKLIST(hist,node);
! 550: MKNODE(TraceList,hist,0);
! 551:
! 552: node = mknode(4,ONE,0,0,ONE);
! 553: chsgnd(s2,(DP *)&ARG2(node));
! 554: MKLIST(hist,node);
! 555: MKNODE(node,hist,TraceList); TraceList = node;
! 556: }
! 557: }
! 558:
! 559: void dp_sp_mod(p1,p2,mod,rp)
! 560: DP p1,p2;
! 561: int mod;
! 562: DP *rp;
! 563: {
! 564: int i,n,td;
! 565: int *w;
! 566: DL d1,d2,d;
! 567: MP m;
! 568: DP t,s,u;
! 569:
! 570: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 571: w = (int *)ALLOCA(n*sizeof(int));
! 572: for ( i = 0, td = 0; i < n; i++ ) {
! 573: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
! 574: }
! 575: NEWDL(d,n); d->td = td - d1->td;
! 576: for ( i = 0; i < n; i++ )
! 577: d->d[i] = w[i] - d1->d[i];
! 578: NEWMP(m); m->dl = d; m->c = (P)BDY(p2)->c; NEXT(m) = 0;
! 579: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
! 580: NEWDL(d,n); d->td = td - d2->td;
! 581: for ( i = 0; i < n; i++ )
! 582: d->d[i] = w[i] - d2->d[i];
! 583: NEWMP(m); m->dl = d; m->c = (P)BDY(p1)->c; NEXT(m) = 0;
! 584: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
! 585: submd(CO,mod,t,u,rp);
! 586: }
! 587:
! 588: void _dp_sp_mod_dup(p1,p2,mod,rp)
! 589: DP p1,p2;
! 590: int mod;
! 591: DP *rp;
! 592: {
! 593: int i,n,td;
! 594: int *w;
! 595: DL d1,d2,d;
! 596: MP m;
! 597: DP t,s,u;
! 598:
! 599: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 600: w = (int *)ALLOCA(n*sizeof(int));
! 601: for ( i = 0, td = 0; i < n; i++ ) {
! 602: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
! 603: }
! 604: _NEWDL(d,n); d->td = td - d1->td;
! 605: for ( i = 0; i < n; i++ )
! 606: d->d[i] = w[i] - d1->d[i];
! 607: _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
! 608: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
! 609: _NEWDL(d,n); d->td = td - d2->td;
! 610: for ( i = 0; i < n; i++ )
! 611: d->d[i] = w[i] - d2->d[i];
! 612: _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
! 613: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
! 614: _addmd_destructive(mod,t,u,rp);
! 615: }
! 616:
! 617: void _dp_sp_mod(p1,p2,mod,rp)
! 618: DP p1,p2;
! 619: int mod;
! 620: DP *rp;
! 621: {
! 622: int i,n,td;
! 623: int *w;
! 624: DL d1,d2,d;
! 625: MP m;
! 626: DP t,s,u;
! 627:
! 628: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 629: w = (int *)ALLOCA(n*sizeof(int));
! 630: for ( i = 0, td = 0; i < n; i++ ) {
! 631: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
! 632: }
! 633: NEWDL(d,n); d->td = td - d1->td;
! 634: for ( i = 0; i < n; i++ )
! 635: d->d[i] = w[i] - d1->d[i];
! 636: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
! 637: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
! 638: NEWDL(d,n); d->td = td - d2->td;
! 639: for ( i = 0; i < n; i++ )
! 640: d->d[i] = w[i] - d2->d[i];
! 641: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
! 642: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
! 643: addmd_destructive(mod,t,u,rp);
! 644: }
! 645:
! 646: /*
! 647: * m-reduction
! 648: *
! 649: */
! 650:
! 651: void dp_red(p0,p1,p2,head,rest,dnp,multp)
! 652: DP p0,p1,p2;
! 653: DP *head,*rest;
! 654: P *dnp;
! 655: DP *multp;
! 656: {
! 657: int i,n;
! 658: DL d1,d2,d;
! 659: MP m;
! 660: DP t,s,r,h;
! 661: Q c,c1,c2;
! 662: N gn,tn;
! 663: P g,a;
! 664:
! 665: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 666: NEWDL(d,n); d->td = d1->td - d2->td;
! 667: for ( i = 0; i < n; i++ )
! 668: d->d[i] = d1->d[i]-d2->d[i];
! 669: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 670: if ( dp_fcoeffs ) {
! 671: /* do nothing */
! 672: } else if ( INT(c1) && INT(c2) ) {
! 673: gcdn(NM(c1),NM(c2),&gn);
! 674: if ( !UNIN(gn) ) {
! 675: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 676: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 677: }
! 678: } else {
! 679: ezgcdpz(CO,(P)c1,(P)c2,&g);
! 680: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 681: }
! 682: NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
! 683: *multp = s;
! 684: muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
! 685: muldc(CO,p0,(P)c2,&h);
! 686: *head = h; *rest = r; *dnp = (P)c2;
! 687: }
! 688:
! 689: void dp_red_mod(p0,p1,p2,mod,head,rest,dnp)
! 690: DP p0,p1,p2;
! 691: int mod;
! 692: DP *head,*rest;
! 693: P *dnp;
! 694: {
! 695: int i,n;
! 696: DL d1,d2,d;
! 697: MP m;
! 698: DP t,s,r,h;
! 699: P c1,c2,g,u;
! 700:
! 701: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 702: NEWDL(d,n); d->td = d1->td - d2->td;
! 703: for ( i = 0; i < n; i++ )
! 704: d->d[i] = d1->d[i]-d2->d[i];
! 705: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
! 706: gcdprsmp(CO,mod,c1,c2,&g);
! 707: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
! 708: if ( NUM(c2) ) {
! 709: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
! 710: }
! 711: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
! 712: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&t);
! 713: if ( NUM(c2) ) {
! 714: addmd(CO,mod,p1,t,&r); h = p0;
! 715: } else {
! 716: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
! 717: }
! 718: *head = h; *rest = r; *dnp = c2;
! 719: }
! 720:
! 721: void _dp_red_mod_destructive(p1,p2,mod,rp)
! 722: DP p1,p2;
! 723: int mod;
! 724: DP *rp;
! 725: {
! 726: int i,n;
! 727: DL d1,d2,d;
! 728: MP m;
! 729: DP t,s;
! 730: int c,c1;
! 731:
! 732: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 733: _NEWDL(d,n); d->td = d1->td - d2->td;
! 734: for ( i = 0; i < n; i++ )
! 735: d->d[i] = d1->d[i]-d2->d[i];
! 736: c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
! 737: _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
! 738: _MKDP(n,m,s); s->sugar = d->td;
! 739: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
! 740: _addmd_destructive(mod,p1,t,rp);
! 741: }
! 742:
! 743: /*
! 744: * normal form computation
! 745: *
! 746: */
1.5 noro 747:
748: void dp_true_nf(b,g,ps,full,rp,dnp)
749: NODE b;
750: DP g;
751: DP *ps;
752: int full;
753: DP *rp;
754: P *dnp;
755: {
756: DP u,p,d,s,t,dmy;
757: NODE l;
758: MP m,mr;
759: int i,n;
760: int *wb;
761: int sugar,psugar;
762: P dn,tdn,tdn1;
763:
764: dn = (P)ONE;
765: if ( !g ) {
766: *rp = 0; *dnp = dn; return;
767: }
768: for ( n = 0, l = b; l; l = NEXT(l), n++ );
769: wb = (int *)ALLOCA(n*sizeof(int));
770: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
771: wb[i] = QTOS((Q)BDY(l));
772: sugar = g->sugar;
773: for ( d = 0; g; ) {
774: for ( u = 0, i = 0; i < n; i++ ) {
775: if ( dp_redble(g,p = ps[wb[i]]) ) {
776: dp_red(d,g,p,&t,&u,&tdn,&dmy);
777: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
778: sugar = MAX(sugar,psugar);
779: if ( !u ) {
780: if ( d )
781: d->sugar = sugar;
782: *rp = d; *dnp = dn; return;
783: } else {
784: d = t;
785: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
786: }
787: break;
788: }
789: }
790: if ( u )
791: g = u;
792: else if ( !full ) {
793: if ( g ) {
794: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
795: }
796: *rp = g; *dnp = dn; return;
797: } else {
798: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
799: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
800: addd(CO,d,t,&s); d = s;
801: dp_rest(g,&t); g = t;
802: }
803: }
804: if ( d )
805: d->sugar = sugar;
806: *rp = d; *dnp = dn;
807: }
808:
809: void dp_nf_ptozp(b,g,ps,full,multiple,rp)
810: NODE b;
811: DP g;
812: DP *ps;
813: int full,multiple;
814: DP *rp;
815: {
816: DP u,p,d,s,t,dmy1;
817: P dmy;
818: NODE l;
819: MP m,mr;
820: int i,n;
821: int *wb;
822: int hmag;
823: int sugar,psugar;
824:
825: if ( !g ) {
826: *rp = 0; return;
827: }
828: for ( n = 0, l = b; l; l = NEXT(l), n++ );
829: wb = (int *)ALLOCA(n*sizeof(int));
830: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
831: wb[i] = QTOS((Q)BDY(l));
832: hmag = multiple*HMAG(g);
833: sugar = g->sugar;
834: for ( d = 0; g; ) {
835: for ( u = 0, i = 0; i < n; i++ ) {
836: if ( dp_redble(g,p = ps[wb[i]]) ) {
837: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
838: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
839: sugar = MAX(sugar,psugar);
840: if ( !u ) {
841: if ( d )
842: d->sugar = sugar;
843: *rp = d; return;
844: }
845: d = t;
846: break;
847: }
848: }
849: if ( u ) {
850: g = u;
851: if ( d ) {
1.7 ! noro 852: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 853: dp_ptozp2(d,g,&t,&u); d = t; g = u;
854: hmag = multiple*HMAG(d);
855: }
856: } else {
1.7 ! noro 857: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 858: dp_ptozp(g,&t); g = t;
859: hmag = multiple*HMAG(g);
860: }
861: }
862: }
863: else if ( !full ) {
864: if ( g ) {
865: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
866: }
867: *rp = g; return;
868: } else {
869: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
870: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
871: addd(CO,d,t,&s); d = s;
872: dp_rest(g,&t); g = t;
873:
874: }
875: }
876: if ( d )
877: d->sugar = sugar;
878: *rp = d;
879: }
880:
881: void dp_nf_mod(b,g,ps,mod,full,rp)
882: NODE b;
883: DP g;
884: DP *ps;
885: int mod,full;
886: DP *rp;
887: {
888: DP u,p,d,s,t;
889: P dmy;
890: NODE l;
891: MP m,mr;
892: int sugar,psugar;
893:
894: if ( !g ) {
895: *rp = 0; return;
896: }
897: sugar = g->sugar;
898: for ( d = 0; g; ) {
899: for ( u = 0, l = b; l; l = NEXT(l) ) {
900: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
901: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
902: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
903: sugar = MAX(sugar,psugar);
904: if ( !u ) {
905: if ( d )
906: d->sugar = sugar;
907: *rp = d; return;
908: }
909: d = t;
910: break;
911: }
912: }
913: if ( u )
914: g = u;
915: else if ( !full ) {
916: if ( g ) {
917: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
918: }
919: *rp = g; return;
920: } else {
921: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
922: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
923: addmd(CO,mod,d,t,&s); d = s;
924: dp_rest(g,&t); g = t;
925: }
926: }
927: if ( d )
928: d->sugar = sugar;
929: *rp = d;
930: }
931:
932: void dp_true_nf_mod(b,g,ps,mod,full,rp,dnp)
933: NODE b;
934: DP g;
935: DP *ps;
936: int mod,full;
937: DP *rp;
938: P *dnp;
939: {
940: DP u,p,d,s,t;
941: NODE l;
942: MP m,mr;
943: int i,n;
944: int *wb;
945: int sugar,psugar;
946: P dn,tdn,tdn1;
947:
948: dn = (P)ONEM;
949: if ( !g ) {
950: *rp = 0; *dnp = dn; 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: sugar = g->sugar;
957: for ( d = 0; g; ) {
958: for ( u = 0, i = 0; i < n; i++ ) {
959: if ( dp_redble(g,p = ps[wb[i]]) ) {
960: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
961: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
962: sugar = MAX(sugar,psugar);
963: if ( !u ) {
964: if ( d )
965: d->sugar = sugar;
966: *rp = d; *dnp = dn; return;
967: } else {
968: d = t;
969: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
970: }
971: break;
972: }
973: }
974: if ( u )
975: g = u;
976: else if ( !full ) {
977: if ( g ) {
978: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
979: }
980: *rp = g; *dnp = dn; return;
981: } else {
982: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
983: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
984: addmd(CO,mod,d,t,&s); d = s;
985: dp_rest(g,&t); g = t;
986: }
987: }
988: if ( d )
989: d->sugar = sugar;
990: *rp = d; *dnp = dn;
991: }
992:
1.7 ! noro 993: void _dp_nf_mod_destructive(b,g,ps,mod,full,rp)
! 994: NODE b;
! 995: DP g;
! 996: DP *ps;
! 997: int mod,full;
! 998: DP *rp;
1.5 noro 999: {
1.7 ! noro 1000: DP u,p,d,s,t;
! 1001: NODE l;
! 1002: MP m,mr,mrd;
! 1003: int sugar,psugar,n,h_reducible,i;
1.5 noro 1004:
1.7 ! noro 1005: if ( !g ) {
! 1006: *rp = 0; return;
1.5 noro 1007: }
1.7 ! noro 1008: sugar = g->sugar;
! 1009: n = g->nv;
! 1010: for ( d = 0; g; ) {
! 1011: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
! 1012: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
! 1013: h_reducible = 1;
! 1014: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1015: _dp_red_mod_destructive(g,p,mod,&u); g = u;
! 1016: sugar = MAX(sugar,psugar);
! 1017: if ( !g ) {
! 1018: if ( d )
! 1019: d->sugar = sugar;
! 1020: _dptodp(d,rp); _free_dp(d); return;
! 1021: }
! 1022: break;
! 1023: }
! 1024: }
! 1025: if ( !h_reducible ) {
! 1026: /* head term is not reducible */
! 1027: if ( !full ) {
! 1028: if ( g )
! 1029: g->sugar = sugar;
! 1030: _dptodp(g,rp); _free_dp(g); return;
! 1031: } else {
! 1032: m = BDY(g);
! 1033: if ( NEXT(m) ) {
! 1034: BDY(g) = NEXT(m); NEXT(m) = 0;
! 1035: } else {
! 1036: _FREEDP(g); g = 0;
! 1037: }
! 1038: if ( d ) {
! 1039: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
! 1040: NEXT(mrd) = m;
! 1041: } else {
! 1042: _MKDP(n,m,d);
! 1043: }
! 1044: }
! 1045: }
1.5 noro 1046: }
1.7 ! noro 1047: if ( d )
! 1048: d->sugar = sugar;
! 1049: _dptodp(d,rp); _free_dp(d);
1.5 noro 1050: }
1051:
1.7 ! noro 1052: void dp_lnf_mod(p1,p2,g,mod,r1p,r2p)
! 1053: DP p1,p2;
! 1054: NODE g;
! 1055: int mod;
! 1056: DP *r1p,*r2p;
1.5 noro 1057: {
1.7 ! noro 1058: DP r1,r2,b1,b2,t,s;
! 1059: P c;
! 1060: MQ c1,c2;
! 1061: NODE l,b;
! 1062: int n;
! 1063:
! 1064: if ( !p1 ) {
! 1065: *r1p = p1; *r2p = p2; return;
! 1066: }
! 1067: n = p1->nv;
! 1068: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
! 1069: if ( !r1 ) {
! 1070: *r1p = r1; *r2p = r2; return;
! 1071: }
! 1072: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
! 1073: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
! 1074: b2 = (DP)BDY(NEXT(b));
! 1075: invmq(mod,(MQ)BDY(b1)->c,&c1);
! 1076: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
! 1077: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
! 1078: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
! 1079: }
! 1080: }
! 1081: *r1p = r1; *r2p = r2;
1.5 noro 1082: }
1083:
1.7 ! noro 1084: void dp_nf_tab_mod(p,tab,mod,rp)
! 1085: DP p;
! 1086: LIST *tab;
! 1087: int mod;
! 1088: DP *rp;
1.5 noro 1089: {
1.7 ! noro 1090: DP s,t,u;
! 1091: MP m;
! 1092: DL h;
! 1093: int i,n;
! 1094:
! 1095: if ( !p ) {
! 1096: *rp = p; return;
! 1097: }
! 1098: n = p->nv;
! 1099: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
! 1100: h = m->dl;
! 1101: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
! 1102: i++;
! 1103: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
! 1104: addmd(CO,mod,s,t,&u); s = u;
! 1105: }
! 1106: *rp = s;
1.5 noro 1107: }
1108:
1.7 ! noro 1109: /*
! 1110: * setting flags
! 1111: *
! 1112: */
! 1113:
! 1114: int create_order_spec(obj,spec)
! 1115: Obj obj;
! 1116: struct order_spec *spec;
1.5 noro 1117: {
1.7 ! noro 1118: int i,j,n,s,row,col;
! 1119: struct order_pair *l;
! 1120: NODE node,t,tn;
! 1121: MAT m;
! 1122: pointer **b;
! 1123: int **w;
1.5 noro 1124:
1.7 ! noro 1125: if ( !obj || NUM(obj) ) {
! 1126: spec->id = 0; spec->obj = obj;
! 1127: spec->ord.simple = QTOS((Q)obj);
! 1128: return 1;
! 1129: } else if ( OID(obj) == O_LIST ) {
! 1130: node = BDY((LIST)obj);
! 1131: for ( n = 0, t = node; t; t = NEXT(t), n++ );
! 1132: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
! 1133: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
! 1134: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
! 1135: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
! 1136: s += l[i].length;
! 1137: }
! 1138: spec->id = 1; spec->obj = obj;
! 1139: spec->ord.block.order_pair = l;
! 1140: spec->ord.block.length = n; spec->nv = s;
! 1141: return 1;
! 1142: } else if ( OID(obj) == O_MAT ) {
! 1143: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
! 1144: w = almat(row,col);
! 1145: for ( i = 0; i < row; i++ )
! 1146: for ( j = 0; j < col; j++ )
! 1147: w[i][j] = QTOS((Q)b[i][j]);
! 1148: spec->id = 2; spec->obj = obj;
! 1149: spec->nv = col; spec->ord.matrix.row = row;
! 1150: spec->ord.matrix.matrix = w;
! 1151: return 1;
! 1152: } else
1.5 noro 1153: return 0;
1154: }
1155:
1.7 ! noro 1156: /*
! 1157: * converters
! 1158: *
! 1159: */
! 1160:
! 1161: void dp_homo(p,rp)
! 1162: DP p;
! 1163: DP *rp;
1.5 noro 1164: {
1.7 ! noro 1165: MP m,mr,mr0;
! 1166: int i,n,nv,td;
! 1167: DL dl,dlh;
1.5 noro 1168:
1.7 ! noro 1169: if ( !p )
! 1170: *rp = 0;
! 1171: else {
! 1172: n = p->nv; nv = n + 1;
! 1173: m = BDY(p); td = sugard(m);
! 1174: for ( mr0 = 0; m; m = NEXT(m) ) {
! 1175: NEXTMP(mr0,mr); mr->c = m->c;
! 1176: dl = m->dl;
! 1177: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
! 1178: dlh->td = td;
! 1179: for ( i = 0; i < n; i++ )
! 1180: dlh->d[i] = dl->d[i];
! 1181: dlh->d[n] = td - dl->td;
! 1182: }
! 1183: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 1184: }
1185: }
1186:
1.7 ! noro 1187: void dp_dehomo(p,rp)
! 1188: DP p;
1.5 noro 1189: DP *rp;
1190: {
1.7 ! noro 1191: MP m,mr,mr0;
! 1192: int i,n,nv;
! 1193: DL dl,dlh;
1.5 noro 1194:
1.7 ! noro 1195: if ( !p )
! 1196: *rp = 0;
! 1197: else {
! 1198: n = p->nv; nv = n - 1;
! 1199: m = BDY(p);
! 1200: for ( mr0 = 0; m; m = NEXT(m) ) {
! 1201: NEXTMP(mr0,mr); mr->c = m->c;
! 1202: dlh = m->dl;
! 1203: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
! 1204: dl->td = dlh->td - dlh->d[nv];
! 1205: for ( i = 0; i < nv; i++ )
! 1206: dl->d[i] = dlh->d[i];
! 1207: }
! 1208: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 1209: }
1.5 noro 1210: }
1211:
1.7 ! noro 1212: void dp_mod(p,mod,subst,rp)
! 1213: DP p;
! 1214: int mod;
! 1215: NODE subst;
1.5 noro 1216: DP *rp;
1217: {
1.7 ! noro 1218: MP m,mr,mr0;
! 1219: P t,s,s1;
! 1220: V v;
! 1221: NODE tn;
1.5 noro 1222:
1.7 ! noro 1223: if ( !p )
! 1224: *rp = 0;
! 1225: else {
! 1226: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 1227: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
! 1228: v = VR((P)BDY(tn)); tn = NEXT(tn);
! 1229: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
! 1230: }
! 1231: ptomp(mod,s,&t);
! 1232: if ( t ) {
! 1233: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
! 1234: }
! 1235: }
! 1236: if ( mr0 ) {
! 1237: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 1238: } else
! 1239: *rp = 0;
! 1240: }
1.5 noro 1241: }
1242:
1.7 ! noro 1243: void dp_rat(p,rp)
! 1244: DP p;
! 1245: DP *rp;
1.5 noro 1246: {
1.7 ! noro 1247: MP m,mr,mr0;
1.5 noro 1248:
1.7 ! noro 1249: if ( !p )
! 1250: *rp = 0;
! 1251: else {
! 1252: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 1253: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
1.5 noro 1254: }
1.7 ! noro 1255: if ( mr0 ) {
! 1256: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 1257: } else
! 1258: *rp = 0;
1.5 noro 1259: }
1260: }
1261:
1262:
1.7 ! noro 1263: void homogenize_order(old,n,new)
! 1264: struct order_spec *old,*new;
! 1265: int n;
1.5 noro 1266: {
1.7 ! noro 1267: struct order_pair *l;
! 1268: int length,nv,row,i,j;
! 1269: int **newm,**oldm;
1.5 noro 1270:
1.7 ! noro 1271: switch ( old->id ) {
! 1272: case 0:
! 1273: switch ( old->ord.simple ) {
! 1274: case 0:
! 1275: new->id = 0; new->ord.simple = 0; break;
! 1276: case 1:
! 1277: l = (struct order_pair *)
! 1278: MALLOC_ATOMIC(2*sizeof(struct order_pair));
! 1279: l[0].length = n; l[0].order = 1;
! 1280: l[1].length = 1; l[1].order = 2;
! 1281: new->id = 1;
! 1282: new->ord.block.order_pair = l;
! 1283: new->ord.block.length = 2; new->nv = n+1;
! 1284: break;
! 1285: case 2:
! 1286: new->id = 0; new->ord.simple = 1; break;
! 1287: case 3: case 4: case 5:
! 1288: new->id = 0; new->ord.simple = old->ord.simple+3;
! 1289: dp_nelim = n-1; break;
! 1290: case 6: case 7: case 8: case 9:
! 1291: new->id = 0; new->ord.simple = old->ord.simple; break;
! 1292: default:
! 1293: error("homogenize_order : invalid input");
! 1294: }
! 1295: break;
! 1296: case 1:
! 1297: length = old->ord.block.length;
! 1298: l = (struct order_pair *)
! 1299: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
! 1300: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
! 1301: l[length].order = 2; l[length].length = 1;
! 1302: new->id = 1; new->nv = n+1;
! 1303: new->ord.block.order_pair = l;
! 1304: new->ord.block.length = length+1;
! 1305: break;
! 1306: case 2:
! 1307: nv = old->nv; row = old->ord.matrix.row;
! 1308: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
! 1309: for ( i = 0; i <= nv; i++ )
! 1310: newm[0][i] = 1;
! 1311: for ( i = 0; i < row; i++ ) {
! 1312: for ( j = 0; j < nv; j++ )
! 1313: newm[i+1][j] = oldm[i][j];
! 1314: newm[i+1][j] = 0;
! 1315: }
! 1316: new->id = 2; new->nv = nv+1;
! 1317: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
! 1318: break;
! 1319: default:
! 1320: error("homogenize_order : invalid input");
1.5 noro 1321: }
1.7 ! noro 1322: }
! 1323:
! 1324: void qltozl(w,n,dvr)
! 1325: Q *w,*dvr;
! 1326: int n;
! 1327: {
! 1328: N nm,dn;
! 1329: N g,l1,l2,l3;
! 1330: Q c,d;
! 1331: int i;
! 1332: struct oVECT v;
1.5 noro 1333:
1334: for ( i = 0; i < n; i++ )
1.7 ! noro 1335: if ( w[i] && !INT(w[i]) )
! 1336: break;
! 1337: if ( i == n ) {
! 1338: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
! 1339: igcdv(&v,dvr); return;
! 1340: }
! 1341: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
! 1342: for ( i = 1; i < n; i++ ) {
! 1343: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
! 1344: gcdn(nm,NM(c),&g); nm = g;
! 1345: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 1346: }
1.7 ! noro 1347: if ( UNIN(dn) )
! 1348: NTOQ(nm,1,d);
! 1349: else
! 1350: NDTOQ(nm,dn,1,d);
! 1351: *dvr = d;
! 1352: }
1.5 noro 1353:
1.7 ! noro 1354: int comp_nm(a,b)
! 1355: Q *a,*b;
! 1356: {
! 1357: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
! 1358: }
! 1359:
! 1360: void sortbynm(w,n)
! 1361: Q *w;
! 1362: int n;
! 1363: {
! 1364: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
! 1365: }
1.5 noro 1366:
1367:
1.7 ! noro 1368: /*
! 1369: * simple operations
! 1370: *
! 1371: */
1.5 noro 1372:
1.7 ! noro 1373: int dp_redble(p1,p2)
! 1374: DP p1,p2;
! 1375: {
! 1376: int i,n;
! 1377: DL d1,d2;
1.5 noro 1378:
1.7 ! noro 1379: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 1380: if ( d1->td < d2->td )
! 1381: return 0;
! 1382: else {
! 1383: for ( i = 0, n = p1->nv; i < n; i++ )
! 1384: if ( d1->d[i] < d2->d[i] )
! 1385: return 0;
! 1386: return 1;
1.5 noro 1387: }
1388: }
1389:
1.7 ! noro 1390: void dp_subd(p1,p2,rp)
1.5 noro 1391: DP p1,p2;
1392: DP *rp;
1393: {
1.7 ! noro 1394: int i,n;
1.5 noro 1395: DL d1,d2,d;
1396: MP m;
1.7 ! noro 1397: DP s;
1.5 noro 1398:
1399: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 ! noro 1400: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 1401: for ( i = 0; i < n; i++ )
1.7 ! noro 1402: d->d[i] = d1->d[i]-d2->d[i];
! 1403: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
! 1404: MKDP(n,m,s); s->sugar = d->td;
! 1405: *rp = s;
! 1406: }
! 1407:
! 1408: void dltod(d,n,rp)
! 1409: DL d;
! 1410: int n;
! 1411: DP *rp;
! 1412: {
! 1413: MP m;
! 1414: DP s;
! 1415:
! 1416: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
! 1417: MKDP(n,m,s); s->sugar = d->td;
! 1418: *rp = s;
1.5 noro 1419: }
1420:
1421: void dp_hm(p,rp)
1422: DP p;
1423: DP *rp;
1424: {
1425: MP m,mr;
1426:
1427: if ( !p )
1428: *rp = 0;
1429: else {
1430: m = BDY(p);
1431: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
1432: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1433: }
1434: }
1435:
1436: void dp_rest(p,rp)
1437: DP p,*rp;
1438: {
1439: MP m;
1440:
1441: m = BDY(p);
1442: if ( !NEXT(m) )
1443: *rp = 0;
1444: else {
1445: MKDP(p->nv,NEXT(m),*rp);
1446: if ( *rp )
1447: (*rp)->sugar = p->sugar;
1448: }
1449: }
1450:
1451: DL lcm_of_DL(nv,dl1,dl2,dl)
1452: int nv;
1453: DL dl1,dl2;
1454: register DL dl;
1455: {
1456: register int n, *d1, *d2, *d, td;
1457:
1458: if ( !dl ) NEWDL(dl,nv);
1459: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1460: for ( td = 0, n = nv; --n >= 0; d1++, d2++, d++ )
1461: td += (*d = *d1 > *d2 ? *d1 : *d2 );
1462: dl->td = td;
1463: return dl;
1464: }
1465:
1466: int dl_equal(nv,dl1,dl2)
1467: int nv;
1468: DL dl1, dl2;
1469: {
1470: register int *d1, *d2, n;
1471:
1472: if ( dl1->td != dl2->td ) return 0;
1473: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
1474: if ( *d1 != *d2 ) return 0;
1475: return 1;
1476: }
1477:
1478: int dp_nt(p)
1479: DP p;
1480: {
1481: int i;
1482: MP m;
1483:
1484: if ( !p )
1485: return 0;
1486: else {
1487: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
1488: return i;
1489: }
1490: }
1491: