![[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.6 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.5 2000/12/05 06:59:15 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.6 ! noro 1540: }
! 1541:
! 1542: void _dp_sp_mod(p1,p2,mod,rp)
! 1543: DP p1,p2;
! 1544: int mod;
! 1545: DP *rp;
! 1546: {
! 1547: int i,n,td;
! 1548: int *w;
! 1549: DL d1,d2,d;
! 1550: MP m;
! 1551: DP t,s,u;
! 1552:
! 1553: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 1554: w = (int *)ALLOCA(n*sizeof(int));
! 1555: for ( i = 0, td = 0; i < n; i++ ) {
! 1556: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
! 1557: }
! 1558: NEWDL(d,n); d->td = td - d1->td;
! 1559: for ( i = 0; i < n; i++ )
! 1560: d->d[i] = w[i] - d1->d[i];
! 1561: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
! 1562: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
! 1563: NEWDL(d,n); d->td = td - d2->td;
! 1564: for ( i = 0; i < n; i++ )
! 1565: d->d[i] = w[i] - d2->d[i];
! 1566: NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
! 1567: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
! 1568: addmd_destructive(mod,t,u,rp);
1.1 noro 1569: }