![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gr.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 gr.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.4 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.3 2000/04/20 02:20:15 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "parse.h"
4: #include "base.h"
5: #include "ox.h"
6:
7: #define ITOS(p) (((unsigned int)(p))&0x7fffffff)
8: #define STOI(i) ((P)((unsigned int)(i)|0x80000000))
9:
10: #define NEXTVL(r,c) \
11: if(!(r)){NEWVL(r);(c)=(r);}else{NEWVL(NEXT(c));(c)=NEXT(c);}
12:
13: #define HMAG(p) (p_mag(BDY(p)->c))
14:
15: struct dp_pairs {
16: int dp1, dp2;
17: DL lcm;
18: int sugar;
19: struct dp_pairs *next;
20: };
21:
22: typedef struct dp_pairs *DP_pairs;
23:
24: #define NEWDPP(a) ((a)=(DP_pairs)MALLOC(sizeof(struct dp_pairs)))
25:
26: #define NEXTDPP(r,c) \
27: if(!(r)){NEWDPP(r);(c)=(r);}else{NEWDPP(NEXT(c));(c)=NEXT(c);}
28:
29: struct oEGT eg_up,eg_sp,eg_spm,eg_nf,eg_nfm;
30: struct oEGT eg_znfm,eg_pz,eg_np,eg_mp,eg_ra,eg_mc,eg_gc;
31: int TP,NBP,NMP,NFP,NDP,ZR,NZR;
32:
33: #define NEWDP_pairs ((DP_pairs)MALLOC(sizeof(struct dp_pairs)))
34:
35: extern int (*cmpdl)();
36:
37: void Pdp_gr_flags(), Pdp_gr_print();
38: void Pdp_gr_main(),Pdp_gr_mod_main(),Pdp_f4_main(),Pdp_f4_mod_main();
39:
40: void print_stat(void);
41: void init_stat(void);
42: int dp_load_t(int,DP *);
43: void dp_load(int,DP *);
44: void dp_save(int,DP);
45: void dp_make_flaglist(LIST *);
46: void dp_set_flag(Obj,Obj);
47: int membercheck(NODE,NODE);
48: int gbcheck(NODE);
49: int dl_redble(DL,DL);
50: NODE remove_reducibles(NODE,int);
51: NODE updbase(NODE,int);
52: DP_pairs criterion_F(DP_pairs);
53: int criterion_2(int,int);
54: static DP_pairs collect_pairs_of_hdlcm(DP_pairs,DP_pairs *);
55: DP_pairs criterion_M(DP_pairs);
56: DP_pairs criterion_B(DP_pairs,int);
57: DP_pairs newpairs(NODE,int);
58: DP_pairs updpairs(DP_pairs,NODE,int);
59: void _dp_nf_ptozp(NODE,DP,DP *,int,int,DP *);
60: void _dp_nf_ptozp_mpi(NODE,DP,DP *,int,int,DP *);
61: void _dp_nf(NODE,DP,DP *,int,DP *);
62: void _dp_nf_mod(NODE,DP,DP *,int,int,DP *);
63: void interreduce_mod(NODE,int,int);
64: NODE gb_mod(NODE,int);
65: NODE gbd(NODE,int,NODE,NODE);
66: NODE gb(NODE,int,NODE);
67: NODE gb_f4(NODE);
68: NODE gb_f4_mod(NODE,int);
69: DP_pairs minp(DP_pairs, DP_pairs *);
70: void minsugar(DP_pairs,DP_pairs *,DP_pairs *);
71: NODE append_one(NODE,int);
72: void reducebase_dehomo(NODE,NODE *);
73: int newps_mod(DP,int);
74: int newps_nosave(DP,int,NODE);
75: int newps(DP,int,NODE);
76: void reduceall_mod(NODE,int,NODE *);
77: void reduceall(NODE,NODE *);
78: NODE NODE_sortbi(NODE,int);
79: NODE NODE_sortbi_insert(int, NODE,int);
80: NODE NODE_sortb(NODE,int);
81: NODE NODE_sortb_insert(DP,NODE,int);
82: void prim_part(DP,int,DP *);
83: void setup_arrays(NODE,int,NODE *);
84: int validhc(P,int,NODE);
85: void vlminus(VL,VL,VL *);
86: void printsubst(NODE);
87: void makesubst(VL,NODE *);
88: void pltovl(LIST,VL *);
89: void printdl(DL);
90: int DPPlength(DP_pairs);
91: void dp_gr_mod_main(LIST,LIST,Num,int,struct order_spec *,LIST *);
92: void dp_gr_main(LIST,LIST,Num,int,struct order_spec *,LIST *);
93: void dp_f4_main(LIST,LIST,struct order_spec *,LIST *);
94: void dp_f4_mod_main(LIST,LIST,int,struct order_spec *,LIST *);
95: double get_rtime();
96: void _dpmod_to_vect(DP,DL *,int *);
97: void dp_to_vect(DP,DL *,Q *);
1.2 noro 98: NODE dp_dllist(DP f),symb_merge(NODE,NODE,int);
1.1 noro 99:
100: extern int dp_nelim;
101: extern int dp_fcoeffs;
102: static DP *ps,*psm;
103: static DL *psh;
104: static P *psc;
105:
106: static int *pss;
107: static int psn,pslen;
108: static int NVars,CNVars,PCoeffs;
109: static VL VC;
110: static int NoSugar = 0;
111: static int NoCriB = 0;
112: static int NoGC = 0;
113: static int NoMC = 0;
114: static int NoRA = 0;
115: int Print = 0;
116: static int PrintShort = 0;
117: static int ShowMag = 0;
118: static int Stat = 0;
119: static int Multiple = 0;
120: static int Denominator = 1;
121: static int Top = 0;
122: static int Reverse = 0;
123: static int InterReduce = 0;
124: static int Max_mag = 0;
125: static char *Demand = 0;
126: static int PtozpRA = 0;
127: LIST Dist = 0;
128: int NoGCD = 0;
129: int doing_f4;
130:
131: void Pdp_gr_flags(arg,rp)
132: NODE arg;
133: LIST *rp;
134: {
135: Obj name,value;
136: NODE n;
137:
138: if ( arg ) {
139: asir_assert(ARG0(arg),O_LIST,"dp_gr_flags");
140: n = BDY((LIST)ARG0(arg));
141: while ( n ) {
142: name = (Obj)BDY(n); n = NEXT(n);
143: if ( !n )
144: break;
145: else {
146: value = (Obj)BDY(n); n = NEXT(n);
147: }
148: dp_set_flag(name,value);
149: }
150: }
151: dp_make_flaglist(rp);
152: }
153:
154: void Pdp_gr_print(arg,rp)
155: NODE arg;
156: Q *rp;
157: {
158: Q q;
159:
160: if ( arg ) {
161: asir_assert(ARG0(arg),O_N,"dp_gr_print");
162: q = (Q)ARG0(arg); Print = QTOS(q);
163: } else
164: STOQ(Print,q);
165: *rp = q;
166: }
167:
168: void Pdp_gr_main(arg,rp)
169: NODE arg;
170: LIST *rp;
171: {
172: LIST f,v;
173: Num homo;
174: Q m;
175: int modular;
176: struct order_spec ord;
177:
178: asir_assert(ARG0(arg),O_LIST,"dp_gr_main");
179: asir_assert(ARG1(arg),O_LIST,"dp_gr_main");
180: asir_assert(ARG2(arg),O_N,"dp_gr_main");
181: asir_assert(ARG3(arg),O_N,"dp_gr_main");
182: f = (LIST)ARG0(arg); v = (LIST)ARG1(arg);
183: homo = (Num)ARG2(arg);
184: m = (Q)ARG3(arg);
185: if ( !m )
186: modular = 0;
187: else if ( PL(NM(m))>1 || (PL(NM(m)) == 1 && BD(NM(m))[0] >= 0x80000000) )
188: error("dp_gr_main : too large modulus");
189: else
190: modular = QTOS(m);
191: create_order_spec(ARG4(arg),&ord);
192: dp_gr_main(f,v,homo,modular,&ord,rp);
193: }
194:
195: void Pdp_f4_main(arg,rp)
196: NODE arg;
197: LIST *rp;
198: {
199: LIST f,v;
200: struct order_spec ord;
201:
202: asir_assert(ARG0(arg),O_LIST,"dp_f4_main");
203: asir_assert(ARG1(arg),O_LIST,"dp_f4_main");
204: f = (LIST)ARG0(arg); v = (LIST)ARG1(arg);
205: create_order_spec(ARG2(arg),&ord);
206: dp_f4_main(f,v,&ord,rp);
207: }
208:
209: void Pdp_f4_mod_main(arg,rp)
210: NODE arg;
211: LIST *rp;
212: {
213: LIST f,v;
214: int m;
215: struct order_spec ord;
216:
217: asir_assert(ARG0(arg),O_LIST,"dp_f4_main");
218: asir_assert(ARG1(arg),O_LIST,"dp_f4_main");
219: asir_assert(ARG2(arg),O_N,"dp_f4_main");
220: f = (LIST)ARG0(arg); v = (LIST)ARG1(arg); m = QTOS((Q)ARG2(arg));
221: create_order_spec(ARG3(arg),&ord);
222: dp_f4_mod_main(f,v,m,&ord,rp);
223: }
224:
225: void Pdp_gr_mod_main(arg,rp)
226: NODE arg;
227: LIST *rp;
228: {
229: LIST f,v;
230: Num homo;
231: int m;
232: struct order_spec ord;
233:
234: asir_assert(ARG0(arg),O_LIST,"dp_gr_mod_main");
235: asir_assert(ARG1(arg),O_LIST,"dp_gr_mod_main");
236: asir_assert(ARG2(arg),O_N,"dp_gr_mod_main");
237: asir_assert(ARG3(arg),O_N,"dp_gr_mod_main");
238: f = (LIST)ARG0(arg); v = (LIST)ARG1(arg);
239: homo = (Num)ARG2(arg); m = QTOS((Q)ARG3(arg));
240: create_order_spec(ARG4(arg),&ord);
241: dp_gr_mod_main(f,v,homo,m,&ord,rp);
242: }
243:
244: int eqdl(nv,dl1,dl2)
245: int nv;
246: DL dl1,dl2;
247: {
248: int i;
249: int *b1,*b2;
250:
251: if ( dl1->td != dl2->td )
252: return 0;
253: for ( i = 0, b1 = dl1->d, b2 = dl2->d; i < nv; i++ )
254: if ( b1[i] != b2[i] )
255: break;
256: if ( i == nv )
257: return 1;
258: else
259: return 0;
260: }
261:
262: void _dpmod_to_vect(f,at,b)
263: DP f;
264: DL *at;
265: int *b;
266: {
267: int i,nv;
268: MP m;
269:
270: nv = f->nv;
271: for ( m = BDY(f), i = 0; m; m = NEXT(m), i++ ) {
272: for ( ; !eqdl(nv,m->dl,at[i]); i++ );
273: b[i] = ITOS(m->c);
274: }
275: }
276:
277: void dp_to_vect(f,at,b)
278: DP f;
279: DL *at;
280: Q *b;
281: {
282: int i,nv;
283: MP m;
284:
285: nv = f->nv;
286: for ( m = BDY(f), i = 0; m; m = NEXT(m), i++ ) {
287: for ( ; !eqdl(nv,m->dl,at[i]); i++ );
288: b[i] =(Q)m->c;
289: }
290: }
291:
1.2 noro 292: NODE dp_dllist(f)
1.1 noro 293: DP f;
294: {
1.2 noro 295: MP m;
296: NODE mp,mp0;
1.1 noro 297:
1.2 noro 298: if ( !f )
299: return 0;
1.1 noro 300: mp0 = 0;
301: for ( m = BDY(f); m; m = NEXT(m) ) {
1.2 noro 302: NEXTNODE(mp0,mp); BDY(mp) = (pointer)m->dl;
1.1 noro 303: }
304: NEXT(mp) = 0;
1.2 noro 305: return mp0;
306: }
307:
308: void pdl(f)
309: NODE f;
310: {
311: while ( f ) {
312: printdl(BDY(f)); f = NEXT(f);
313: }
314: fflush(stdout);
315: printf("\n");
1.1 noro 316: }
317:
318: void dp_gr_main(f,v,homo,modular,ord,rp)
319: LIST f,v;
320: Num homo;
321: int modular;
322: struct order_spec *ord;
323: LIST *rp;
324: {
325: int i,mindex,m,nochk;
326: struct order_spec ord1;
327: VL fv,vv,vc;
328: NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
329:
330: mindex = 0; nochk = 0; dp_fcoeffs = 0;
331: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
332: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
333: CNVars = homo ? NVars+1 : NVars;
334: if ( ord->id && NVars != ord->nv )
335: error("dp_gr_main : invalid order specification");
336: initd(ord);
337: if ( homo ) {
338: homogenize_order(ord,NVars,&ord1);
339: for ( fd0 = fi0 = 0, t = BDY(f); t; t = NEXT(t) ) {
340: NEXTNODE(fd0,fd); NEXTNODE(fi0,fi);
341: ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fi)); dp_homo((DP)BDY(fi),(DP *)&BDY(fd));
342: }
343: if ( fd0 ) NEXT(fd) = 0;
344: if ( fi0 ) NEXT(fi) = 0;
345: initd(&ord1);
346: } else {
347: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
348: NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
349: }
350: if ( fd0 ) NEXT(fd) = 0;
351: fi0 = fd0;
352: }
353: if ( modular < 0 ) {
354: modular = -modular; nochk = 1;
355: }
356: if ( modular )
357: m = modular > 1 ? modular : lprime[mindex];
358: else
359: m = 0;
360: makesubst(vc,&subst);
361: setup_arrays(fd0,0,&s);
362: init_stat();
363: while ( 1 ) {
364: if ( homo ) {
365: initd(&ord1); CNVars = NVars+1;
366: }
367: if ( Print && modular ) {
368: fprintf(asir_out,"mod= %d, eval = ",m); printsubst(subst);
369: }
370: x = gb(s,m,subst);
371: if ( x ) {
372: if ( homo ) {
373: reducebase_dehomo(x,&xx); x = xx;
374: initd(ord); CNVars = NVars;
375: }
376: reduceall(x,&xx); x = xx;
377: if ( modular ) {
378: if ( nochk || (membercheck(fi0,x) && gbcheck(x)) )
379: break;
380: } else
381: break;
382: }
383: if ( modular )
384: if ( modular > 1 ) {
385: *rp = 0; return;
386: } else
387: m = lprime[++mindex];
388: makesubst(vc,&subst);
389: psn = length(s);
390: for ( i = psn; i < pslen; i++ ) {
391: pss[i] = 0; psh[i] = 0; psc[i] = 0; ps[i] = 0;
392: }
393: }
394: for ( r0 = 0; x; x = NEXT(x) ) {
395: NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
396: dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
397: }
398: if ( r0 ) NEXT(r) = 0;
399: MKLIST(*rp,r0);
400: print_stat();
401: if ( ShowMag )
402: fprintf(asir_out,"\nMax_mag=%d\n",Max_mag);
403: }
404:
405: void dp_gr_mod_main(f,v,homo,m,ord,rp)
406: LIST f,v;
407: Num homo;
408: int m;
409: struct order_spec *ord;
410: LIST *rp;
411: {
412: struct order_spec ord1;
413: VL fv,vv,vc;
414: NODE fd,fd0,r,r0,t,x,s,xx;
415: DP a,b,c;
416:
417: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
418: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
419: CNVars = homo ? NVars+1 : NVars;
420: if ( ord->id && NVars != ord->nv )
421: error("dp_gr_mod_main : invalid order specification");
422: initd(ord);
423: if ( homo ) {
424: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
425: ptod(CO,vv,(P)BDY(t),&a); dp_homo(a,&b);
426: if ( PCoeffs )
427: dp_mod(b,m,0,&c);
428: else
429: _dp_mod(b,m,(NODE)0,&c);
430: if ( c ) {
431: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
432: }
433: }
434: homogenize_order(ord,NVars,&ord1); initd(&ord1);
435: } else {
436: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
437: ptod(CO,vv,(P)BDY(t),&b);
438: if ( PCoeffs )
439: dp_mod(b,m,0,&c);
440: else
441: _dp_mod(b,m,0,&c);
442: if ( c ) {
443: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
444: }
445: }
446: }
447: if ( fd0 ) NEXT(fd) = 0;
448: setup_arrays(fd0,m,&s);
449: init_stat();
450: if ( homo ) {
451: initd(&ord1); CNVars = NVars+1;
452: }
453: x = gb_mod(s,m);
454: if ( homo ) {
455: reducebase_dehomo(x,&xx); x = xx;
456: initd(ord); CNVars = NVars;
457: }
458: reduceall_mod(x,m,&xx); x = xx;
459: if ( PCoeffs )
460: for ( r0 = 0; x; x = NEXT(x) ) {
461: NEXTNODE(r0,r); mdtop(CO,m,vv,ps[(int)BDY(x)],(P *)&BDY(r));
462: }
463: else
464: for ( r0 = 0; x; x = NEXT(x) ) {
465: NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
466: }
467: print_stat();
468: if ( r0 ) NEXT(r) = 0;
469: MKLIST(*rp,r0);
470: }
471:
472: void dp_f4_main(f,v,ord,rp)
473: LIST f,v;
474: struct order_spec *ord;
475: LIST *rp;
476: {
477: int i,mindex,m,nochk;
478: struct order_spec ord1;
479: VL fv,vv,vc;
480: NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
481:
482: dp_fcoeffs = 0;
483: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
484: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
485: CNVars = NVars;
486: if ( ord->id && NVars != ord->nv )
487: error("dp_f4_main : invalid order specification");
488: initd(ord);
489: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
490: NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
491: }
492: if ( fd0 ) NEXT(fd) = 0;
493: setup_arrays(fd0,0,&s);
494: x = gb_f4(s);
495: reduceall(x,&xx); x = xx;
496: for ( r0 = 0; x; x = NEXT(x) ) {
497: NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
498: dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
499: }
500: if ( r0 ) NEXT(r) = 0;
501: MKLIST(*rp,r0);
502: }
503:
504: void dp_f4_mod_main(f,v,m,ord,rp)
505: LIST f,v;
506: int m;
507: struct order_spec *ord;
508: LIST *rp;
509: {
510: int i;
511: struct order_spec ord1;
512: VL fv,vv,vc;
513: DP b,c;
514: NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
515:
516: dp_fcoeffs = 0;
517: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
518: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
519: CNVars = NVars;
520: if ( ord->id && NVars != ord->nv )
521: error("dp_f4_mod_main : invalid order specification");
522: initd(ord);
523: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
524: ptod(CO,vv,(P)BDY(t),&b);
525: _dp_mod(b,m,0,&c);
526: if ( c ) {
527: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
528: }
529: }
530: if ( fd0 ) NEXT(fd) = 0;
531: setup_arrays(fd0,m,&s);
532: x = gb_f4_mod(s,m);
533: reduceall_mod(x,m,&xx); x = xx;
534: for ( r0 = 0; x; x = NEXT(x) ) {
535: NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
536: }
537: if ( r0 ) NEXT(r) = 0;
538: MKLIST(*rp,r0);
539: }
540:
541: NODE gb_f4(f)
542: NODE f;
543: {
544: int i,j,k,nh,row,col,nv;
545: NODE r,g,gall;
1.2 noro 546: NODE s,s0;
1.1 noro 547: DP_pairs d,dm,dr,t;
1.2 noro 548: DP h,nf,nf1,f1,f2,f21,f21r,sp,sp1,sd,sdm,tdp;
1.1 noro 549: MP mp,mp0;
550: NODE blist,bt,nt;
551: DL *ht,*at;
552: MAT mat,nm;
553: int *colstat;
554: int *rind,*cind;
555: int rank,nred;
556: Q dn;
557: struct oEGT tmp0,tmp1,tmp2,eg_split_symb,eg_split_elim;
558: extern struct oEGT eg_mod,eg_elim,eg_chrem,eg_gschk,eg_intrat,eg_symb;
559:
560: init_eg(&eg_mod); init_eg(&eg_elim); init_eg(&eg_chrem);
561: init_eg(&eg_gschk); init_eg(&eg_intrat); init_eg(&eg_symb);
562:
563: doing_f4 = 1;
564: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
565: i = (int)BDY(r);
566: d = updpairs(d,g,i);
567: g = updbase(g,i);
568: gall = append_one(gall,i);
569: }
570: if ( gall )
571: nv = ((DP)ps[(int)BDY(gall)])->nv;
572: while ( d ) {
573: get_eg(&tmp0);
574: minsugar(d,&dm,&dr); d = dr;
575: if ( Print )
576: fprintf(asir_out,"sugar=%d\n",dm->sugar);
577: blist = 0; s0 = 0;
578: /* asph : sum of all head terms of spoly */
579: for ( t = dm; t; t = NEXT(t) ) {
580: dp_sp(ps[t->dp1],ps[t->dp2],&sp);
581: if ( sp ) {
582: MKNODE(bt,sp,blist); blist = bt;
1.2 noro 583: s0 = symb_merge(s0,dp_dllist(sp),nv);
1.1 noro 584: }
585: }
586: /* s0 : all the terms appeared in symbolic redunction */
1.2 noro 587: for ( s = s0, nred = 0; s; s = NEXT(s) ) {
1.1 noro 588: for ( r = gall; r; r = NEXT(r) )
1.2 noro 589: if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
1.1 noro 590: break;
591: if ( r ) {
1.2 noro 592: dltod(BDY(s),nv,&tdp);
593: dp_subd(tdp,ps[(int)BDY(r)],&sd);
1.4 ! noro 594: muld(CO,ps[(int)BDY(r)],sd,&f2);
1.1 noro 595: MKNODE(bt,f2,blist); blist = bt;
1.2 noro 596: s = symb_merge(s,dp_dllist(f2),nv);
1.1 noro 597: nred++;
598: }
599: }
600:
601: /* the first nred polys in blist are reducers */
602: /* row = the number of all the polys */
603: for ( r = blist, row = 0; r; r = NEXT(r), row++ );
604: ht = (DL *)MALLOC(nred*sizeof(DL));
605: for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
606: ht[i] = BDY((DP)BDY(r))->dl;
1.2 noro 607: for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.1 noro 608: at = (DL *)MALLOC(col*sizeof(DL));
1.2 noro 609: for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
610: at[i] = (DL)BDY(s);
1.1 noro 611: MKMAT(mat,row,col);
612: for ( i = 0, r = blist; i < row; r = NEXT(r), i++ )
613: dp_to_vect(BDY(r),at,(Q *)mat->body[i]);
614: get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
615: init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
616: if ( Print ) {
617: print_eg("Symb",&eg_split_symb);
618: fprintf(asir_out,"mat : %d x %d",row,col);
619: fflush(asir_out);
620: }
1.3 noro 621: #if 0
622: rank = generic_gauss_elim_hensel(mat,&nm,&dn,&rind,&cind);
623: #else
1.1 noro 624: rank = generic_gauss_elim(mat,&nm,&dn,&rind,&cind);
1.3 noro 625: #endif
1.1 noro 626: if ( Print )
627: fprintf(asir_out,"done rank = %d\n",rank,row,col);
628: for ( i = 0; i < rank; i++ ) {
629: for ( k = 0; k < nred; k++ )
630: if ( !cmpdl(nv,at[rind[i]],ht[k]) )
631: break;
632: if ( k == nred ) {
633: /* this is a new base */
634: mp0 = 0;
635: NEXTMP(mp0,mp); mp->dl = at[rind[i]]; mp->c = (P)dn;
636: for ( k = 0; k < col-rank; k++ )
637: if ( nm->body[i][k] ) {
638: NEXTMP(mp0,mp); mp->dl = at[cind[k]];
639: mp->c = (P)nm->body[i][k];
640: }
641: NEXT(mp) = 0;
642: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
643: dp_ptozp(nf,&nf1);
644: nh = newps(nf1,0,0);
645: d = updpairs(d,g,nh);
646: g = updbase(g,nh);
647: gall = append_one(gall,nh);
648: }
649: }
650: }
651: if ( Print ) {
652: print_eg("Symb",&eg_symb);
653: print_eg("Mod",&eg_mod); print_eg("GaussElim",&eg_elim);
654: print_eg("ChRem",&eg_chrem); print_eg("IntToRat",&eg_intrat);
655: print_eg("Check",&eg_gschk);
656: }
657: return g;
658: }
659:
660: NODE gb_f4_mod(f,m)
661: NODE f;
662: int m;
663: {
664: int i,j,k,nh,row,col,nv;
665: NODE r,g,gall;
1.2 noro 666: NODE s,s0;
1.1 noro 667: DP_pairs d,dm,dr,t;
1.2 noro 668: DP h,nf,f1,f2,f21,f21r,sp,sp1,sd,sdm,tdp;
1.1 noro 669: MP mp,mp0;
670: NODE blist,bt,nt;
671: DL *ht,*at;
672: int **mat;
673: int *colstat;
674: int rank,nred,nonzero;
675: struct oEGT tmp0,tmp1,tmp2,eg_split_symb,eg_split_elim;
676: extern struct oEGT eg_symb,eg_elim;
677:
678: init_eg(&eg_symb); init_eg(&eg_elim);
679: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
680: i = (int)BDY(r);
681: d = updpairs(d,g,i);
682: g = updbase(g,i);
683: gall = append_one(gall,i);
684: }
685: if ( gall )
686: nv = ((DP)ps[(int)BDY(gall)])->nv;
687: while ( d ) {
688: get_eg(&tmp0);
689: minsugar(d,&dm,&dr); d = dr;
690: if ( Print )
691: fprintf(asir_out,"sugar=%d\n",dm->sugar);
692: blist = 0; s0 = 0;
693: /* asph : sum of all head terms of spoly */
694: for ( t = dm; t; t = NEXT(t) ) {
695: _dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp);
696: if ( sp ) {
697: MKNODE(bt,sp,blist); blist = bt;
1.2 noro 698: s0 = symb_merge(s0,dp_dllist(sp),nv);
1.1 noro 699: }
700: }
701: /* s0 : all the terms appeared in symbolic redunction */
1.2 noro 702: for ( s = s0, nred = 0; s; s = NEXT(s) ) {
1.1 noro 703: for ( r = gall; r; r = NEXT(r) )
1.2 noro 704: if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
1.1 noro 705: break;
706: if ( r ) {
1.2 noro 707: dltod(BDY(s),nv,&tdp);
708: dp_subd(tdp,ps[(int)BDY(r)],&sd);
1.1 noro 709: _dp_mod(sd,m,0,&sdm);
1.4 ! noro 710: _mulmd(CO,m,ps[(int)BDY(r)],sdm,&f2);
1.1 noro 711: MKNODE(bt,f2,blist); blist = bt;
1.2 noro 712: s = symb_merge(s,dp_dllist(f2),nv);
1.1 noro 713: nred++;
714: }
715: }
716:
717: /* the first nred polys in blist are reducers */
718: /* row = the number of all the polys */
719: for ( r = blist, row = 0; r; r = NEXT(r), row++ );
720: ht = (DL *)MALLOC(nred*sizeof(DL));
721: for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
722: ht[i] = BDY((DP)BDY(r))->dl;
1.2 noro 723: for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.1 noro 724: at = (DL *)MALLOC(col*sizeof(DL));
1.2 noro 725: for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
726: at[i] = (DL)BDY(s);
1.1 noro 727: mat = almat(row,col);
728: for ( i = 0, r = blist; i < row; r = NEXT(r), i++ )
729: _dpmod_to_vect(BDY(r),at,mat[i]);
730: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
731: for ( i = 0, nonzero=0; i < row; i++ )
732: for ( j = 0; j < col; j++ )
733: if ( mat[i][j] )
734: nonzero++;
735: if ( Print )
736: fprintf(asir_out,"mat : %d x %d (nonzero=%f%%)...",
737: row,col,((double)nonzero*100)/(row*col));
738: get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
739: init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
740: rank = generic_gauss_elim_mod(mat,row,col,m,colstat);
741: get_eg(&tmp2); add_eg(&eg_elim,&tmp1,&tmp2);
742: init_eg(&eg_split_elim); add_eg(&eg_split_elim,&tmp1,&tmp2);
743: if ( Print ) {
744: fprintf(asir_out,"done rank = %d\n",rank,row,col);
745: print_eg("Symb",&eg_split_symb);
746: print_eg("Elim",&eg_split_elim);
747: fprintf(asir_out,"\n");
748: }
749: for ( j = 0, i = 0; j < col; j++ )
750: if ( colstat[j] ) {
751: for ( k = 0; k < nred; k++ )
752: if ( !cmpdl(nv,at[j],ht[k]) )
753: break;
754: if ( k == nred ) {
755: /* this is a new base */
756: mp0 = 0;
757: NEXTMP(mp0,mp); mp->dl = at[j]; mp->c = STOI(1);
758: for ( k = j+1; k < col; k++ )
759: if ( !colstat[k] && mat[i][k] ) {
760: NEXTMP(mp0,mp); mp->dl = at[k];
761: mp->c = STOI(mat[i][k]);
762: }
763: NEXT(mp) = 0;
764: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
765: nh = newps_mod(nf,m);
766: d = updpairs(d,g,nh);
767: g = updbase(g,nh);
768: gall = append_one(gall,nh);
769: }
770: i++;
771: }
772: }
773: if ( Print ) {
774: print_eg("Symb",&eg_symb);
775: print_eg("Elim",&eg_elim);
776: fflush(asir_out);
777: }
778: return g;
779: }
780:
781: int DPPlength(n)
782: DP_pairs n;
783: {
784: int i;
785:
786: for ( i = 0; n; n = NEXT(n), i++ );
787: return i;
788: }
789:
790: void printdl(dl)
791: DL dl;
792: {
793: int i;
794:
795: fprintf(asir_out,"<<");
796: for ( i = 0; i < CNVars-1; i++ )
797: fprintf(asir_out,"%d,",dl->d[i]);
798: fprintf(asir_out,"%d>>",dl->d[i]);
799: }
800:
801: void pltovl(l,vl)
802: LIST l;
803: VL *vl;
804: {
805: NODE n;
806: VL r,r0;
807:
808: n = BDY(l);
809: for ( r0 = 0; n; n = NEXT(n) ) {
810: NEXTVL(r0,r); r->v = VR((P)BDY(n));
811: }
812: if ( r0 ) NEXT(r) = 0;
813: *vl = r0;
814: }
815:
816: void makesubst(v,s)
817: VL v;
818: NODE *s;
819: {
820: NODE r,r0;
821: Q q;
822: unsigned int n;
823:
824: for ( r0 = 0; v; v = NEXT(v) ) {
825: NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
826: #if defined(_PA_RISC1_1)
827: n = mrand48()&BMASK; UTOQ(n,q);
828: #else
829: n = random(); UTOQ(n,q);
830: #endif
831: NEXTNODE(r0,r); BDY(r) = (pointer)q;
832: }
833: if ( r0 ) NEXT(r) = 0;
834: *s = r0;
835: }
836:
837: void printsubst(s)
838: NODE s;
839: {
840: fputc('[',asir_out);
841: while ( s ) {
842: printv(CO,(V)BDY(s)); s = NEXT(s);
843: fprintf(asir_out,"->%d",QTOS((Q)BDY(s)));
844: if ( NEXT(s) ) {
845: fputc(',',asir_out); s = NEXT(s);
846: } else
847: break;
848: }
849: fprintf(asir_out,"]\n"); return;
850: }
851:
852: void vlminus(v,w,d)
853: VL v,w,*d;
854: {
855: int i,j,n,m;
856: V *va,*wa;
857: V a;
858: VL r,r0;
859: VL t;
860:
861: for ( n = 0, t = v; t; t = NEXT(t), n++ );
862: va = (V *)ALLOCA(n*sizeof(V));
863: for ( i = 0, t = v; t; t = NEXT(t), i++ )
864: va[i] = t->v;
865: for ( m = 0, t = w; t; t = NEXT(t), m++ );
866: wa = (V *)ALLOCA(m*sizeof(V));
867: for ( i = 0, t = w; t; t = NEXT(t), i++ )
868: wa[i] = t->v;
869: for ( i = 0; i < n; i++ ) {
870: a = va[i];
871: for ( j = 0; j < m; j++ )
872: if ( a == wa[j] )
873: break;
874: if ( j < m )
875: va[i] = 0;
876: }
877: for ( r0 = 0, i = 0; i < n; i++ )
878: if ( va[i] ) { NEXTVL(r0,r); r->v = va[i]; }
879: if ( r0 ) NEXT(r) = 0;
880: *d = r0;
881: }
882:
883: int validhc(a,m,s)
884: P a;
885: int m;
886: NODE s;
887: {
888: P c,c1;
889: V v;
890:
891: if ( !a )
892: return 0;
893: for ( c = a; s; s = NEXT(s) ) {
894: v = (V)BDY(s); s = NEXT(s);
895: substp(CO,c,v,(P)BDY(s),&c1); c = c1;
896: }
897: ptomp(m,c,&c1);
898: return c1 ? 1 : 0;
899: }
900:
901: void setup_arrays(f,m,r)
902: NODE f,*r;
903: int m;
904: {
905: int i;
906: NODE s,s0;
907:
908: f = NODE_sortb(f,1);
909: psn = length(f); pslen = 2*psn;
910: ps = (DP *)MALLOC(pslen*sizeof(DP));
911: psh = (DL *)MALLOC(pslen*sizeof(DL));
912: pss = (int *)MALLOC(pslen*sizeof(int));
913: psc = (P *)MALLOC(pslen*sizeof(P));
914: for ( i = 0; i < psn; i++, f = NEXT(f) ) {
915: prim_part((DP)BDY(f),m,&ps[i]);
916: if ( Demand )
917: dp_save(i,ps[i]);
918: psh[i] = BDY(ps[i])->dl;
919: pss[i] = ps[i]->sugar;
920: psc[i] = BDY(ps[i])->c;
921: }
922: for ( s0 = 0, i = 0; i < psn; i++ ) {
923: NEXTNODE(s0,s); BDY(s) = (pointer)i;
924: }
925: if ( s0 ) NEXT(s) = 0;
926: *r = s0;
927: }
928:
929: void prim_part(f,m,r)
930: DP f,*r;
931: int m;
932: {
933: if ( m > 0 ) {
934: if ( PCoeffs )
935: dp_prim_mod(f,m,r);
936: else
937: *r = f;
938: } else {
939: if ( dp_fcoeffs )
940: *r = f;
941: else if ( PCoeffs )
942: dp_prim(f,r);
943: else
944: dp_ptozp(f,r);
945: }
946: }
947:
948: NODE /* of DP */ NODE_sortb_insert( newdp, nd, dec )
949: DP newdp;
950: NODE /* of DP */ nd;
951: int dec;
952: {
953: register NODE last, p;
954: register DL newdl = BDY(newdp)->dl;
955: register int (*cmpfun)() = cmpdl, nv = CNVars;
956: NODE newnd;
957: int sgn = dec ? 1 : -1;
958: MKNODE( newnd, newdp, 0 );
959: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(last))->dl ) > 0 ) {
960: NEXT(newnd) = last;
961: return newnd;
962: }
963: for ( ; p = NEXT(last); last = p )
964: if ( sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(p))->dl ) > 0 ) break;
965: if ( p ) NEXT(NEXT(last) = newnd) = p;
966: else NEXT(last) = newnd;
967: return nd;
968: }
969:
970: NODE NODE_sortb( node, dec )
971: NODE node;
972: int dec;
973: {
974: register NODE nd, ans;
975:
976: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
977: ans = NODE_sortb_insert( (DP) BDY(nd), ans, dec );
978: return ans;
979: }
980:
981: NODE /* of index */ NODE_sortbi_insert( newdpi, nd, dec )
982: int newdpi;
983: NODE /* of index */ nd;
984: int dec;
985: {
986: register NODE last, p;
987: register DL newdl = psh[newdpi];
988: register int (*cmpfun)() = cmpdl, nv = CNVars;
989: NODE newnd;
990: int sgn = dec ? 1 : -1;
991: MKNODE( newnd, newdpi, 0 );
992: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(last)] ) > 0 ) {
993: NEXT(newnd) = last;
994: return newnd;
995: }
996: for ( ; p = NEXT(last); last = p )
997: if ( sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(p)] ) > 0 ) break;
998: if ( p ) NEXT(NEXT(last) = newnd) = p;
999: else NEXT(last) = newnd;
1000: return nd;
1001: }
1002:
1003: NODE NODE_sortbi( node, dec )
1004: NODE node;
1005: int dec;
1006: {
1007: register NODE nd, ans;
1008:
1009: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
1010: ans = NODE_sortbi_insert( (int) BDY(nd), ans, dec );
1011: return ans;
1012: }
1013:
1014: void reduceall(in,h)
1015: NODE in;
1016: NODE *h;
1017: {
1018: NODE r,t,top;
1019: int n,i,j;
1020: int *w;
1021: DP g,g1;
1022: struct oEGT tmp0,tmp1;
1023:
1024: if ( NoRA ) {
1025: *h = in; return;
1026: }
1027: if ( Print || PrintShort ) {
1028: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1029: }
1030: r = NODE_sortbi(in,0);
1031: n = length(r);
1032: w = (int *)ALLOCA(n*sizeof(int));
1033: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1034: w[i] = (int)BDY(t);
1035: for ( i = 0; i < n; i++ ) {
1036: for ( top = 0, j = n-1; j >= 0; j-- )
1037: if ( j != i ) {
1038: MKNODE(t,(pointer)w[j],top); top = t;
1039: }
1040: get_eg(&tmp0);
1041: dp_load(w[i],&ps[w[i]]);
1042: _dp_nf(top,ps[w[i]],ps,1,&g);
1043:
1044: if ( !PtozpRA || !Multiple )
1045: _dp_nf(top,ps[w[i]],ps,1,&g);
1046: else
1047: #if MPI
1048: _dp_nf_ptozp_mpi(top,ps[w[i]],ps,1,Multiple,&g);
1049: #else
1050: _dp_nf_ptozp(top,ps[w[i]],ps,1,Multiple,&g);
1051: #endif
1052: prim_part(g,0,&g1);
1053: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1054: if ( Print || PrintShort ) {
1055: fprintf(asir_out,"."); fflush(asir_out);
1056: }
1057: w[i] = newps(g1,0,(NODE)0);
1058: }
1059: for ( top = 0, j = n-1; j >= 0; j-- ) {
1060: MKNODE(t,(pointer)w[j],top); top = t;
1061: }
1062: *h = top;
1063: if ( Print || PrintShort )
1064: fprintf(asir_out,"\n");
1065: }
1066:
1067: void reduceall_mod(in,m,h)
1068: NODE in;
1069: int m;
1070: NODE *h;
1071: {
1072: NODE r,t,top;
1073: int n,i,j;
1074: int *w;
1075: DP g;
1076: struct oEGT tmp0,tmp1;
1077:
1078: if ( NoRA ) {
1079: *h = in; return;
1080: }
1081: if ( Print || PrintShort ) {
1082: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1083: }
1084: r = NODE_sortbi(in,0);
1085: n = length(r);
1086: w = (int *)ALLOCA(n*sizeof(int));
1087: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1088: w[i] = (int)BDY(t);
1089: for ( i = 0; i < n; i++ ) {
1090: for ( top = 0, j = n-1; j >= 0; j-- )
1091: if ( j != i ) {
1092: MKNODE(t,(pointer)w[j],top); top = t;
1093: }
1094: get_eg(&tmp0);
1095: if ( PCoeffs )
1096: dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1097: else
1098: _dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1099: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1100: if ( Print || PrintShort ) {
1101: fprintf(asir_out,"."); fflush(asir_out);
1102: }
1103: w[i] = newps_mod(g,m);
1104: }
1105: for ( top = 0, j = n-1; j >= 0; j-- ) {
1106: MKNODE(t,(pointer)w[j],top); top = t;
1107: }
1108: *h = top;
1109: if ( Print || PrintShort )
1110: fprintf(asir_out,"\n");
1111: }
1112:
1113: int newps(a,m,subst)
1114: DP a;
1115: int m;
1116: NODE subst;
1117: {
1118: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1119: return -1;
1120: if ( psn == pslen ) {
1121: pslen *= 2;
1122: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1123: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1124: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1125: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1126: if ( m )
1127: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1128: }
1129: if ( Demand ) {
1130: if ( doing_f4 )
1131: ps[psn] = a;
1132: else
1133: ps[psn] = 0;
1134: dp_save(psn,a);
1135: } else
1136: ps[psn] = a;
1137: psh[psn] = BDY(a)->dl;
1138: pss[psn] = a->sugar;
1139: psc[psn] = BDY(a)->c;
1140: if ( m )
1141: _dp_mod(a,m,subst,&psm[psn]);
1142: return psn++;
1143: }
1144:
1145: int newps_nosave(a,m,subst)
1146: DP a;
1147: int m;
1148: NODE subst;
1149: {
1150: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1151: return -1;
1152: if ( psn == pslen ) {
1153: pslen *= 2;
1154: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1155: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1156: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1157: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1158: if ( m )
1159: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1160: }
1161: ps[psn] = 0;
1162: psh[psn] = BDY(a)->dl;
1163: pss[psn] = a->sugar;
1164: psc[psn] = BDY(a)->c;
1165: if ( m )
1166: _dp_mod(a,m,subst,&psm[psn]);
1167: return psn++;
1168: }
1169:
1170: int newps_mod(a,m)
1171: DP a;
1172: int m;
1173: {
1174: if ( psn == pslen ) {
1175: pslen *= 2;
1176: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1177: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1178: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1179: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P)); /* XXX */
1180: }
1181: ps[psn] = a;
1182: psh[psn] = BDY(ps[psn])->dl;
1183: pss[psn] = ps[psn]->sugar;
1184: return psn++;
1185: }
1186:
1187: void reducebase_dehomo(f,g)
1188: NODE f,*g;
1189: {
1190: int n,i,j,k;
1191: int *r;
1192: DL *w,d;
1193: DP u;
1194: NODE t,top;
1195:
1196: n = length(f);
1197: w = (DL *)ALLOCA(n*sizeof(DL));
1198: r = (int *)ALLOCA(n*sizeof(int));
1199: for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {
1200: r[i] = (int)BDY(t); w[i] = psh[r[i]];
1201: }
1202: for ( i = 0; i < n; i++ ) {
1203: for ( j = 0, d = w[i]; j < n; j++ ) {
1204: if ( j != i ) {
1205: for ( k = 0; k < NVars; k++ )
1206: if ( d->d[k] < w[j]->d[k] )
1207: break;
1208: if ( k == NVars )
1209: break;
1210: }
1211: }
1212: if ( j != n )
1213: r[i] = -1;
1214: }
1215: for ( top = 0, i = n-1; i >= 0; i-- )
1216: if ( r[i] >= 0 ) {
1217: dp_load(r[i],&ps[r[i]]); dp_dehomo(ps[r[i]],&u); j = newps(u,0,0);
1218: MKNODE(t,j,top); top = t;
1219: }
1220: *g = top;
1221: }
1222:
1223: NODE append_one(f,n)
1224: NODE f;
1225: int n;
1226: {
1227: NODE t;
1228:
1229: if ( Reverse || !f ) {
1230: MKNODE(t,(pointer)n,f); return t;
1231: } else {
1232: for ( t = f; NEXT(t); t = NEXT(t) );
1233: MKNODE(NEXT(t),(pointer)n,0);
1234: return f;
1235: }
1236: }
1237:
1238: DP_pairs minp( d, prest )
1239: DP_pairs d, *prest;
1240: {
1241: register DP_pairs m, ml, p, l;
1242: register DL lcm;
1243: register int s, nv = CNVars;
1244: register int (*cmpfun)() = cmpdl;
1245:
1246: if ( !(p = NEXT(m = d)) ) {
1247: *prest = p;
1248: NEXT(m) = 0;
1249: return m;
1250: }
1251: for ( lcm = m->lcm, s = m->sugar, ml = 0, l = m; p; p = NEXT(l = p) )
1252: if ( NoSugar ? (*cmpfun)( nv, lcm, p->lcm ) >= 0 :
1253: (s > p->sugar || s == p->sugar && (*cmpfun)( nv, lcm, p->lcm ) >= 0) )
1254: ml = l, lcm = (m = p)->lcm, s = p->sugar;
1255: if ( !ml ) *prest = NEXT(m);
1256: else {
1257: NEXT(ml) = NEXT(m);
1258: *prest = d;
1259: }
1260: NEXT(m) = 0;
1261: return m;
1262: }
1263:
1264: void minsugar(d,dm,dr)
1265: DP_pairs d;
1266: DP_pairs *dm,*dr;
1267: {
1268: int msugar;
1269: DP_pairs t,dm0,dr0,dmt,drt;
1270:
1271: for ( msugar = d->sugar, t = NEXT(d); t; t = NEXT(t) )
1272: if ( t->sugar < msugar )
1273: msugar = t->sugar;
1274: dm0 = 0; dr0 = 0;
1275: for ( t = d; t; t = NEXT(t) ) {
1276: if ( t->sugar == msugar ) {
1277: NEXTDPP(dm0,dmt);
1278: dmt->dp1 = t->dp1; dmt->dp2 = t->dp2;
1279: dmt->lcm = t->lcm; dmt->sugar = t->sugar;
1280: } else {
1281: NEXTDPP(dr0,drt);
1282: drt->dp1 = t->dp1; drt->dp2 = t->dp2;
1283: drt->lcm = t->lcm; drt->sugar = t->sugar;
1284: }
1285: }
1286: if ( dm0 ) NEXT(dmt) = 0;
1287: if ( dr0 ) NEXT(drt) = 0;
1288: *dm = dm0; *dr = dr0;
1289: }
1290:
1291: NODE gb(f,m,subst)
1292: NODE f;
1293: int m;
1294: NODE subst;
1295: {
1296: int i,nh,prev,mag;
1297: NODE r,g,gall;
1298: DP_pairs d,d1;
1299: DP_pairs l;
1300: DP h,nf,nfm,dp1,dp2;
1301: MP mp;
1302: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tsp0,tsp1,tspm0,tspm1,tnp0,tnp1,tmp0,tmp1;
1303: int skip_nf_flag;
1304: double t_0;
1305: static prev_sugar = -1;
1306:
1307: Max_mag = 0;
1308: prev = 1;
1309: doing_f4 = 0;
1310: if ( m ) {
1311: psm = (DP *)MALLOC(pslen*sizeof(DP));
1312: for ( i = 0; i < psn; i++ )
1313: if ( psh[i] && !validhc(psc[i],m,subst) )
1314: return 0;
1315: else
1316: _dp_mod(ps[i],m,subst,&psm[i]);
1317: }
1318: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1319: i = (int)BDY(r);
1320: d = updpairs(d,g,i);
1321: g = updbase(g,i);
1322: gall = append_one(gall,i);
1323: }
1324: while ( d ) {
1325: get_eg(&tmp0);
1326: l = minp(d,&d1); d = d1;
1327: get_eg(&tmp1); add_eg(&eg_mp,&tmp0,&tmp1);
1328: if ( m ) {
1329: get_eg(&tspm0);
1330: _dp_sp_mod(psm[l->dp1],psm[l->dp2],m,&h);
1331: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1);
1332: get_eg(&tnfm0);
1333: _dp_nf_mod(gall,h,psm,m,0,&nfm);
1334: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1335: } else
1336: nfm = (DP)1;
1337: if ( nfm ) {
1338: get_eg(&tsp0);
1339: if ( Demand ) {
1340: if ( dp_load_t(psn,&nf) ) {
1341: skip_nf_flag = 1;
1342: get_eg(&tsp1); add_eg(&eg_nf,&tsp0,&tsp1);
1343: tnf0=tsp0; tnf1=tsp1;
1344: goto skip_nf;
1345: } else {
1346: skip_nf_flag = 0;
1347: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1348: dp_sp(dp1,dp2,&h);
1349: }
1350: } else
1351: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1352: get_eg(&tsp1); add_eg(&eg_sp,&tsp0,&tsp1);
1353: get_eg(&tnf0);
1354: t_0 = get_rtime();
1355: if ( !Multiple )
1356: _dp_nf(gall,h,ps,!Top,&nf);
1357: else
1358: #if MPI
1359: _dp_nf_ptozp_mpi(gall,h,ps,!Top,Multiple,&nf);
1360: #else
1361: _dp_nf_ptozp(gall,h,ps,!Top,Multiple,&nf);
1362: #endif
1363: if ( Print )
1364: fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
1365: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1366: } else
1367: nf = 0;
1368: skip_nf:
1369: if ( nf ) {
1370: NZR++;
1371: get_eg(&tpz0);
1372: prim_part(nf,0,&h);
1373: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1374: get_eg(&tnp0);
1375: if ( Demand && skip_nf_flag )
1376: nh = newps_nosave(h,m,subst);
1377: else
1378: nh = newps(h,m,subst);
1379: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1380: if ( nh < 0 )
1381: return 0;
1382: d = updpairs(d,g,nh);
1383: g = updbase(g,nh);
1384: gall = append_one(gall,nh);
1385: if ( !dp_fcoeffs && ShowMag ) {
1386: for ( mag = 0, mp = BDY(h); mp; mp = NEXT(mp) )
1387: mag += p_mag((P)mp->c);
1388: Max_mag = MAX(Max_mag,mag);
1389: }
1390: if ( Print ) {
1391: if ( !prev )
1392: fprintf(asir_out,"\n");
1393: print_split_e(&tnf0,&tnf1); print_split_e(&tpz0,&tpz1);
1394: printdl(psh[nh]);
1395: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",
1396: l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
1397: pss[nh]);
1398: if ( ShowMag )
1399: fprintf(asir_out,",mag=%d",mag);
1400: fprintf(asir_out,"\n"); fflush(asir_out);
1401: } else if ( PrintShort ) {
1402: fprintf(asir_out,"+"); fflush(asir_out);
1403: }
1404: prev = 1;
1405: } else {
1406: if ( m )
1407: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1408: ZR++;
1409: if ( Print || PrintShort ) {
1410: if ( h && (h->sugar != prev_sugar) ) {
1411: fprintf(asir_out,"[%d]",h->sugar);
1412: prev_sugar = h->sugar;
1413: }
1414: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1415: }
1416: }
1417: }
1418: if ( Print || PrintShort )
1419: fprintf(asir_out,"gb done\n");
1420: return g;
1421: }
1422:
1423: NODE gbd(f,m,subst,dlist)
1424: NODE f;
1425: int m;
1426: NODE subst;
1427: NODE dlist;
1428: {
1429: int i,nh,prev;
1430: NODE r,g,gall;
1431: struct dp_pairs ol;
1432: DP_pairs l;
1433: NODE pair;
1434: DP h,nf,nfm,dp1,dp2;
1435: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tsp0,tsp1,tspm0,tspm1,tnp0,tnp1;
1436:
1437: prev = 1;
1438: l = &ol;
1439: if ( m ) {
1440: psm = (DP *)MALLOC(pslen*sizeof(DP));
1441: for ( i = 0; i < psn; i++ )
1442: if ( psh[i] && !validhc(psc[i],m,subst) )
1443: return 0;
1444: else
1445: _dp_mod(ps[i],m,subst,&psm[i]);
1446: }
1447: for ( gall = g = 0, r = f; r; r = NEXT(r) ) {
1448: i = (int)BDY(r);
1449: g = updbase(g,i);
1450: gall = append_one(gall,i);
1451: }
1452: while ( dlist ) {
1453: pair = BDY((LIST)BDY(dlist)); dlist = NEXT(dlist);
1454: l->dp1 = QTOS((Q)BDY(pair)); pair = NEXT(pair);
1455: l->dp2 = QTOS((Q)BDY(pair));
1456: if ( m ) {
1457: get_eg(&tspm0);
1458: _dp_sp_mod(psm[l->dp1],psm[l->dp2],m,&h);
1459: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1);
1460: get_eg(&tnfm0);
1461: _dp_nf_mod(gall,h,psm,m,0,&nfm);
1462: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1463: } else
1464: nfm = (DP)1;
1465: if ( nfm ) {
1466: get_eg(&tsp0);
1467: if ( Demand ) {
1468: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1469: dp_sp(dp1,dp2,&h);
1470: } else
1471: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1472: get_eg(&tsp1); add_eg(&eg_sp,&tsp0,&tsp1);
1473: get_eg(&tnfm0);
1474: get_eg(&tnf0);
1475: if ( !Multiple )
1476: _dp_nf(gall,h,ps,!Top,&nf);
1477: else
1478: _dp_nf_ptozp(gall,h,ps,!Top,Multiple,&nf);
1479: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1480: } else
1481: nf = 0;
1482: if ( nf ) {
1483: NZR++;
1484: get_eg(&tpz0);
1485: prim_part(nf,0,&h);
1486: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1487: get_eg(&tnp0);
1488: nh = newps(h,m,subst);
1489: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1490: if ( nh < 0 )
1491: return 0;
1492: g = updbase(g,nh);
1493: gall = append_one(gall,nh);
1494: if ( Print ) {
1495: if ( !prev )
1496: fprintf(asir_out,"\n");
1497: print_split_eg(&tnf0,&tnf1); fflush(asir_out);
1498: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",l->dp1,l->dp2,length(g),length(gall),length(dlist),pss[nh]);
1499: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1500: }
1501: prev = 1;
1502: } else {
1503: if ( m )
1504: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1505: ZR++;
1506: if ( Print ) {
1507: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1508: }
1509: }
1510: }
1511: if ( Print )
1512: fprintf(asir_out,"gb done\n");
1513: return g;
1514: }
1515:
1516: NODE gb_mod(f,m)
1517: NODE f;
1518: int m;
1519: {
1520: int i,nh,prev;
1521: NODE r,g,gall;
1522: DP_pairs d,d1;
1523: DP_pairs l;
1524: DP h,nf;
1525: struct oEGT tnfm0,tnfm1,tspm0,tspm1,tmp0,tmp1,tpz0,tpz1;
1526:
1527: prev = 1;
1528: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1529: i = (int)BDY(r);
1530: d = updpairs(d,g,i);
1531: g = updbase(g,i);
1532: gall = append_one(gall,i);
1533: }
1534: while ( d ) {
1535: get_eg(&tmp0);
1536: l = minp(d,&d1); d = d1;
1537: get_eg(&tmp1); add_eg(&eg_mp,&tmp0,&tmp1); get_eg(&tspm0);
1538: if ( PCoeffs ) {
1539: dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1540: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1); get_eg(&tnfm0);
1541: dp_nf_mod(gall,h,ps,m,!Top,&nf);
1542: } else {
1543: _dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1544: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1); get_eg(&tnfm0);
1545: _dp_nf_mod(gall,h,ps,m,!Top,&nf);
1546: }
1547: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1548: if ( nf ) {
1549: NZR++;
1550: get_eg(&tpz0);
1551: prim_part(nf,m,&h);
1552: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1553: nh = newps_mod(h,m);
1554: if ( nh < 0 )
1555: return 0;
1556: d = updpairs(d,g,nh);
1557: g = updbase(g,nh);
1558: if ( InterReduce )
1559: interreduce_mod(g,nh,m);
1560: gall = append_one(gall,nh);
1561: if ( Print ) {
1562: if ( !prev )
1563: fprintf(asir_out,"\n");
1564: print_split_eg(&tnfm0,&tnfm1); fflush(asir_out);
1565: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",l->dp1,l->dp2,length(g),length(gall),DPPlength(d),pss[nh]);
1566: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1567: }
1568: prev = 1;
1569: } else {
1570: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1571: ZR++;
1572: if ( Print ) {
1573: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1574: }
1575: }
1576: }
1577: if ( Print )
1578: fprintf(asir_out,"gb_mod done\n");
1579: return g;
1580: }
1581:
1582: void interreduce_mod(g,n,m)
1583: NODE g;
1584: int n,m;
1585: {
1586: DP nf;
1587: NODE gn,t;
1588: int i;
1589:
1590: MKNODE(gn,(pointer)n,0);
1591: if ( PCoeffs )
1592: for ( t = g; t; t = NEXT(t) ) {
1593: i = (int)BDY(t);
1594: if ( i != n ) {
1595: dp_nf_mod(gn,ps[i],ps,m,1,&nf); ps[i] = nf;
1596: }
1597: }
1598: else
1599: for ( t = g; t; t = NEXT(t) ) {
1600: i = (int)BDY(t);
1601: if ( i != n ) {
1602: _dp_nf_mod(gn,ps[i],ps,m,1,&nf); ps[i] = nf;
1603: }
1604: }
1605: }
1606:
1607: void _dp_nf_mod(b,g,ps,mod,full,rp)
1608: NODE b;
1609: DP g;
1610: DP *ps;
1611: int mod,full;
1612: DP *rp;
1613: {
1614: DP u,p,d,s,t;
1615: NODE l;
1616: MP m,mr;
1617: int sugar,psugar;
1618:
1619: if ( !g ) {
1620: *rp = 0; return;
1621: }
1622: sugar = g->sugar;
1623: for ( d = 0; g; ) {
1624: for ( u = 0, l = b; l; l = NEXT(l) ) {
1625: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1626: _dp_red_mod(g,p,mod,&u);
1627: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1628: sugar = MAX(sugar,psugar);
1629: if ( !u ) {
1630: if ( d )
1631: d->sugar = sugar;
1632: *rp = d; return;
1633: }
1634: break;
1635: }
1636: }
1637: if ( u )
1638: g = u;
1639: else if ( !full ) {
1640: if ( g ) {
1641: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1642: }
1643: *rp = g; return;
1644: } else {
1645: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1646: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1647: _addmd(CO,mod,d,t,&s); d = s;
1648: dp_rest(g,&t); g = t;
1649: }
1650: }
1651: if ( d )
1652: d->sugar = sugar;
1653: *rp = d;
1654: }
1655:
1656: void _dp_nf(b,g,ps,full,rp)
1657: NODE b;
1658: DP g;
1659: DP *ps;
1660: int full;
1661: DP *rp;
1662: {
1663: DP u,p,d,s,t;
1664: P dmy;
1665: NODE l;
1666: MP m,mr;
1667: int sugar,psugar;
1668:
1669: if ( !g ) {
1670: *rp = 0; return;
1671: }
1672: sugar = g->sugar;
1673: for ( d = 0; g; ) {
1674: for ( u = 0, l = b; l; l = NEXT(l) ) {
1675: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
1676: dp_load((int)BDY(l),&p);
1677: dp_red(d,g,p,&t,&u,&dmy);
1678: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1679: sugar = MAX(sugar,psugar);
1680: if ( !u ) {
1681: if ( d )
1682: d->sugar = sugar;
1683: *rp = d; return;
1684: }
1685: d = t;
1686: break;
1687: }
1688: }
1689: if ( u )
1690: g = u;
1691: else if ( !full ) {
1692: if ( g ) {
1693: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1694: }
1695: *rp = g; return;
1696: } else {
1697: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1698: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1699: addd(CO,d,t,&s); d = s;
1700: dp_rest(g,&t); g = t;
1701: }
1702: }
1703: if ( d )
1704: d->sugar = sugar;
1705: *rp = d;
1706: }
1707:
1708: #define SAFENM(q) ((q)?NM(q):0)
1709:
1710: double pz_t_e, pz_t_d, pz_t_d1, pz_t_c, im_t_s, im_t_r;
1711:
1712: void _dp_nf_ptozp(b,g,ps,full,multiple,r)
1713: NODE b;
1714: DP g;
1715: DP *ps;
1716: int full,multiple;
1717: DP *r;
1718: {
1719: DP u,dp,rp,t,t1,t2,red,shift;
1720: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
1721: N gn,tn,cn;
1722: NODE l;
1723: MP m,mr;
1724: int hmag,denom;
1725: int sugar,psugar;
1726: NODE dist;
1727: STRING imul;
1728: int ndist;
1729: int kara_bit;
1730: extern int kara_mag;
1731: double get_rtime();
1732: double t_0,t_00,tt,ttt,t_p,t_m,t_m1,t_m2,t_s,t_g,t_a;
1733:
1734: if ( !g ) {
1735: *r = 0; return;
1736: }
1737: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1738: t_p = t_m = t_m1 = t_m2 = t_s = t_g = t_a = 0;
1739:
1740: denom = Denominator?Denominator:1;
1741: hmag = multiple*HMAG(g)/denom;
1742: kara_bit = kara_mag*27; /* XXX */
1743: if ( Dist ) {
1744: dist = BDY(Dist);
1745: ndist = length(dist);
1746: }
1747: sugar = g->sugar;
1748:
1749: dc = 0; dp = 0; rc = ONE; rp = g;
1750: MKSTR(imul,"dp_imul_index");
1751:
1752: for ( ; rp; ) {
1753: for ( u = 0, l = b; l; l = NEXT(l) ) {
1754: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
1755: t_0 = get_rtime();
1756: dp_load((int)BDY(l),&red);
1757: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
1758: gcdn(NM(hr),NM(hred),&gn);
1759: divsn(NM(hred),gn,&tn); NTOQ(tn,SGN(hred),cr);
1760: divsn(NM(hr),gn,&tn); NTOQ(tn,SGN(hr),cred);
1761: tt = get_rtime()-t_0; t_p += tt;
1762:
1763: t_0 = get_rtime();
1764: dp_subd(rp,red,&shift);
1765: #if INET
1766: if ( Dist && ndist && HMAG(red) > kara_bit ) {
1767: NODE n0,n1,n2,n3;
1768: int i,s,id;
1769: Obj dmy;
1770: Q ind;
1771:
1772: if ( Print )
1773: fprintf(asir_out,"d");
1774: i = (int)BDY(l); STOQ(i,ind);
1775: chsgnp((P)cred,(P *)&mcred);
1776:
1777: MKNODE(n3,ind,0); MKNODE(n2,mcred,n3);
1778: MKNODE(n1,imul,n2); MKNODE(n0,BDY(dist),n1);
1779: Pox_rpc(n0,&dmy);
1780: muldc(CO,rp,(P)cr,&t);
1781: NEXT(n0)=0;
1782: Pox_pop_local(n0,&t1);
1783: } else
1784: #endif
1785: {
1786: /*
1787: if ( Print )
1788: fprintf(asir_out,"l");
1789: */
1790: t_00 = get_rtime();
1791: muldc(CO,rp,(P)cr,&t);
1792: ttt = get_rtime()-t_00; t_m1 += ttt/dp_nt(rp);
1793: t_00 = get_rtime();
1794: chsgnp((P)cred,(P *)&mcred);
1795: muldc(CO,red,(P)mcred,&t1);
1796: ttt = get_rtime()-t_00; t_m2 += ttt/dp_nt(red);
1797: }
1798: t_00 = get_rtime();
1.4 ! noro 1799: muld(CO,t1,shift,&t2);
1.1 noro 1800: addd(CO,t,t2,&u);
1801: tt = get_rtime(); t_m += tt-t_0;
1802: ttt = get_rtime(); t_s += ttt-t_00;
1803:
1804: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
1805: sugar = MAX(sugar,psugar);
1806: if ( !u ) {
1807: if ( dp )
1808: dp->sugar = sugar;
1809: *r = dp;
1810: goto final;
1811: }
1812: break;
1813: }
1814: }
1815: if ( u ) {
1816: if ( HMAG(u) > hmag ) {
1817: t_0 = get_rtime();
1818: #if INET
1819: if ( Dist && HMAG(u) > kara_bit ) {
1820: if ( Print )
1821: fprintf(asir_out,"D");
1822: dp_ptozp_d(dist,ndist,u,&t);
1823: } else {
1824: if ( Print )
1825: fprintf(asir_out,"L");
1826: dp_ptozp_d(0,0,u,&t);
1827: }
1828: #else
1829: dp_ptozp(u,&t); /* XXX */
1830: #endif
1831: tt = get_rtime()-t_0; t_g += tt;
1832: t_0 = get_rtime();
1833: divsn(NM((Q)BDY(u)->c),NM((Q)BDY(t)->c),&cn); NTOQ(cn,1,cont);
1834: if ( !dp_fcoeffs && Print ) {
1835: fprintf(asir_out,"(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1836: fflush(asir_out);
1837: }
1838: mulq(cr,dc,&dcq);
1839: mulq(cont,rc,&rcq);
1840: gcdn(SAFENM(dcq),SAFENM(rcq),&gn);
1841: divsn(SAFENM(dcq),gn,&tn); NTOQ(tn,SGN(dcq),dc);
1842: divsn(SAFENM(rcq),gn,&tn); NTOQ(tn,SGN(rcq),rc);
1843: tt = get_rtime()-t_0; t_a += tt;
1844: rp = t;
1845: hmag = multiple*HMAG(rp)/denom;
1846: } else {
1847: t_0 = get_rtime();
1848: mulq(cr,dc,&dcq); dc = dcq;
1849: tt = get_rtime()-t_0; t_a += tt;
1850: rp = u;
1851: }
1852: } else if ( !full ) {
1853: if ( rp ) {
1854: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
1855: }
1856: *r = rp;
1857: goto final;
1858: } else {
1859: t_0 = get_rtime();
1860: mulq((Q)BDY(rp)->c,rc,&c);
1861: gcdn(SAFENM(dc),SAFENM(c),&gn); NTOQ(gn,1,gcd);
1862: divsn(SAFENM(dc),gn,&tn); NTOQ(tn,SGN(dc),dcq);
1863: divsn(SAFENM(c),gn,&tn); NTOQ(tn,SGN(c),cq);
1864: muldc(CO,dp,(P)dcq,&t1);
1865: m = BDY(rp); NEWMP(mr); mr->dl = m->dl; mr->c = (P)cq;
1866: NEXT(mr) = 0; MKDP(rp->nv,mr,t); t->sugar = mr->dl->td;
1867: addd(CO,t,t1,&dp);
1868: dc = gcd;
1869:
1870: dp_rest(rp,&t); rp = t;
1871: tt = get_rtime()-t_0; t_a += tt;
1872: }
1873: }
1874: if ( dp )
1875: dp->sugar = sugar;
1876: *r = dp;
1877: final:
1878: if ( Print )
1879: fprintf(asir_out,"(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
1880: t_p,t_m,t_m1,t_m2,t_s,
1881: t_g,t_a,
1882: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
1883: }
1884:
1885: void dp_imul_d();
1886: void imulv();
1887:
1888: void _dp_nf_ptozp_mpi(b,g,ps,full,multiple,r)
1889: NODE b;
1890: DP g;
1891: DP *ps;
1892: int full,multiple;
1893: DP *r;
1894: {
1895: Obj dmy;
1896: DP u,dp,rp,t,t1,t2,red,shift;
1897: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
1898: N gn,tn,cn;
1899: NODE l,n0,n1,n2,n3;
1900: MP m,mr;
1901: int i,n;
1902: int hmag,denom;
1903: int sugar,psugar;
1904: NODE dist;
1905: STRING imul;
1906: int ndist;
1907: int kara_bit;
1908: extern int kara_mag;
1909: extern int mpi_mag;
1910: Q ind;
1911: int id;
1912: double get_rtime();
1913: double t_0,t_00,tt,ttt,t_p,t_m,t_m1,t_m2,t_s,t_g,t_a;
1914:
1915: if ( !g ) {
1916: *r = 0; return;
1917: }
1918: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1919: t_p = t_m = t_m1 = t_m2 = t_s = t_g = t_a = 0;
1920:
1921: denom = Denominator?Denominator:1;
1922: hmag = multiple*HMAG(g)/denom;
1923: kara_bit = kara_mag*27; /* XXX */
1924: if ( Dist ) {
1925: dist = BDY(Dist);
1926: ndist = length(dist);
1927: }
1928: sugar = g->sugar;
1929:
1930: dc = 0; dp = 0; rc = ONE; rp = g;
1931: MKSTR(imul,"dp_imul_index");
1932:
1933: for ( ; rp; ) {
1934: for ( u = 0, l = b; l; l = NEXT(l) ) {
1935: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
1936: t_0 = get_rtime();
1937: dp_load((int)BDY(l),&red);
1938: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
1939: gcdn(NM(hr),NM(hred),&gn);
1940: divsn(NM(hred),gn,&tn); NTOQ(tn,SGN(hred),cr);
1941: divsn(NM(hr),gn,&tn); NTOQ(tn,SGN(hr),cred);
1942: tt = get_rtime()-t_0; t_p += tt;
1943:
1944: t_0 = get_rtime();
1945: dp_subd(rp,red,&shift);
1946: t_00 = get_rtime();
1947: if ( Dist && ndist
1948: && HMAG(rp) > mpi_mag && p_mag((P)cr) > mpi_mag ) {
1949: if ( Print ) fprintf(asir_out,"~");
1950: dp_imul_d(rp,cr,&t);
1951: } else {
1952: if ( Print ) fprintf(asir_out,"_");
1953: muldc(CO,rp,(P)cr,&t);
1954: }
1955: ttt = get_rtime()-t_00; t_m1 += ttt/dp_nt(rp);
1956:
1957: t_00 = get_rtime();
1958: chsgnp((P)cred,(P *)&mcred);
1959: if ( Dist && ndist
1960: && HMAG(red) > mpi_mag && p_mag((P)mcred) > mpi_mag ) {
1961: if ( Print ) fprintf(asir_out,"=");
1962: dp_imul_d(red,mcred,&t1);
1963: } else {
1964: if ( Print ) fprintf(asir_out,"_");
1965: muldc(CO,red,(P)mcred,&t1);
1966: }
1967: ttt = get_rtime()-t_00; t_m2 += ttt/dp_nt(red);
1968:
1969: t_00 = get_rtime();
1.4 ! noro 1970: muld(CO,t1,shift,&t2);
1.1 noro 1971: addd(CO,t,t2,&u);
1972: tt = get_rtime(); t_m += tt-t_0;
1973: ttt = get_rtime(); t_s += ttt-t_00;
1974:
1975: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
1976: sugar = MAX(sugar,psugar);
1977: if ( !u ) {
1978: if ( dp )
1979: dp->sugar = sugar;
1980: *r = dp;
1981: goto final;
1982: }
1983: break;
1984: }
1985: }
1986: if ( u ) {
1987: if ( HMAG(u) > hmag ) {
1988: t_0 = get_rtime();
1989: if ( Dist && HMAG(u) > mpi_mag ) {
1990: if ( Print )
1991: fprintf(asir_out,"D");
1992: dp_ptozp_d(dist,ndist,u,&t);
1993: } else {
1994: if ( Print )
1995: fprintf(asir_out,"L");
1996: dp_ptozp_d(0,0,u,&t);
1997: }
1998: tt = get_rtime()-t_0; t_g += tt;
1999: t_0 = get_rtime();
2000: divsn(NM((Q)BDY(u)->c),NM((Q)BDY(t)->c),&cn); NTOQ(cn,1,cont);
2001: if ( Print ) {
2002: fprintf(asir_out,"(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
2003: fflush(asir_out);
2004: }
2005: mulq(cr,dc,&dcq);
2006: mulq(cont,rc,&rcq);
2007: gcdn(SAFENM(dcq),SAFENM(rcq),&gn);
2008: divsn(SAFENM(dcq),gn,&tn); NTOQ(tn,SGN(dcq),dc);
2009: divsn(SAFENM(rcq),gn,&tn); NTOQ(tn,SGN(rcq),rc);
2010: tt = get_rtime()-t_0; t_a += tt;
2011: rp = t;
2012: hmag = multiple*HMAG(rp)/denom;
2013: } else {
2014: t_0 = get_rtime();
2015: mulq(cr,dc,&dcq); dc = dcq;
2016: tt = get_rtime()-t_0; t_a += tt;
2017: rp = u;
2018: }
2019: } else if ( !full ) {
2020: if ( rp ) {
2021: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
2022: }
2023: *r = rp;
2024: goto final;
2025: } else {
2026: t_0 = get_rtime();
2027: mulq((Q)BDY(rp)->c,(Q)rc,(Q *)&c);
2028: gcdn(SAFENM(dc),SAFENM(c),&gn); NTOQ(gn,1,gcd);
2029: divsn(SAFENM(dc),gn,&tn); NTOQ(tn,SGN(dc),dcq);
2030: divsn(SAFENM(c),gn,&tn); NTOQ(tn,SGN(c),cq);
2031: muldc(CO,dp,(P)dcq,&t1);
2032: m = BDY(rp); NEWMP(mr); mr->dl = m->dl; mr->c = (P)cq;
2033: NEXT(mr) = 0; MKDP(rp->nv,mr,t); t->sugar = mr->dl->td;
2034: addd(CO,t,t1,&dp);
2035: dc = gcd;
2036:
2037: dp_rest(rp,&t); rp = t;
2038: tt = get_rtime()-t_0; t_a += tt;
2039: }
2040: }
2041: if ( dp )
2042: dp->sugar = sugar;
2043: *r = dp;
2044: final:
2045: if ( Print )
2046: fprintf(asir_out,"(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
2047: t_p,t_m,t_m1,t_m2,t_s,
2048: t_g,t_a,
2049: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
2050: }
2051:
2052: void dp_imul_d(p,q,rp)
2053: DP p;
2054: Q q;
2055: DP *rp;
2056: {
2057: int nsep,ndist,i,j,k,l,n;
2058: double t0,t1,t2;
2059: Q *s;
2060: pointer *b;
2061: VECT c,cs,ri;
2062: VECT *r;
2063: MP m;
2064: NODE tn,dist,n0,n1,n2;
2065: Obj dmy;
2066: STRING imul;
2067:
2068: extern LIST Dist;
2069:
2070: if ( !p || !q ) {
2071: *rp = 0; return;
2072: }
2073: dist = BDY(Dist);
2074: for ( tn = dist, ndist = 0; tn; tn = NEXT(tn), ndist++ );
2075: nsep = ndist + 1;
2076: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
2077: if ( n <= nsep ) {
2078: muldc(CO,p,(P)q,rp); return;
2079: }
2080: MKSTR(imul,"imulv");
2081: t0 = get_rtime();
2082: dp_dtov(p,&c);
2083: sepvect(c,nsep,&cs);
2084: r = (VECT *)CALLOC(nsep,sizeof(VECT *));
2085: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
2086: n0 = mknode(4,BDY(tn),imul,b[i],q);
2087: Pox_rpc(n0,&dmy);
2088: }
2089: t1 = get_rtime();
2090: im_t_s += t1 - t0;
2091: imulv(b[i],q,&r[i]);
2092: t1 = get_rtime();
2093: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
2094: MKNODE(n0,BDY(tn),0);
2095: Pox_pop_local(n0,&r[i]);
2096: if ( OID(r[i]) == O_ERR ) {
2097: printexpr(CO,(Obj)r[i]);
2098: error("dp_imul_d : aborted");
2099: }
2100: }
2101: t2 = get_rtime();
2102: im_t_r += t2 - t1;
2103: s = (Q *)CALLOC(n,sizeof(Q));
2104: for ( i = j = 0; i < nsep; i++ ) {
2105: for ( k = 0, ri = r[i], l = ri->len; k < l; k++, j++ ) {
2106: s[j] = (Q)BDY(ri)[k];
2107: }
2108: }
2109: dp_vtod(s,p,rp);
2110: }
2111:
2112: void imulv(w,c,rp)
2113: VECT w;
2114: Q c;
2115: VECT *rp;
2116: {
2117: int n,i;
2118: VECT r;
2119:
2120: n = w->len;
2121: MKVECT(r,n); *rp = r;
2122: for ( i = 0; i < n; i++ )
2123: mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
2124: }
2125:
2126: DP_pairs updpairs( d, g, t)
2127: DP_pairs d;
2128: NODE /* of index */ g;
2129: int t;
2130: {
2131: register DP_pairs d1, dd, nd;
2132: struct oEGT tup0,tup1;
2133: int dl,dl1;
2134:
2135: if ( !g ) return d;
2136: get_eg(&tup0);
2137: if ( !NoCriB && d ) {
2138: dl = DPPlength(d);
2139: d = criterion_B( d, t );
2140: dl -= DPPlength(d); NBP += dl;
2141: }
2142: d1 = newpairs( g, t );
2143: if ( NEXT(d1) ) {
2144: dl = DPPlength(d1); TP += dl;
2145: d1 = criterion_M( d1 );
2146: dl1 = DPPlength(d1); NMP += (dl-dl1); dl = dl1;
2147: d1 = criterion_F( d1 );
2148: dl1 = DPPlength(d1); NFP += (dl-dl1); dl = dl1;
2149: } else
2150: dl = 1;
2151: for ( dd = 0; d1; d1 = nd ) {
2152: nd = NEXT(d1);
2153: if ( !criterion_2( d1->dp1, d1->dp2 ) ) {
2154: NEXT(d1) = dd;
2155: dd = d1;
2156: }
2157: }
2158: dl1 = DPPlength(dd); NDP += (dl-dl1);
2159: get_eg(&tup1);
2160: add_eg(&eg_up,&tup0,&tup1);
2161: if ( !(nd = d) ) return dd;
2162: while ( nd = NEXT(d1 = nd) ) ;
2163: NEXT(d1) = dd;
2164: return d;
2165: }
2166:
2167: DP_pairs newpairs( g, t )
2168: NODE /* of index */ g;
2169: register int t;
2170: {
2171: register NODE r;
2172: register DL tdl = psh[t];
2173: register int ts;
2174: register DP_pairs p, last;
2175: int dp;
2176: register DL dl;
2177: register int s;
2178:
2179: ts = pss[t] - tdl->td;
2180: for ( last = 0, r = g; r; r = NEXT(r) ) {
2181: NEXT(p = NEWDP_pairs) = last;
2182: last = p;
2183: dp = p->dp1 = (int)BDY(r); p->dp2 = t;
2184: p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
2185: p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
2186: }
2187: return last;
2188: }
2189:
2190: DP_pairs criterion_B( d, s )
2191: register DP_pairs d;
2192: int s;
2193: {
2194: register DP_pairs dd, p;
2195: register DL tij, t = psh[s], dltmp;
2196:
2197: if ( !d ) return 0;
2198: NEWDL( dltmp, CNVars );
2199: for ( dd = 0; d; d = p ) {
2200: p = NEXT(d),
2201: tij = d->lcm;
2202: if ( tij->td != lcm_of_DL(CNVars, tij, t, dltmp )->td
2203: || !dl_equal(CNVars, tij, dltmp )
2204: || (tij->td == lcm_of_DL(CNVars, psh[d->dp1], t, dltmp )->td
2205: && dl_equal(CNVars, dltmp, tij ))
2206: || (tij->td == lcm_of_DL(CNVars, psh[d->dp2], t, dltmp )->td
2207: && dl_equal(CNVars, dltmp, tij )) ) {
2208: NEXT(d) = dd;
2209: dd = d;
2210: }
2211: }
2212: return dd;
2213: }
2214:
2215: DP_pairs criterion_M( d1 )
2216: DP_pairs d1;
2217: {
2218: register DP_pairs dd, e, d3, d2, p;
2219: register DL itdl, jtdl;
2220: register int itdltd, jtdltd;
2221:
2222: for ( dd = 0, e = d1; e; e = d3 ) {
2223: if ( !(d2 = NEXT(e)) ) {
2224: NEXT(e) = dd;
2225: return e;
2226: }
2227: itdltd = (itdl = e->lcm)->td;
2228: for ( d3 = 0; d2; d2 = p ) {
2229: p = NEXT(d2),
2230: jtdltd = (jtdl = d2->lcm)->td;
2231: if ( jtdltd == itdltd )
2232: if ( dl_equal(CNVars, itdl, jtdl ) ) ;
2233: else if ( dl_redble( jtdl, itdl ) ) continue;
2234: else if ( dl_redble( itdl, jtdl ) ) goto delit;
2235: else ;
2236: else if ( jtdltd > itdltd )
2237: if ( dl_redble( jtdl, itdl ) ) continue;
2238: else ;
2239: else if ( dl_redble( itdl, jtdl ) ) goto delit;
2240: NEXT(d2) = d3;
2241: d3 = d2;
2242: }
2243: NEXT(e) = dd;
2244: dd = e;
2245: continue;
2246: /**/
2247: delit: NEXT(d2) = d3;
2248: d3 = d2;
2249: for ( ; p; p = d2 ) {
2250: d2 = NEXT(p);
2251: NEXT(p) = d3;
2252: d3 = p;
2253: }
2254: }
2255: return dd;
2256: }
2257:
2258: static DP_pairs collect_pairs_of_hdlcm( d1, prest )
2259: DP_pairs d1, *prest;
2260: {
2261: register DP_pairs w, p, r, s;
2262: register DL ti;
2263: register int td;
2264:
2265: td = (ti = (w = d1)->lcm)->td;
2266: s = NEXT(w);
2267: NEXT(w) = 0;
2268: for ( r = 0; s; s = p ) {
2269: p = NEXT(s);
2270: if ( td == s->lcm->td && dl_equal(CNVars, ti, s->lcm ) )
2271: {
2272: NEXT(s) = w;
2273: w = s;
2274: } else {
2275: NEXT(s) = r;
2276: r = s;
2277: }
2278: }
2279: *prest = r;
2280: return w;
2281: }
2282:
2283: int criterion_2( dp1, dp2 )
2284: int dp1, dp2;
2285: {
2286: register int i, *d1, *d2;
2287:
2288: d1 = psh[dp1]->d, d2 = psh[dp2]->d;
2289: for ( i = CNVars; --i >= 0; d1++, d2++ )
2290: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
2291: return 1;
2292: }
2293:
2294: DP_pairs criterion_F( d1 )
2295: DP_pairs d1;
2296: {
2297: DP_pairs rest, head;
2298: register DP_pairs last, p, r, w;
2299: register int s;
2300:
2301: for ( head = last = 0, p = d1; NEXT(p); ) {
2302: s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
2303: while ( w = NEXT(w) )
2304: if ( criterion_2( w->dp1, w->dp2 ) ) {
2305: r = w;
2306: break;
2307: } else if ( w->sugar < s ) s = (r = w)->sugar;
2308: if ( last ) NEXT(last) = r;
2309: else head = r;
2310: NEXT(last = r) = 0;
2311: if ( !(p = rest) ) return head;
2312: }
2313: if ( !last ) return p;
2314: NEXT(last) = p;
2315: return head;
2316: }
2317:
2318: NODE updbase(g,t)
2319: NODE g;
2320: int t;
2321: {
2322: g = remove_reducibles(g,t);
2323: g = append_one(g,t);
2324: return g;
2325: }
2326:
2327: NODE /* of index */ remove_reducibles( nd, newdp )
2328: NODE /* of index */ nd;
2329: int newdp;
2330: {
2331: register DL dl, dln;
2332: register NODE last, p, head;
2333: register int td;
2334:
2335: dl = psh[newdp];
2336: td = dl->td;
2337: for ( head = last = 0, p = nd; p; ) {
2338: dln = psh[(int)BDY(p)];
2339: if ( dln->td >= td && dl_redble( dln, dl ) ) {
2340: p = NEXT(p);
2341: if ( last ) NEXT(last) = p;
2342: } else {
2343: if ( !last ) head = p;
2344: p = NEXT(last = p);
2345: }
2346: }
2347: return head;
2348: }
2349:
2350: int dl_redble( dl1, dl2 )
2351: DL dl1, dl2;
2352: {
2353: register int n, *d1, *d2;
2354:
2355: for ( d1 = dl1->d, d2 = dl2->d, n = CNVars; --n >= 0; d1++, d2++ )
2356: if ( *d1 < *d2 ) return 0;
2357: return 1;
2358: }
2359:
2360: int gbcheck(f)
2361: NODE f;
2362: {
2363: int i;
2364: NODE r,g,gall;
2365: DP_pairs d,l;
2366: DP h,nf,dp1,dp2;
2367: struct oEGT tmp0,tmp1;
2368:
2369: if ( NoGC )
2370: return 1;
2371: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
2372: i = (int)BDY(r);
2373: d = updpairs(d,g,i);
2374: g = updbase(g,i);
2375: gall = append_one(gall,i);
2376: }
2377: if ( Print || PrintShort ) {
2378: fprintf(asir_out,"gbcheck total %d pairs\n",DPPlength(d)); fflush(asir_out);
2379: }
2380: while ( d ) {
2381: l = d; d = NEXT(d);
2382: get_eg(&tmp0);
2383: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2); dp_sp(dp1,dp2,&h);
2384: _dp_nf(gall,h,ps,1,&nf);
2385: get_eg(&tmp1); add_eg(&eg_gc,&tmp0,&tmp1);
2386: if ( Print || PrintShort ) {
2387: fprintf(asir_out,"."); fflush(asir_out);
2388: }
2389: if ( nf )
2390: return 0;
2391: }
2392: if ( Print || PrintShort )
2393: fprintf(asir_out,"\n");
2394: return 1;
2395: }
2396:
2397: int membercheck(f,x)
2398: NODE f,x;
2399: {
2400: DP g;
2401: struct oEGT tmp0,tmp1;
2402:
2403: if ( NoMC )
2404: return 1;
2405: if ( Print || PrintShort ) {
2406: fprintf(asir_out,"membercheck\n"); fflush(asir_out);
2407: }
2408: for ( ; f; f = NEXT(f) ) {
2409: get_eg(&tmp0);
2410: _dp_nf(x,(DP)BDY(f),ps,1,&g);
2411: get_eg(&tmp1); add_eg(&eg_mc,&tmp0,&tmp1);
2412: if ( Print ) {
2413: print_split_eg(&tmp0,&tmp1); fflush(asir_out);
2414: } else if ( PrintShort ) {
2415: fprintf(asir_out,"."); fflush(asir_out);
2416: }
2417: if ( g )
2418: return 0;
2419: }
2420: if ( Print || PrintShort )
2421: fprintf(asir_out,"\n");
2422: return 1;
2423: }
2424:
2425: void dp_set_flag(name,value)
2426: Obj name,value;
2427: {
2428: char *n;
2429: int v;
2430:
2431: if ( OID(name) != O_STR )
2432: return;
2433: n = BDY((STRING)name);
2434: if ( !strcmp(n,"Demand") ) {
2435: Demand = value ? BDY((STRING)value) : 0; return;
2436: }
2437: if ( !strcmp(n,"Dist") ) {
2438: Dist = (LIST)value; return;
2439: }
2440: if ( value && OID(value) != O_N )
2441: return;
2442: v = QTOS((Q)value);
2443: if ( !strcmp(n,"NoSugar") )
2444: NoSugar = v;
2445: else if ( !strcmp(n,"NoCriB") )
2446: NoCriB = v;
2447: else if ( !strcmp(n,"NoGC") )
2448: NoGC = v;
2449: else if ( !strcmp(n,"NoMC") )
2450: NoMC = v;
2451: else if ( !strcmp(n,"NoRA") )
2452: NoRA = v;
2453: else if ( !strcmp(n,"NoGCD") )
2454: NoGCD = v;
2455: else if ( !strcmp(n,"Top") )
2456: Top = v;
2457: else if ( !strcmp(n,"ShowMag") )
2458: ShowMag = v;
2459: else if ( !strcmp(n,"PrintShort") )
2460: PrintShort = v;
2461: else if ( !strcmp(n,"Print") )
2462: Print = v;
2463: else if ( !strcmp(n,"Stat") )
2464: Stat = v;
2465: else if ( !strcmp(n,"Reverse") )
2466: Reverse = v;
2467: else if ( !strcmp(n,"InterReduce") )
2468: InterReduce = v;
2469: else if ( !strcmp(n,"Multiple") )
2470: Multiple = v;
2471: else if ( !strcmp(n,"Denominator") )
2472: Denominator = v;
2473: else if ( !strcmp(n,"PtozpRA") )
2474: PtozpRA = v;
2475: }
2476:
2477: void dp_make_flaglist(list)
2478: LIST *list;
2479: {
2480: Q v;
2481: STRING name,path;
2482: NODE n,n1;
2483:
2484: STOQ(Multiple,v); MKNODE(n,v,0); MKSTR(name,"Multiple"); MKNODE(n1,name,n); n = n1;
2485: STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
2486: MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
2487: STOQ(InterReduce,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"InterReduce"); MKNODE(n1,name,n); n = n1;
2488: STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
2489: STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
2490: STOQ(Print,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Print"); MKNODE(n1,name,n); n = n1;
2491: STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
2492: STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
2493: STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
2494: STOQ(NoGCD,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGCD"); MKNODE(n1,name,n); n = n1;
2495: STOQ(NoRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoRA"); MKNODE(n1,name,n); n = n1;
2496: STOQ(NoMC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoMC"); MKNODE(n1,name,n); n = n1;
2497: STOQ(NoGC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGC"); MKNODE(n1,name,n); n = n1;
2498: STOQ(NoCriB,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoCriB"); MKNODE(n1,name,n); n = n1;
2499: STOQ(NoSugar,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoSugar"); MKNODE(n1,name,n); n = n1;
2500: if ( Demand )
2501: MKSTR(path,Demand);
2502: else
2503: path = 0;
2504: MKNODE(n1,path,n); n = n1; MKSTR(name,"Demand"); MKNODE(n1,name,n); n = n1;
2505: MKLIST(*list,n);
2506: }
2507:
2508: #define DELIM '/'
2509:
2510: #if defined(THINK_C)
2511: #undef DELIM
2512: #define DELIM ':'
2513: #endif
2514:
2515: void dp_save(index,p)
2516: int index;
2517: DP p;
2518: {
2519: FILE *fp;
2520: char path[BUFSIZ];
2521:
2522: sprintf(path,"%s%c%d",Demand,DELIM,index);
2523: #if defined(VISUAL) || defined(THINK_C)
2524: if ( !(fp = fopen(path,"wb") ) )
2525: #else
2526: if ( !(fp = fopen(path,"w") ) )
2527: #endif
2528: error("dp_save : cannot open a file");
2529: savevl(fp,VC); saveobj(fp,(Obj)p); fclose(fp);
2530: }
2531:
2532: void dp_load(index,p)
2533: int index;
2534: DP *p;
2535: {
2536: FILE *fp;
2537: char path[BUFSIZ];
2538:
2539: if ( !Demand || ps[index] )
2540: *p = ps[index];
2541: else {
2542: sprintf(path,"%s%c%d",Demand,DELIM,index);
2543: #if defined(VISUAL) || defined(THINK_C)
2544: if ( !(fp = fopen(path,"rb") ) )
2545: #else
2546: if ( !(fp = fopen(path,"r") ) )
2547: #endif
2548: error("dp_load : cannot open a file");
2549: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp);
2550: }
2551: }
2552:
2553: int dp_load_t(index,p)
2554: int index;
2555: DP *p;
2556: {
2557: FILE *fp;
2558: char path[BUFSIZ];
2559:
2560: sprintf(path,"%s%c%d",Demand,DELIM,index);
2561: #if defined(VISUAL) || defined(THINK_C)
2562: if ( !(fp = fopen(path,"rb") ) )
2563: #else
2564: if ( !(fp = fopen(path,"r") ) )
2565: #endif
2566: return 0;
2567: else {
2568: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp); return 1;
2569: }
2570: }
2571:
2572: void init_stat() {
2573: init_eg(&eg_up); init_eg(&eg_sp); init_eg(&eg_spm);
2574: init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
2575: init_eg(&eg_pz); init_eg(&eg_np); init_eg(&eg_mp);
2576: init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
2577: ZR = NZR = TP = NBP = NMP = NFP = NDP = 0;
2578: }
2579:
2580: void print_stat() {
2581: if ( !Print && !Stat )
2582: return;
2583: print_eg("UP",&eg_up); print_eg("SP",&eg_sp); print_eg("SPM",&eg_spm);
2584: print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
2585: print_eg("PZ",&eg_pz); print_eg("NP",&eg_np); print_eg("MP",&eg_mp);
2586: print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
2587: fprintf(asir_out,"T=%d,B=%d M=%d F=%d D=%d ZR=%d NZR=%d\n",TP,NBP,NMP,NFP,NDP,ZR,NZR);
2588: }
2589:
2590: /*
2591: * Old codes
2592: */
2593:
2594: #if 0
2595: void _dp_nf_ptozp(b,g,ps,full,multiple,rp)
2596: NODE b;
2597: DP g;
2598: DP *ps;
2599: int full,multiple;
2600: DP *rp;
2601: {
2602: DP u,p,d,s,t,dmy;
2603: NODE l;
2604: MP m,mr;
2605: int i,n;
2606: int hmag,denom;
2607: int sugar,psugar;
2608: NODE dist;
2609: int ndist;
2610: int kara_bit;
2611: extern int kara_mag;
2612:
2613: if ( !g ) {
2614: *rp = 0; return;
2615: }
2616: denom = Denominator?Denominator:1;
2617: hmag = multiple*HMAG(g)/denom;
2618: kara_bit = kara_mag*27; /* XXX */
2619: if ( Dist ) {
2620: dist = BDY(Dist);
2621: ndist = length(dist);
2622: }
2623: sugar = g->sugar;
2624: for ( d = 0; g; ) {
2625: for ( u = 0, l = b; l; l = NEXT(l) ) {
2626: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
2627: dp_load((int)BDY(l),&p);
2628: dp_red(d,g,p,&t,&u,&dmy);
2629: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
2630: sugar = MAX(sugar,psugar);
2631: if ( !u ) {
2632: if ( d )
2633: d->sugar = sugar;
2634: *rp = d; return;
2635: }
2636: d = t;
2637: break;
2638: }
2639: }
2640: if ( u ) {
2641: g = u;
2642: if ( d ) {
2643: if ( HMAG(d) > hmag ) {
2644: if ( Dist && HMAG(g) > kara_bit )
2645: dp_ptozp2_d(dist,ndist,d,g,&t,&u);
2646: else
2647: dp_ptozp2(d,g,&t,&u);
2648: d = t; g = u;
2649: hmag = multiple*HMAG(d)/denom;
2650: }
2651: } else {
2652: if ( HMAG(g) > hmag ) {
2653: #if INET
2654: if ( Dist && HMAG(g) > kara_bit ) {
2655: dp_ptozp_d(dist,ndist,g,&t);
2656: } else
2657: #endif
2658: dp_ptozp(g,&t);
2659: g = t;
2660: hmag = multiple*HMAG(g)/denom;
2661: }
2662: }
2663: }
2664: else if ( !full ) {
2665: if ( g ) {
2666: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
2667: }
2668: *rp = g; return;
2669: } else {
2670: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
2671: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
2672: addd(CO,d,t,&s); d = s;
2673: dp_rest(g,&t); g = t;
2674: }
2675: }
2676: if ( d )
2677: d->sugar = sugar;
2678: *rp = d;
2679: }
2680:
2681: int criterion_2(dp1,dp2)
2682: {
2683: DP g1,g2,g,s1,s2;
2684:
2685: monomial_gcd(ps[dp1],&g1); monomial_gcd(ps[dp2],&g2);
2686: dp_gcd(g1,g2,&g);
2687: dp_subd(ps[dp1],g,&s1); dp_subd(ps[dp2],g,&s2);
2688: return _criterion_2(s1,s2);
2689: }
2690:
2691: int _criterion_2( dp1, dp2 )
2692: DP dp1, dp2;
2693: {
2694: register int i, *d1, *d2;
2695:
2696: d1 = BDY(dp1)->dl->d; d2 = BDY(dp2)->dl->d;
2697: for ( i = CNVars; --i >= 0; d1++, d2++ )
2698: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
2699: return 1;
2700: }
2701:
2702: void dp_gcd(p1,p2,rp)
2703: DP p1,p2;
2704: DP *rp;
2705: {
2706: int i,n,td;
2707: DL d1,d2,d;
2708: MP m;
2709: DP s;
2710:
2711: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
2712: NEWDL(d,n);
2713: for ( i = 0, td = 0; i < n; i++ ) {
2714: d->d[i] = MIN(d1->d[i],d2->d[i]);
2715: td += d->d[i];
2716: }
2717: d->td = td;
2718: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
2719: *rp = s;
2720: }
2721:
2722: void monomial_gcd(p,rp)
2723: DP p;
2724: DP *rp;
2725: {
2726: int n,i,td;
2727: DL d,d1;
2728: MP m;
2729: DP s;
2730:
2731: n = p->nv; m = BDY(p); d = m->dl;
2732: NEWDL(d1,n);
2733: for ( i = 0; i < n; i++ )
2734: d1->d[i] = d->d[i];
2735: for ( m = NEXT(m); m; m = NEXT(m) ) {
2736: d = m->dl;
2737: for ( i = 0; i < n; i++ )
2738: d1->d[i] = MIN(d1->d[i],d->d[i]);
2739: }
2740: for ( i = 0, td = 0; i < n; i++ )
2741: td += d1->d[i];
2742: NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
2743: *rp = s;
2744: }
2745: #endif