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