![[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.1 noro 1: /* $OpenXM: OpenXM/src/asir99/builtin/dp-supp.c,v 1.2 1999/11/18 05:42:01 noro Exp $ */
2: #include "ca.h"
3: #include "base.h"
4: #include "parse.h"
5: #include "ox.h"
6:
7: extern int (*cmpdl)();
8: double pz_t_e,pz_t_d,pz_t_d1,pz_t_c;
9:
10: void Pdp_mbase(),Pdp_lnf_mod(),Pdp_nf_tab_mod(),Pdp_mdtod();
11: void Pdp_vtoe(), Pdp_etov(), Pdp_dtov(), Pdp_idiv(), Pdp_sep();
12: void Pdp_cont(), Pdp_igcdv_hist(), Pdp_idivv_hist();
13:
14: struct ftab dp_supp_tab[] = {
15: {"dp_mbase",Pdp_mbase,1},
16: {"dp_lnf_mod",Pdp_lnf_mod,3},
17: {"dp_nf_tab_mod",Pdp_nf_tab_mod,3},
18: {"dp_etov",Pdp_etov,1},
19: {"dp_vtoe",Pdp_vtoe,1},
20: {"dp_dtov",Pdp_dtov,1},
21: {"dp_idiv",Pdp_idiv,2},
22: {"dp_cont",Pdp_cont,1},
23: {"dp_sep",Pdp_sep,2},
24: {"dp_mdtod",Pdp_mdtod,1},
25: {"dp_igcdv_hist",Pdp_igcdv_hist,1},
26: {"dp_idivv_hist",Pdp_idivv_hist,1},
27: {0,0,0}
28: };
29:
30: void Pdp_mdtod(arg,rp)
31: NODE arg;
32: DP *rp;
33: {
34: MP m,mr,mr0;
35: DP p;
36: P t;
37:
38: p = (DP)ARG0(arg);
39: if ( !p )
40: *rp = 0;
41: else {
42: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
43: mptop(m->c,&t); NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
44: }
45: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
46: }
47: }
48:
49: void Pdp_sep(arg,rp)
50: NODE arg;
51: VECT *rp;
52: {
53: DP p,r;
54: MP m,t;
55: MP *w0,*w;
56: int i,n,d,nv,sugar;
57: VECT v;
58: pointer *pv;
59:
60: p = (DP)ARG0(arg); m = BDY(p);
61: d = QTOS((Q)ARG1(arg));
62: for ( t = m, n = 0; t; t = NEXT(t), n++ );
63: if ( d > n )
64: d = n;
65: MKVECT(v,d); *rp = v;
66: pv = BDY(v); nv = p->nv; sugar = p->sugar;
67: w0 = (MP *)MALLOC(d*sizeof(MP)); bzero(w0,d*sizeof(MP));
68: w = (MP *)MALLOC(d*sizeof(MP)); bzero(w,d*sizeof(MP));
69: for ( t = BDY(p), i = 0; t; t = NEXT(t), i++, i %= d ) {
70: NEXTMP(w0[i],w[i]); w[i]->c = t->c; w[i]->dl = t->dl;
71: }
72: for ( i = 0; i < d; i++ ) {
73: NEXT(w[i]) = 0; MKDP(nv,w0[i],r); r->sugar = sugar;
74: pv[i] = (pointer)r;
75: }
76: }
77:
78: void Pdp_idiv(arg,rp)
79: NODE arg;
80: DP *rp;
81: {
82: dp_idiv((DP)ARG0(arg),(Q)ARG1(arg),rp);
83: }
84:
85: void Pdp_cont(arg,rp)
86: NODE arg;
87: Q *rp;
88: {
89: dp_cont((DP)ARG0(arg),rp);
90: }
91:
92: void Pdp_dtov(arg,rp)
93: NODE arg;
94: VECT *rp;
95: {
96: dp_dtov((DP)ARG0(arg),rp);
97: }
98:
99: void Pdp_mbase(arg,rp)
100: NODE arg;
101: LIST *rp;
102: {
103: NODE mb;
104:
105: asir_assert(ARG0(arg),O_LIST,"dp_mbase");
106: dp_mbase(BDY((LIST)ARG0(arg)),&mb);
107: MKLIST(*rp,mb);
108: }
109:
110: void Pdp_etov(arg,rp)
111: NODE arg;
112: VECT *rp;
113: {
114: DP dp;
115: int n,i;
116: int *d;
117: VECT v;
118: Q t;
119:
120: dp = (DP)ARG0(arg);
121: asir_assert(dp,O_DP,"dp_etov");
122: n = dp->nv; d = BDY(dp)->dl->d;
123: MKVECT(v,n);
124: for ( i = 0; i < n; i++ ) {
125: STOQ(d[i],t); v->body[i] = (pointer)t;
126: }
127: *rp = v;
128: }
129:
130: void Pdp_vtoe(arg,rp)
131: NODE arg;
132: DP *rp;
133: {
134: DP dp;
135: DL dl;
136: MP m;
137: int n,i,td;
138: int *d;
139: VECT v;
140:
141: v = (VECT)ARG0(arg);
142: asir_assert(v,O_VECT,"dp_vtoe");
143: n = v->len;
144: NEWDL(dl,n); d = dl->d;
145: for ( i = 0, td = 0; i < n; i++ ) {
146: d[i] = QTOS((Q)(v->body[i])); td += d[i];
147: }
148: dl->td = td;
149: NEWMP(m); m->dl = dl; m->c = (P)ONE; NEXT(m) = 0;
150: MKDP(n,m,dp); dp->sugar = td;
151: *rp = dp;
152: }
153:
154: void Pdp_lnf_mod(arg,rp)
155: NODE arg;
156: LIST *rp;
157: {
158: DP r1,r2;
159: NODE b,g,n;
160: int mod;
161:
162: asir_assert(ARG0(arg),O_LIST,"dp_lnf_mod");
163: asir_assert(ARG1(arg),O_LIST,"dp_lnf_mod");
164: asir_assert(ARG2(arg),O_N,"dp_lnf_mod");
165: b = BDY((LIST)ARG0(arg)); g = BDY((LIST)ARG1(arg));
166: mod = QTOS((Q)ARG2(arg));
167: dp_lnf_mod((DP)BDY(b),(DP)BDY(NEXT(b)),g,mod,&r1,&r2);
168: NEWNODE(n); BDY(n) = (pointer)r1;
169: NEWNODE(NEXT(n)); BDY(NEXT(n)) = (pointer)r2;
170: NEXT(NEXT(n)) = 0; MKLIST(*rp,n);
171: }
172:
173: void Pdp_nf_tab_mod(arg,rp)
174: NODE arg;
175: DP *rp;
176: {
177: asir_assert(ARG0(arg),O_DP,"dp_nf_tab_mod");
178: asir_assert(ARG1(arg),O_VECT,"dp_nf_tab_mod");
179: asir_assert(ARG2(arg),O_N,"dp_nf_tab_mod");
180: dp_nf_tab_mod((DP)ARG0(arg),(LIST *)BDY((VECT)ARG1(arg)),
181: QTOS((Q)ARG2(arg)),rp);
182: }
183:
184: void Pdp_igcdv_hist(arg,rp)
185: NODE arg;
186: Q *rp;
187: {
188: dp_igcdv_hist((DP)ARG0(arg),rp);
189: }
190:
191: void Pdp_idivv_hist(arg,rp)
192: NODE arg;
193: DP *rp;
194: {
195: dp_idivv_hist((Q)ARG0(arg),rp);
196: }
197:
198: void dp_idiv(p,c,rp)
199: DP p;
200: Q c;
201: DP *rp;
202: {
203: Q t;
204: N nm,q;
205: int sgn,s;
206: MP mr0,m,mr;
207:
208: if ( !p )
209: *rp = 0;
210: else if ( MUNIQ((Q)c) )
211: *rp = p;
212: else if ( MUNIQ((Q)c) )
213: chsgnd(p,rp);
214: else {
215: nm = NM(c); sgn = SGN(c);
216: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
217: NEXTMP(mr0,mr);
218:
219: divsn(NM((Q)(m->c)),nm,&q);
220: s = sgn*SGN((Q)(m->c));
221: NTOQ(q,s,t);
222: mr->c = (P)t;
223: mr->dl = m->dl;
224: }
225: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
226: if ( *rp )
227: (*rp)->sugar = p->sugar;
228: }
229: }
230:
231: void dp_cont(p,rp)
232: DP p;
233: Q *rp;
234: {
235: VECT v;
236:
237: dp_dtov(p,&v); igcdv(v,rp);
238: }
239:
240: void dp_dtov(dp,rp)
241: DP dp;
242: VECT *rp;
243: {
244: MP m,t;
245: int i,n;
246: VECT v;
247: pointer *p;
248:
249: m = BDY(dp);
250: for ( t = m, n = 0; t; t = NEXT(t), n++ );
251: MKVECT(v,n);
252: for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
253: p[i] = (pointer)(t->c);
254: *rp = v;
255: }
256:
257: void dp_mbase(hlist,mbase)
258: NODE hlist;
259: NODE *mbase;
260: {
261: DL *dl;
262: DL d;
263: int i,j,n,nvar,td;
264:
265: n = length(hlist); nvar = ((DP)BDY(hlist))->nv;
266: dl = (DL *)MALLOC(n*sizeof(DL));
267: for ( i = 0; i < n; i++, hlist = NEXT(hlist) )
268: dl[i] = BDY((DP)BDY(hlist))->dl;
269: NEWDL(d,nvar); *mbase = 0;
270: while ( 1 ) {
271: insert_to_node(d,mbase,nvar);
272: for ( i = nvar-1; i >= 0; ) {
273: d->d[i]++; d->td++;
274: for ( j = 0; j < n; j++ ) {
275: if ( _dl_redble(dl[j],d,nvar) )
276: break;
277: }
278: if ( j < n ) {
279: for ( j = nvar-1; j >= i; j-- )
280: d->d[j] = 0;
281: for ( j = 0, td = 0; j < i; j++ )
282: td += d->d[j];
283: d->td = td;
284: i--;
285: } else
286: break;
287: }
288: if ( i < 0 )
289: break;
290: }
291: }
292:
293: int _dl_redble(d1,d2,nvar)
294: DL d1,d2;
295: int nvar;
296: {
297: int i;
298:
299: if ( d1->td > d2->td )
300: return 0;
301: for ( i = 0; i < nvar; i++ )
302: if ( d1->d[i] > d2->d[i] )
303: break;
304: if ( i < nvar )
305: return 0;
306: else
307: return 1;
308: }
309:
310: void insert_to_node(d,n,nvar)
311: DL d;
312: NODE *n;
313: int nvar;
314: {
315: DL d1;
316: MP m;
317: DP dp;
318: NODE n0,n1,n2;
319:
320: NEWDL(d1,nvar); d1->td = d->td;
321: bcopy((char *)d->d,(char *)d1->d,nvar*sizeof(int));
322: NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0;
323: MKDP(nvar,m,dp); dp->sugar = d->td;
324: if ( !(*n) ) {
325: MKNODE(n1,dp,0); *n = n1;
326: } else {
327: for ( n1 = *n, n0 = 0; n1; n0 = n1, n1 = NEXT(n1) )
328: if ( (*cmpdl)(nvar,d,BDY((DP)BDY(n1))->dl) > 0 ) {
329: MKNODE(n2,dp,n1);
330: if ( !n0 )
331: *n = n2;
332: else
333: NEXT(n0) = n2;
334: break;
335: }
336: if ( !n1 ) {
337: MKNODE(n2,dp,0); NEXT(n0) = n2;
338: }
339: }
340: }
341:
342: void dp_lnf_mod(p1,p2,g,mod,r1p,r2p)
343: DP p1,p2;
344: NODE g;
345: int mod;
346: DP *r1p,*r2p;
347: {
348: DP r1,r2,b1,b2,t,s;
349: P c;
350: MQ c1,c2;
351: NODE l,b;
352: int n;
353:
354: if ( !p1 ) {
355: *r1p = p1; *r2p = p2; return;
356: }
357: n = p1->nv;
358: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
359: if ( !r1 ) {
360: *r1p = r1; *r2p = r2; return;
361: }
362: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
363: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
364: b2 = (DP)BDY(NEXT(b));
365: invmq(mod,(MQ)BDY(b1)->c,&c1);
366: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
367: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
368: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
369: }
370: }
371: *r1p = r1; *r2p = r2;
372: }
373:
374: void dp_nf_tab_mod(p,tab,mod,rp)
375: DP p;
376: LIST *tab;
377: int mod;
378: DP *rp;
379: {
380: DP s,t,u;
381: MP m;
382: DL h;
383: int i,n;
384:
385: if ( !p ) {
386: *rp = p; return;
387: }
388: n = p->nv;
389: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
390: h = m->dl;
391: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
392: i++;
393: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
394: addmd(CO,mod,s,t,&u); s = u;
395: }
396: *rp = s;
397: }
398:
399: void dp_ptozp(p,rp)
400: DP p,*rp;
401: {
402: MP m,mr,mr0;
403: int i,n;
404: Q *w;
405: Q dvr;
406: P t;
407:
408: if ( !p )
409: *rp = 0;
410: else {
411: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
412: w = (Q *)ALLOCA(n*sizeof(Q));
413: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
414: if ( NUM(m->c) )
415: w[i] = (Q)m->c;
416: else
417: ptozp(m->c,1,&w[i],&t);
418: sortbynm(w,n);
419: qltozl(w,n,&dvr);
420: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
421: NEXTMP(mr0,mr); divsp(CO,m->c,(P)dvr,&mr->c); mr->dl = m->dl;
422: }
423: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
424: }
425: }
426:
427: void dp_ptozp2(p0,p1,hp,rp)
428: DP p0,p1;
429: DP *hp,*rp;
430: {
431: DP t,s,h,r;
432: MP m,mr,mr0,m0;
433:
434: addd(CO,p0,p1,&t); dp_ptozp(t,&s);
435: if ( !p0 ) {
436: h = 0; r = s;
437: } else if ( !p1 ) {
438: h = s; r = 0;
439: } else {
440: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
441: m = NEXT(m), m0 = NEXT(m0) ) {
442: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
443: }
444: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
445: }
446: if ( h )
447: h->sugar = p0->sugar;
448: if ( r )
449: r->sugar = p1->sugar;
450: *hp = h; *rp = r;
451: }
452:
453: struct igcdv_hist {
454: int len;
455: DP dp;
456: Q d,d1;
457: Q *q,*q1;
458: };
459:
460: static struct igcdv_hist current_igcdv_hist;
461:
462: void dp_igcdv_hist(p,rp)
463: DP p;
464: Q *rp;
465: {
466: VECT c;
467: Q *a,*q,*r,*q1,*r1;
468: int i,n,nz;
469: Q d,d1;
470: N qn,rn,gn;
471: struct oVECT v;
472:
473: dp_dtov(p,&c); a = (Q *)c->body;
474: n = c->len;
475: q = (Q *)MALLOC(n*sizeof(Q)); r = (Q *)MALLOC(n*sizeof(Q));
476: igcdv_estimate(c,&d);
477: for ( i = 0, nz = 0; i < n; i++ ) {
478: divn(NM(a[i]),NM(d),&qn,&rn);
479: NTOQ(qn,SGN(a[i]),q[i]); NTOQ(rn,SGN(a[i]),r[i]);
480: if ( r[i] )
481: nz = 1;
482: }
483: current_igcdv_hist.len = n;
484: current_igcdv_hist.dp = p;
485: current_igcdv_hist.d = d;
486: current_igcdv_hist.q = q;
487: if ( nz ) {
488: v.id = O_VECT; v.len = n; v.body = (pointer *)r;
489: igcdv(&v,&d1);
490: q1 = (Q *)MALLOC(n*sizeof(Q)); r1 = (Q *)MALLOC(n*sizeof(Q));
491: for ( i = 0; i < n; i++ ) {
492: if ( r[i] )
493: divn(NM(r[i]),NM(d1),&qn,&rn);
494: else {
495: qn = rn = 0;
496: }
497: NTOQ(qn,SGN(r[i]),q1[i]); NTOQ(rn,SGN(r[i]),r1[i]);
498: }
499: gcdn(NM(d),NM(d1),&gn); NTOQ(gn,1,*rp);
500: } else {
501: d1 = 0; q1 = 0; *rp = d;
502: }
503: current_igcdv_hist.d1 = d1;
504: current_igcdv_hist.q1 = q1;
505: }
506:
507: void dp_idivv_hist(g,rp)
508: Q g;
509: DP *rp;
510: {
511: int i,n;
512: Q *q,*q1,*s;
513: Q t,u,m;
514: N qn;
515:
516: n = current_igcdv_hist.len;
517: q = current_igcdv_hist.q; q1 = current_igcdv_hist.q1;
518: divsn(NM(current_igcdv_hist.d),NM(g),&qn); NTOQ(qn,1,m);
519: s = (Q *)MALLOC(n*sizeof(Q));
520: for ( i = 0; i < n; i++ )
521: mulq(m,q[i],&s[i]);
522: if ( q1 ) {
523: divsn(NM(current_igcdv_hist.d1),NM(g),&qn); NTOQ(qn,1,m);
524: for ( i = 0; i < n; i++ ) {
525: mulq(m,q1[i],&t); addq(t,s[i],&u); s[i] = u;
526: }
527: }
528: dp_vtod(s,current_igcdv_hist.dp,rp);
529: }
530:
531: void dp_vtod(c,p,rp)
532: Q *c;
533: DP p;
534: DP *rp;
535: {
536: MP mr0,m,mr;
537: int i;
538:
539: if ( !p )
540: *rp = 0;
541: else {
542: for ( mr0 = 0, m = BDY(p), i = 0; m; m = NEXT(m), i++ ) {
543: NEXTMP(mr0,mr); mr->c = (P)c[i]; mr->dl = m->dl;
544: }
545: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
546: (*rp)->sugar = p->sugar;
547: }
548: }
549:
550: #if INET
551: void dp_ptozp_d_old(dist,ndist,p,rp)
552: NODE dist;
553: int ndist;
554: DP p,*rp;
555: {
556: DP d,s,u,su;
557: MP m,mr,mr0;
558: MP *w0,*w;
559: Q *g;
560: int i,n,nv,nsep;
561: P t;
562: NODE tn,n0,n1,n2;
563: struct oVECT v;
564: int id,sindex;
565: Obj dmy;
566: Q gcd;
567: STRING cont,div;
568:
569: if ( !p )
570: *rp = 0;
571: else {
572: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
573: nsep = ndist + 1;
574: if ( n <= nsep ) {
575: dp_ptozp(p,rp); return;
576: }
577: nv = p->nv;
578: w0 = (MP *)MALLOC(nsep*sizeof(MP)); bzero(w0,nsep*sizeof(MP));
579: w = (MP *)MALLOC(nsep*sizeof(MP)); bzero(w,nsep*sizeof(MP));
580: g = (Q *)MALLOC(nsep*sizeof(Q)); bzero(g,nsep*sizeof(Q));
581: for ( m = BDY(p), i = 0; m; m = NEXT(m), i++, i %= nsep ) {
582: NEXTMP(w0[i],w[i]); w[i]->c = m->c; w[i]->dl = m->dl;
583: }
584: MKSTR(cont,"reg_dp_cont"); MKSTR(div,"reg_dp_idiv");
585:
586: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
587: NEXT(w[i]) = 0; MKDP(nv,w0[i],d); d->sugar = p->sugar;
588: MKNODE(n2,d,0); MKNODE(n1,cont,n2); MKNODE(n0,BDY(tn),n1);
589: Pox_rpc(n0,&dmy);
590: }
591: NEXT(w[i]) = 0; MKDP(nv,w0[i],d); d->sugar = p->sugar;
592: dp_cont(d,&g[i]);
593:
594: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) )
595: Pox_pop_local(tn,&g[i]);
596: v.id = O_VECT; v.len = nsep; v.body = (pointer *)g; igcdv(&v,&gcd);
597:
598: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
599: MKNODE(n2,gcd,0); MKNODE(n1,div,n2); MKNODE(n0,BDY(tn),n1);
600: Pox_rpc(n0,&dmy);
601: }
602: dp_idiv(d,gcd,&s);
603:
604: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
605: Pox_pop_local(tn,&u);
606: addd(CO,s,u,&su); s = su;
607: }
608: *rp = s;
609: }
610: }
611:
612: void dp_ptozp_d(dist,ndist,p,rp)
613: NODE dist;
614: int ndist;
615: DP p,*rp;
616: {
617: int i,j,k,l,n,nsep;
618: MP m;
619: NODE tn,n0,n1,n2,n3;
620: struct oVECT v;
621: VECT c,cs;
622: VECT qi,ri;
623: LIST *qr;
624: int s,id;
625: Obj dmy;
626: Q d0,d1,gcd,a,u,u1;
627: Q *q,*r;
628: STRING iqr_v;
629: pointer *b;
630: N qn,gn;
631: double get_rtime();
632: int blen;
633: double t0;
634: double t_e,t_d,t_d1,t_c;
635:
636: if ( !p )
637: *rp = 0;
638: else {
639: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
640: nsep = ndist + 1;
641: if ( n <= nsep ) {
642: dp_ptozp(p,rp); return;
643: }
644: t0 = get_rtime();
645: dp_dtov(p,&c);
646: igcdv_estimate(c,&d0);
647: t_e = get_rtime()-t0;
648: t0 = get_rtime();
649: dp_dtov(p,&c);
650: sepvect(c,nsep,&cs);
651: MKSTR(iqr_v,"iqr");
652: qr = (LIST *)CALLOC(nsep,sizeof(LIST));
653: q = (Q *)CALLOC(n,sizeof(Q));
654: r = (Q *)CALLOC(n,sizeof(Q));
655: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
656: MKNODE(n3,d0,0); MKNODE(n2,b[i],n3);
657: MKNODE(n1,iqr_v,n2); MKNODE(n0,BDY(tn),n1);
658: Pox_rpc(n0,&dmy);
659: }
660: iqrv(b[i],d0,&qr[i]);
661: dp_dtov(p,&c);
662: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
663: Pox_pop_local(tn,&qr[i]);
664: if ( OID(qr[i]) == O_ERR ) {
665: printexpr(CO,(Obj)qr[i]);
666: error("dp_ptozp_d : aborted");
667: }
668: }
669: t_d = get_rtime()-t0;
670: t_d1 = t_d/n;
671: t0 = get_rtime();
672: for ( i = j = 0; i < nsep; i++ ) {
673: tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
674: for ( k = 0, l = qi->len; k < l; k++, j++ ) {
675: q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
676: }
677: }
678: v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
679: if ( d1 ) {
680: gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
681: divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
682: for ( i = 0; i < n; i++ ) {
683: mulq(a,q[i],&u);
684: if ( r[i] ) {
685: divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
686: addq(u,u1,&q[i]);
687: } else
688: q[i] = u;
689: }
690: } else
691: gcd = d0;
692: dp_vtod(q,p,rp);
693: t_c = get_rtime()-t0;
694: blen=p_mag((P)gcd);
695: pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
696: if ( 0 )
697: fprintf(stderr,"(%d,%d)",p_mag((P)d0)-blen,blen);
698: }
699: }
700:
701: void dp_ptozp2_d(dist,ndist,p0,p1,hp,rp)
702: NODE dist;
703: int ndist;
704: DP p0,p1;
705: DP *hp,*rp;
706: {
707: DP t,s,h,r;
708: MP m,mr,mr0,m0;
709:
710: addd(CO,p0,p1,&t); dp_ptozp_d(dist,ndist,t,&s);
711: if ( !p0 ) {
712: h = 0; r = s;
713: } else if ( !p1 ) {
714: h = s; r = 0;
715: } else {
716: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
717: m = NEXT(m), m0 = NEXT(m0) ) {
718: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
719: }
720: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
721: }
722: if ( h )
723: h->sugar = p0->sugar;
724: if ( r )
725: r->sugar = p1->sugar;
726: *hp = h; *rp = r;
727: }
728: #endif
729:
730: /*
731: * Old codes
732: */
733:
734: #if 0
735: void Pdp_sep(arg,rp)
736: NODE arg;
737: VECT *rp;
738: {
739: DP p,s;
740: MP m,t,mr,mr0;
741: int i,j,n,d,r,q,q1,nv;
742: VECT w;
743: pointer *pw;
744:
745: p = (DP)ARG0(arg); m = BDY(p);
746: d = QTOS((Q)ARG1(arg));
747: for ( t = m, n = 0; t; t = NEXT(t), n++ );
748: if ( d > n )
749: d = n;
750: q = n/d; r = n%d; q1 = q+1;
751: MKVECT(w,d); *rp = w;
752: t = m; pw = BDY(w); nv = p->nv;
753: for ( i = 0; i < r; i++ ) {
754: for ( mr0 = 0, j = 0; j < q1; j++, t = NEXT(t) ) {
755: NEXTMP(mr0,mr);
756: mr->c = t->c;
757: mr->dl = t->dl;
758: }
759: NEXT(mr) = 0; MKDP(nv,mr0,s); s->sugar = p->sugar;
760: pw[i] = (pointer)s;
761: }
762: for ( ; i < d; i++ ) {
763: for ( mr0 = 0, j = 0; j < q; j++, t = NEXT(t) ) {
764: NEXTMP(mr0,mr);
765: mr->c = t->c;
766: mr->dl = t->dl;
767: }
768: NEXT(mr) = 0; MKDP(nv,mr0,s); s->sugar = p->sugar;
769: pw[i] = (pointer)s;
770: }
771: }
772:
773: void dp_ptozp_mpi(dist,ndist,p,rp)
774: NODE dist;
775: int ndist;
776: DP p,*rp;
777: {
778: int i,j,k,l,n,nsep;
779: MP m;
780: NODE tn,n0,n1,n2,n3;
781: struct oVECT v;
782: VECT c,cs;
783: VECT qi,ri;
784: LIST *qr;
785: int id;
786: Obj dmy;
787: Q d0,d1,gcd,a,u,u1;
788: Q *q,*r;
789: STRING iqr_v;
790: pointer *b;
791: N qn,gn;
792: double get_rtime();
793: int blen;
794: double t0;
795: double t_e,t_d,t_d1,t_c;
796:
797: if ( !p )
798: *rp = 0;
799: else {
800: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
801: nsep = ndist + 1;
802: if ( n <= nsep ) {
803: dp_ptozp(p,rp); return;
804: }
805: t0 = get_rtime();
806: dp_dtov(p,&c);
807: igcdv_estimate(c,&d0);
808: t_e = get_rtime()-t0;
809: t0 = get_rtime();
810: dp_dtov(p,&c);
811: sepvect(c,nsep,&cs);
812: MKSTR(iqr_v,"iqr");
813: qr = (LIST *)CALLOC(nsep,sizeof(LIST));
814: q = (Q *)CALLOC(n,sizeof(Q));
815: r = (Q *)CALLOC(n,sizeof(Q));
816: for ( i = 0, tn = dist, b = BDY(cs); i < ndist;
817: i++, tn = NEXT(tn) ) {
818: div_vect_mpi(QTOS((Q)BDY(tn)),b[i],d0);
819: }
820: iqrv(b[i],d0,&qr[i]);
821: dp_dtov(p,&c);
822: for ( i = j = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
823: get_div_vect_mpi(QTOS((Q)BDY(tn)),&qr[i]);
824: }
825: t_d = get_rtime()-t0;
826: t_d1 = t_d/n;
827: t0 = get_rtime();
828: for ( i = j = 0; i < nsep; i++ ) {
829: tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
830: for ( k = 0, l = qi->len; k < l; k++, j++ ) {
831: q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
832: }
833: }
834: v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
835: if ( d1 ) {
836: gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
837: divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
838: for ( i = 0; i < n; i++ ) {
839: mulq(a,q[i],&u);
840: if ( r[i] ) {
841: divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
842: addq(u,u1,&q[i]);
843: } else
844: q[i] = u;
845: }
846: } else
847: gcd = d0;
848: dp_vtod(q,p,rp);
849: t_c = get_rtime()-t0;
850: blen=p_mag(gcd);
851: pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
852: if ( 0 )
853: fprintf(stderr,"(%d,%d)",p_mag(d0)-blen,blen);
854: }
855: }
856:
857: void dp_ptozp_d(dist,ndist,p,rp)
858: NODE dist;
859: int ndist;
860: DP p,*rp;
861: {
862: DP d,s,u,su;
863: MP m,mr,mr0;
864: MP *w0,*w;
865: Q *g;
866: int i,n,nv,nsep;
867: P t;
868: NODE tn,n0,n1,n2;
869: struct oVECT v;
870: int sindex,id;
871: Obj dmy;
872: Q gcd;
873: STRING gcd_h,div_h;
874:
875: if ( !p )
876: *rp = 0;
877: else {
878: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
879: nsep = ndist + 1;
880: if ( n <= nsep ) {
881: dp_ptozp(p,rp); return;
882: }
883: nv = p->nv;
884: w0 = (MP *)MALLOC(nsep*sizeof(MP)); bzero(w0,nsep*sizeof(MP));
885: w = (MP *)MALLOC(nsep*sizeof(MP)); bzero(w,nsep*sizeof(MP));
886: g = (Q *)MALLOC(nsep*sizeof(Q)); bzero(g,nsep*sizeof(Q));
887: for ( m = BDY(p), i = 0; m; m = NEXT(m), i++, i %= nsep ) {
888: NEXTMP(w0[i],w[i]); w[i]->c = m->c; w[i]->dl = m->dl;
889: }
890: MKSTR(gcd_h,"dp_igcdv_hist"); MKSTR(div_h,"dp_idivv_hist");
891:
892: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
893: NEXT(w[i]) = 0; MKDP(nv,w0[i],d); d->sugar = p->sugar;
894: MKNODE(n2,d,0); MKNODE(n1,gcd_h,n2); MKNODE(n0,BDY(tn),n1);
895: Pox_rpc(n0,&dmy);
896: }
897: NEXT(w[i]) = 0; MKDP(nv,w0[i],d); d->sugar = p->sugar;
898: dp_igcdv_hist(d,&g[i]);
899:
900: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) )
901: Pox_pop_local(tn,&g[i]);
902: v.id = O_VECT; v.len = nsep; v.body = (pointer *)g; igcdv(&v,&gcd);
903:
904: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
905: MKNODE(n2,gcd,0); MKNODE(n1,div_h,n2); MKNODE(n0,BDY(tn),n1);
906: Pox_rpc(n0,&dmy);
907: }
908: dp_idivv_hist(gcd,&s);
909:
910: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
911: Pox_pop_local(tn,&u);
912: addd(CO,s,u,&su); s = su;
913: }
914: *rp = s;
915: }
916: }
917: #endif