![[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.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.1.1.1 1999/12/03 07:39:07 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.1 noro 594: muld(CO,sd,ps[(int)BDY(r)],&f2);
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: }
621: rank = generic_gauss_elim(mat,&nm,&dn,&rind,&cind);
622: if ( Print )
623: fprintf(asir_out,"done rank = %d\n",rank,row,col);
624: for ( i = 0; i < rank; i++ ) {
625: for ( k = 0; k < nred; k++ )
626: if ( !cmpdl(nv,at[rind[i]],ht[k]) )
627: break;
628: if ( k == nred ) {
629: /* this is a new base */
630: mp0 = 0;
631: NEXTMP(mp0,mp); mp->dl = at[rind[i]]; mp->c = (P)dn;
632: for ( k = 0; k < col-rank; k++ )
633: if ( nm->body[i][k] ) {
634: NEXTMP(mp0,mp); mp->dl = at[cind[k]];
635: mp->c = (P)nm->body[i][k];
636: }
637: NEXT(mp) = 0;
638: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
639: dp_ptozp(nf,&nf1);
640: nh = newps(nf1,0,0);
641: d = updpairs(d,g,nh);
642: g = updbase(g,nh);
643: gall = append_one(gall,nh);
644: }
645: }
646: }
647: if ( Print ) {
648: print_eg("Symb",&eg_symb);
649: print_eg("Mod",&eg_mod); print_eg("GaussElim",&eg_elim);
650: print_eg("ChRem",&eg_chrem); print_eg("IntToRat",&eg_intrat);
651: print_eg("Check",&eg_gschk);
652: }
653: return g;
654: }
655:
656: NODE gb_f4_mod(f,m)
657: NODE f;
658: int m;
659: {
660: int i,j,k,nh,row,col,nv;
661: NODE r,g,gall;
1.2 ! noro 662: NODE s,s0;
1.1 noro 663: DP_pairs d,dm,dr,t;
1.2 ! noro 664: DP h,nf,f1,f2,f21,f21r,sp,sp1,sd,sdm,tdp;
1.1 noro 665: MP mp,mp0;
666: NODE blist,bt,nt;
667: DL *ht,*at;
668: int **mat;
669: int *colstat;
670: int rank,nred,nonzero;
671: struct oEGT tmp0,tmp1,tmp2,eg_split_symb,eg_split_elim;
672: extern struct oEGT eg_symb,eg_elim;
673:
674: init_eg(&eg_symb); init_eg(&eg_elim);
675: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
676: i = (int)BDY(r);
677: d = updpairs(d,g,i);
678: g = updbase(g,i);
679: gall = append_one(gall,i);
680: }
681: if ( gall )
682: nv = ((DP)ps[(int)BDY(gall)])->nv;
683: while ( d ) {
684: get_eg(&tmp0);
685: minsugar(d,&dm,&dr); d = dr;
686: if ( Print )
687: fprintf(asir_out,"sugar=%d\n",dm->sugar);
688: blist = 0; s0 = 0;
689: /* asph : sum of all head terms of spoly */
690: for ( t = dm; t; t = NEXT(t) ) {
691: _dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp);
692: if ( sp ) {
693: MKNODE(bt,sp,blist); blist = bt;
1.2 ! noro 694: s0 = symb_merge(s0,dp_dllist(sp),nv);
1.1 noro 695: }
696: }
697: /* s0 : all the terms appeared in symbolic redunction */
1.2 ! noro 698: for ( s = s0, nred = 0; s; s = NEXT(s) ) {
1.1 noro 699: for ( r = gall; r; r = NEXT(r) )
1.2 ! noro 700: if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
1.1 noro 701: break;
702: if ( r ) {
1.2 ! noro 703: dltod(BDY(s),nv,&tdp);
! 704: dp_subd(tdp,ps[(int)BDY(r)],&sd);
1.1 noro 705: _dp_mod(sd,m,0,&sdm);
706: _mulmd(CO,m,sdm,ps[(int)BDY(r)],&f2);
707: MKNODE(bt,f2,blist); blist = bt;
1.2 ! noro 708: s = symb_merge(s,dp_dllist(f2),nv);
1.1 noro 709: nred++;
710: }
711: }
712:
713: /* the first nred polys in blist are reducers */
714: /* row = the number of all the polys */
715: for ( r = blist, row = 0; r; r = NEXT(r), row++ );
716: ht = (DL *)MALLOC(nred*sizeof(DL));
717: for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
718: ht[i] = BDY((DP)BDY(r))->dl;
1.2 ! noro 719: for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.1 noro 720: at = (DL *)MALLOC(col*sizeof(DL));
1.2 ! noro 721: for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
! 722: at[i] = (DL)BDY(s);
1.1 noro 723: mat = almat(row,col);
724: for ( i = 0, r = blist; i < row; r = NEXT(r), i++ )
725: _dpmod_to_vect(BDY(r),at,mat[i]);
726: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
727: for ( i = 0, nonzero=0; i < row; i++ )
728: for ( j = 0; j < col; j++ )
729: if ( mat[i][j] )
730: nonzero++;
731: if ( Print )
732: fprintf(asir_out,"mat : %d x %d (nonzero=%f%%)...",
733: row,col,((double)nonzero*100)/(row*col));
734: get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
735: init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
736: rank = generic_gauss_elim_mod(mat,row,col,m,colstat);
737: get_eg(&tmp2); add_eg(&eg_elim,&tmp1,&tmp2);
738: init_eg(&eg_split_elim); add_eg(&eg_split_elim,&tmp1,&tmp2);
739: if ( Print ) {
740: fprintf(asir_out,"done rank = %d\n",rank,row,col);
741: print_eg("Symb",&eg_split_symb);
742: print_eg("Elim",&eg_split_elim);
743: fprintf(asir_out,"\n");
744: }
745: for ( j = 0, i = 0; j < col; j++ )
746: if ( colstat[j] ) {
747: for ( k = 0; k < nred; k++ )
748: if ( !cmpdl(nv,at[j],ht[k]) )
749: break;
750: if ( k == nred ) {
751: /* this is a new base */
752: mp0 = 0;
753: NEXTMP(mp0,mp); mp->dl = at[j]; mp->c = STOI(1);
754: for ( k = j+1; k < col; k++ )
755: if ( !colstat[k] && mat[i][k] ) {
756: NEXTMP(mp0,mp); mp->dl = at[k];
757: mp->c = STOI(mat[i][k]);
758: }
759: NEXT(mp) = 0;
760: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
761: nh = newps_mod(nf,m);
762: d = updpairs(d,g,nh);
763: g = updbase(g,nh);
764: gall = append_one(gall,nh);
765: }
766: i++;
767: }
768: }
769: if ( Print ) {
770: print_eg("Symb",&eg_symb);
771: print_eg("Elim",&eg_elim);
772: fflush(asir_out);
773: }
774: return g;
775: }
776:
777: int DPPlength(n)
778: DP_pairs n;
779: {
780: int i;
781:
782: for ( i = 0; n; n = NEXT(n), i++ );
783: return i;
784: }
785:
786: void printdl(dl)
787: DL dl;
788: {
789: int i;
790:
791: fprintf(asir_out,"<<");
792: for ( i = 0; i < CNVars-1; i++ )
793: fprintf(asir_out,"%d,",dl->d[i]);
794: fprintf(asir_out,"%d>>",dl->d[i]);
795: }
796:
797: void pltovl(l,vl)
798: LIST l;
799: VL *vl;
800: {
801: NODE n;
802: VL r,r0;
803:
804: n = BDY(l);
805: for ( r0 = 0; n; n = NEXT(n) ) {
806: NEXTVL(r0,r); r->v = VR((P)BDY(n));
807: }
808: if ( r0 ) NEXT(r) = 0;
809: *vl = r0;
810: }
811:
812: void makesubst(v,s)
813: VL v;
814: NODE *s;
815: {
816: NODE r,r0;
817: Q q;
818: unsigned int n;
819:
820: for ( r0 = 0; v; v = NEXT(v) ) {
821: NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
822: #if defined(_PA_RISC1_1)
823: n = mrand48()&BMASK; UTOQ(n,q);
824: #else
825: n = random(); UTOQ(n,q);
826: #endif
827: NEXTNODE(r0,r); BDY(r) = (pointer)q;
828: }
829: if ( r0 ) NEXT(r) = 0;
830: *s = r0;
831: }
832:
833: void printsubst(s)
834: NODE s;
835: {
836: fputc('[',asir_out);
837: while ( s ) {
838: printv(CO,(V)BDY(s)); s = NEXT(s);
839: fprintf(asir_out,"->%d",QTOS((Q)BDY(s)));
840: if ( NEXT(s) ) {
841: fputc(',',asir_out); s = NEXT(s);
842: } else
843: break;
844: }
845: fprintf(asir_out,"]\n"); return;
846: }
847:
848: void vlminus(v,w,d)
849: VL v,w,*d;
850: {
851: int i,j,n,m;
852: V *va,*wa;
853: V a;
854: VL r,r0;
855: VL t;
856:
857: for ( n = 0, t = v; t; t = NEXT(t), n++ );
858: va = (V *)ALLOCA(n*sizeof(V));
859: for ( i = 0, t = v; t; t = NEXT(t), i++ )
860: va[i] = t->v;
861: for ( m = 0, t = w; t; t = NEXT(t), m++ );
862: wa = (V *)ALLOCA(m*sizeof(V));
863: for ( i = 0, t = w; t; t = NEXT(t), i++ )
864: wa[i] = t->v;
865: for ( i = 0; i < n; i++ ) {
866: a = va[i];
867: for ( j = 0; j < m; j++ )
868: if ( a == wa[j] )
869: break;
870: if ( j < m )
871: va[i] = 0;
872: }
873: for ( r0 = 0, i = 0; i < n; i++ )
874: if ( va[i] ) { NEXTVL(r0,r); r->v = va[i]; }
875: if ( r0 ) NEXT(r) = 0;
876: *d = r0;
877: }
878:
879: int validhc(a,m,s)
880: P a;
881: int m;
882: NODE s;
883: {
884: P c,c1;
885: V v;
886:
887: if ( !a )
888: return 0;
889: for ( c = a; s; s = NEXT(s) ) {
890: v = (V)BDY(s); s = NEXT(s);
891: substp(CO,c,v,(P)BDY(s),&c1); c = c1;
892: }
893: ptomp(m,c,&c1);
894: return c1 ? 1 : 0;
895: }
896:
897: void setup_arrays(f,m,r)
898: NODE f,*r;
899: int m;
900: {
901: int i;
902: NODE s,s0;
903:
904: f = NODE_sortb(f,1);
905: psn = length(f); pslen = 2*psn;
906: ps = (DP *)MALLOC(pslen*sizeof(DP));
907: psh = (DL *)MALLOC(pslen*sizeof(DL));
908: pss = (int *)MALLOC(pslen*sizeof(int));
909: psc = (P *)MALLOC(pslen*sizeof(P));
910: for ( i = 0; i < psn; i++, f = NEXT(f) ) {
911: prim_part((DP)BDY(f),m,&ps[i]);
912: if ( Demand )
913: dp_save(i,ps[i]);
914: psh[i] = BDY(ps[i])->dl;
915: pss[i] = ps[i]->sugar;
916: psc[i] = BDY(ps[i])->c;
917: }
918: for ( s0 = 0, i = 0; i < psn; i++ ) {
919: NEXTNODE(s0,s); BDY(s) = (pointer)i;
920: }
921: if ( s0 ) NEXT(s) = 0;
922: *r = s0;
923: }
924:
925: void prim_part(f,m,r)
926: DP f,*r;
927: int m;
928: {
929: if ( m > 0 ) {
930: if ( PCoeffs )
931: dp_prim_mod(f,m,r);
932: else
933: *r = f;
934: } else {
935: if ( dp_fcoeffs )
936: *r = f;
937: else if ( PCoeffs )
938: dp_prim(f,r);
939: else
940: dp_ptozp(f,r);
941: }
942: }
943:
944: NODE /* of DP */ NODE_sortb_insert( newdp, nd, dec )
945: DP newdp;
946: NODE /* of DP */ nd;
947: int dec;
948: {
949: register NODE last, p;
950: register DL newdl = BDY(newdp)->dl;
951: register int (*cmpfun)() = cmpdl, nv = CNVars;
952: NODE newnd;
953: int sgn = dec ? 1 : -1;
954: MKNODE( newnd, newdp, 0 );
955: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(last))->dl ) > 0 ) {
956: NEXT(newnd) = last;
957: return newnd;
958: }
959: for ( ; p = NEXT(last); last = p )
960: if ( sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(p))->dl ) > 0 ) break;
961: if ( p ) NEXT(NEXT(last) = newnd) = p;
962: else NEXT(last) = newnd;
963: return nd;
964: }
965:
966: NODE NODE_sortb( node, dec )
967: NODE node;
968: int dec;
969: {
970: register NODE nd, ans;
971:
972: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
973: ans = NODE_sortb_insert( (DP) BDY(nd), ans, dec );
974: return ans;
975: }
976:
977: NODE /* of index */ NODE_sortbi_insert( newdpi, nd, dec )
978: int newdpi;
979: NODE /* of index */ nd;
980: int dec;
981: {
982: register NODE last, p;
983: register DL newdl = psh[newdpi];
984: register int (*cmpfun)() = cmpdl, nv = CNVars;
985: NODE newnd;
986: int sgn = dec ? 1 : -1;
987: MKNODE( newnd, newdpi, 0 );
988: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(last)] ) > 0 ) {
989: NEXT(newnd) = last;
990: return newnd;
991: }
992: for ( ; p = NEXT(last); last = p )
993: if ( sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(p)] ) > 0 ) break;
994: if ( p ) NEXT(NEXT(last) = newnd) = p;
995: else NEXT(last) = newnd;
996: return nd;
997: }
998:
999: NODE NODE_sortbi( node, dec )
1000: NODE node;
1001: int dec;
1002: {
1003: register NODE nd, ans;
1004:
1005: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
1006: ans = NODE_sortbi_insert( (int) BDY(nd), ans, dec );
1007: return ans;
1008: }
1009:
1010: void reduceall(in,h)
1011: NODE in;
1012: NODE *h;
1013: {
1014: NODE r,t,top;
1015: int n,i,j;
1016: int *w;
1017: DP g,g1;
1018: struct oEGT tmp0,tmp1;
1019:
1020: if ( NoRA ) {
1021: *h = in; return;
1022: }
1023: if ( Print || PrintShort ) {
1024: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1025: }
1026: r = NODE_sortbi(in,0);
1027: n = length(r);
1028: w = (int *)ALLOCA(n*sizeof(int));
1029: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1030: w[i] = (int)BDY(t);
1031: for ( i = 0; i < n; i++ ) {
1032: for ( top = 0, j = n-1; j >= 0; j-- )
1033: if ( j != i ) {
1034: MKNODE(t,(pointer)w[j],top); top = t;
1035: }
1036: get_eg(&tmp0);
1037: dp_load(w[i],&ps[w[i]]);
1038: _dp_nf(top,ps[w[i]],ps,1,&g);
1039:
1040: if ( !PtozpRA || !Multiple )
1041: _dp_nf(top,ps[w[i]],ps,1,&g);
1042: else
1043: #if MPI
1044: _dp_nf_ptozp_mpi(top,ps[w[i]],ps,1,Multiple,&g);
1045: #else
1046: _dp_nf_ptozp(top,ps[w[i]],ps,1,Multiple,&g);
1047: #endif
1048: prim_part(g,0,&g1);
1049: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1050: if ( Print || PrintShort ) {
1051: fprintf(asir_out,"."); fflush(asir_out);
1052: }
1053: w[i] = newps(g1,0,(NODE)0);
1054: }
1055: for ( top = 0, j = n-1; j >= 0; j-- ) {
1056: MKNODE(t,(pointer)w[j],top); top = t;
1057: }
1058: *h = top;
1059: if ( Print || PrintShort )
1060: fprintf(asir_out,"\n");
1061: }
1062:
1063: void reduceall_mod(in,m,h)
1064: NODE in;
1065: int m;
1066: NODE *h;
1067: {
1068: NODE r,t,top;
1069: int n,i,j;
1070: int *w;
1071: DP g;
1072: struct oEGT tmp0,tmp1;
1073:
1074: if ( NoRA ) {
1075: *h = in; return;
1076: }
1077: if ( Print || PrintShort ) {
1078: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1079: }
1080: r = NODE_sortbi(in,0);
1081: n = length(r);
1082: w = (int *)ALLOCA(n*sizeof(int));
1083: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1084: w[i] = (int)BDY(t);
1085: for ( i = 0; i < n; i++ ) {
1086: for ( top = 0, j = n-1; j >= 0; j-- )
1087: if ( j != i ) {
1088: MKNODE(t,(pointer)w[j],top); top = t;
1089: }
1090: get_eg(&tmp0);
1091: if ( PCoeffs )
1092: dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1093: else
1094: _dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1095: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1096: if ( Print || PrintShort ) {
1097: fprintf(asir_out,"."); fflush(asir_out);
1098: }
1099: w[i] = newps_mod(g,m);
1100: }
1101: for ( top = 0, j = n-1; j >= 0; j-- ) {
1102: MKNODE(t,(pointer)w[j],top); top = t;
1103: }
1104: *h = top;
1105: if ( Print || PrintShort )
1106: fprintf(asir_out,"\n");
1107: }
1108:
1109: int newps(a,m,subst)
1110: DP a;
1111: int m;
1112: NODE subst;
1113: {
1114: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1115: return -1;
1116: if ( psn == pslen ) {
1117: pslen *= 2;
1118: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1119: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1120: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1121: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1122: if ( m )
1123: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1124: }
1125: if ( Demand ) {
1126: if ( doing_f4 )
1127: ps[psn] = a;
1128: else
1129: ps[psn] = 0;
1130: dp_save(psn,a);
1131: } else
1132: ps[psn] = a;
1133: psh[psn] = BDY(a)->dl;
1134: pss[psn] = a->sugar;
1135: psc[psn] = BDY(a)->c;
1136: if ( m )
1137: _dp_mod(a,m,subst,&psm[psn]);
1138: return psn++;
1139: }
1140:
1141: int newps_nosave(a,m,subst)
1142: DP a;
1143: int m;
1144: NODE subst;
1145: {
1146: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1147: return -1;
1148: if ( psn == pslen ) {
1149: pslen *= 2;
1150: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1151: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1152: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1153: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1154: if ( m )
1155: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1156: }
1157: ps[psn] = 0;
1158: psh[psn] = BDY(a)->dl;
1159: pss[psn] = a->sugar;
1160: psc[psn] = BDY(a)->c;
1161: if ( m )
1162: _dp_mod(a,m,subst,&psm[psn]);
1163: return psn++;
1164: }
1165:
1166: int newps_mod(a,m)
1167: DP a;
1168: int m;
1169: {
1170: if ( psn == pslen ) {
1171: pslen *= 2;
1172: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1173: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1174: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1175: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P)); /* XXX */
1176: }
1177: ps[psn] = a;
1178: psh[psn] = BDY(ps[psn])->dl;
1179: pss[psn] = ps[psn]->sugar;
1180: return psn++;
1181: }
1182:
1183: void reducebase_dehomo(f,g)
1184: NODE f,*g;
1185: {
1186: int n,i,j,k;
1187: int *r;
1188: DL *w,d;
1189: DP u;
1190: NODE t,top;
1191:
1192: n = length(f);
1193: w = (DL *)ALLOCA(n*sizeof(DL));
1194: r = (int *)ALLOCA(n*sizeof(int));
1195: for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {
1196: r[i] = (int)BDY(t); w[i] = psh[r[i]];
1197: }
1198: for ( i = 0; i < n; i++ ) {
1199: for ( j = 0, d = w[i]; j < n; j++ ) {
1200: if ( j != i ) {
1201: for ( k = 0; k < NVars; k++ )
1202: if ( d->d[k] < w[j]->d[k] )
1203: break;
1204: if ( k == NVars )
1205: break;
1206: }
1207: }
1208: if ( j != n )
1209: r[i] = -1;
1210: }
1211: for ( top = 0, i = n-1; i >= 0; i-- )
1212: if ( r[i] >= 0 ) {
1213: dp_load(r[i],&ps[r[i]]); dp_dehomo(ps[r[i]],&u); j = newps(u,0,0);
1214: MKNODE(t,j,top); top = t;
1215: }
1216: *g = top;
1217: }
1218:
1219: NODE append_one(f,n)
1220: NODE f;
1221: int n;
1222: {
1223: NODE t;
1224:
1225: if ( Reverse || !f ) {
1226: MKNODE(t,(pointer)n,f); return t;
1227: } else {
1228: for ( t = f; NEXT(t); t = NEXT(t) );
1229: MKNODE(NEXT(t),(pointer)n,0);
1230: return f;
1231: }
1232: }
1233:
1234: DP_pairs minp( d, prest )
1235: DP_pairs d, *prest;
1236: {
1237: register DP_pairs m, ml, p, l;
1238: register DL lcm;
1239: register int s, nv = CNVars;
1240: register int (*cmpfun)() = cmpdl;
1241:
1242: if ( !(p = NEXT(m = d)) ) {
1243: *prest = p;
1244: NEXT(m) = 0;
1245: return m;
1246: }
1247: for ( lcm = m->lcm, s = m->sugar, ml = 0, l = m; p; p = NEXT(l = p) )
1248: if ( NoSugar ? (*cmpfun)( nv, lcm, p->lcm ) >= 0 :
1249: (s > p->sugar || s == p->sugar && (*cmpfun)( nv, lcm, p->lcm ) >= 0) )
1250: ml = l, lcm = (m = p)->lcm, s = p->sugar;
1251: if ( !ml ) *prest = NEXT(m);
1252: else {
1253: NEXT(ml) = NEXT(m);
1254: *prest = d;
1255: }
1256: NEXT(m) = 0;
1257: return m;
1258: }
1259:
1260: void minsugar(d,dm,dr)
1261: DP_pairs d;
1262: DP_pairs *dm,*dr;
1263: {
1264: int msugar;
1265: DP_pairs t,dm0,dr0,dmt,drt;
1266:
1267: for ( msugar = d->sugar, t = NEXT(d); t; t = NEXT(t) )
1268: if ( t->sugar < msugar )
1269: msugar = t->sugar;
1270: dm0 = 0; dr0 = 0;
1271: for ( t = d; t; t = NEXT(t) ) {
1272: if ( t->sugar == msugar ) {
1273: NEXTDPP(dm0,dmt);
1274: dmt->dp1 = t->dp1; dmt->dp2 = t->dp2;
1275: dmt->lcm = t->lcm; dmt->sugar = t->sugar;
1276: } else {
1277: NEXTDPP(dr0,drt);
1278: drt->dp1 = t->dp1; drt->dp2 = t->dp2;
1279: drt->lcm = t->lcm; drt->sugar = t->sugar;
1280: }
1281: }
1282: if ( dm0 ) NEXT(dmt) = 0;
1283: if ( dr0 ) NEXT(drt) = 0;
1284: *dm = dm0; *dr = dr0;
1285: }
1286:
1287: NODE gb(f,m,subst)
1288: NODE f;
1289: int m;
1290: NODE subst;
1291: {
1292: int i,nh,prev,mag;
1293: NODE r,g,gall;
1294: DP_pairs d,d1;
1295: DP_pairs l;
1296: DP h,nf,nfm,dp1,dp2;
1297: MP mp;
1298: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tsp0,tsp1,tspm0,tspm1,tnp0,tnp1,tmp0,tmp1;
1299: int skip_nf_flag;
1300: double t_0;
1301: static prev_sugar = -1;
1302:
1303: Max_mag = 0;
1304: prev = 1;
1305: doing_f4 = 0;
1306: if ( m ) {
1307: psm = (DP *)MALLOC(pslen*sizeof(DP));
1308: for ( i = 0; i < psn; i++ )
1309: if ( psh[i] && !validhc(psc[i],m,subst) )
1310: return 0;
1311: else
1312: _dp_mod(ps[i],m,subst,&psm[i]);
1313: }
1314: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1315: i = (int)BDY(r);
1316: d = updpairs(d,g,i);
1317: g = updbase(g,i);
1318: gall = append_one(gall,i);
1319: }
1320: while ( d ) {
1321: get_eg(&tmp0);
1322: l = minp(d,&d1); d = d1;
1323: get_eg(&tmp1); add_eg(&eg_mp,&tmp0,&tmp1);
1324: if ( m ) {
1325: get_eg(&tspm0);
1326: _dp_sp_mod(psm[l->dp1],psm[l->dp2],m,&h);
1327: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1);
1328: get_eg(&tnfm0);
1329: _dp_nf_mod(gall,h,psm,m,0,&nfm);
1330: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1331: } else
1332: nfm = (DP)1;
1333: if ( nfm ) {
1334: get_eg(&tsp0);
1335: if ( Demand ) {
1336: if ( dp_load_t(psn,&nf) ) {
1337: skip_nf_flag = 1;
1338: get_eg(&tsp1); add_eg(&eg_nf,&tsp0,&tsp1);
1339: tnf0=tsp0; tnf1=tsp1;
1340: goto skip_nf;
1341: } else {
1342: skip_nf_flag = 0;
1343: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1344: dp_sp(dp1,dp2,&h);
1345: }
1346: } else
1347: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1348: get_eg(&tsp1); add_eg(&eg_sp,&tsp0,&tsp1);
1349: get_eg(&tnf0);
1350: t_0 = get_rtime();
1351: if ( !Multiple )
1352: _dp_nf(gall,h,ps,!Top,&nf);
1353: else
1354: #if MPI
1355: _dp_nf_ptozp_mpi(gall,h,ps,!Top,Multiple,&nf);
1356: #else
1357: _dp_nf_ptozp(gall,h,ps,!Top,Multiple,&nf);
1358: #endif
1359: if ( Print )
1360: fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
1361: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1362: } else
1363: nf = 0;
1364: skip_nf:
1365: if ( nf ) {
1366: NZR++;
1367: get_eg(&tpz0);
1368: prim_part(nf,0,&h);
1369: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1370: get_eg(&tnp0);
1371: if ( Demand && skip_nf_flag )
1372: nh = newps_nosave(h,m,subst);
1373: else
1374: nh = newps(h,m,subst);
1375: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1376: if ( nh < 0 )
1377: return 0;
1378: d = updpairs(d,g,nh);
1379: g = updbase(g,nh);
1380: gall = append_one(gall,nh);
1381: if ( !dp_fcoeffs && ShowMag ) {
1382: for ( mag = 0, mp = BDY(h); mp; mp = NEXT(mp) )
1383: mag += p_mag((P)mp->c);
1384: Max_mag = MAX(Max_mag,mag);
1385: }
1386: if ( Print ) {
1387: if ( !prev )
1388: fprintf(asir_out,"\n");
1389: print_split_e(&tnf0,&tnf1); print_split_e(&tpz0,&tpz1);
1390: printdl(psh[nh]);
1391: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",
1392: l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
1393: pss[nh]);
1394: if ( ShowMag )
1395: fprintf(asir_out,",mag=%d",mag);
1396: fprintf(asir_out,"\n"); fflush(asir_out);
1397: } else if ( PrintShort ) {
1398: fprintf(asir_out,"+"); fflush(asir_out);
1399: }
1400: prev = 1;
1401: } else {
1402: if ( m )
1403: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1404: ZR++;
1405: if ( Print || PrintShort ) {
1406: if ( h && (h->sugar != prev_sugar) ) {
1407: fprintf(asir_out,"[%d]",h->sugar);
1408: prev_sugar = h->sugar;
1409: }
1410: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1411: }
1412: }
1413: }
1414: if ( Print || PrintShort )
1415: fprintf(asir_out,"gb done\n");
1416: return g;
1417: }
1418:
1419: NODE gbd(f,m,subst,dlist)
1420: NODE f;
1421: int m;
1422: NODE subst;
1423: NODE dlist;
1424: {
1425: int i,nh,prev;
1426: NODE r,g,gall;
1427: struct dp_pairs ol;
1428: DP_pairs l;
1429: NODE pair;
1430: DP h,nf,nfm,dp1,dp2;
1431: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tsp0,tsp1,tspm0,tspm1,tnp0,tnp1;
1432:
1433: prev = 1;
1434: l = &ol;
1435: if ( m ) {
1436: psm = (DP *)MALLOC(pslen*sizeof(DP));
1437: for ( i = 0; i < psn; i++ )
1438: if ( psh[i] && !validhc(psc[i],m,subst) )
1439: return 0;
1440: else
1441: _dp_mod(ps[i],m,subst,&psm[i]);
1442: }
1443: for ( gall = g = 0, r = f; r; r = NEXT(r) ) {
1444: i = (int)BDY(r);
1445: g = updbase(g,i);
1446: gall = append_one(gall,i);
1447: }
1448: while ( dlist ) {
1449: pair = BDY((LIST)BDY(dlist)); dlist = NEXT(dlist);
1450: l->dp1 = QTOS((Q)BDY(pair)); pair = NEXT(pair);
1451: l->dp2 = QTOS((Q)BDY(pair));
1452: if ( m ) {
1453: get_eg(&tspm0);
1454: _dp_sp_mod(psm[l->dp1],psm[l->dp2],m,&h);
1455: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1);
1456: get_eg(&tnfm0);
1457: _dp_nf_mod(gall,h,psm,m,0,&nfm);
1458: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1459: } else
1460: nfm = (DP)1;
1461: if ( nfm ) {
1462: get_eg(&tsp0);
1463: if ( Demand ) {
1464: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1465: dp_sp(dp1,dp2,&h);
1466: } else
1467: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1468: get_eg(&tsp1); add_eg(&eg_sp,&tsp0,&tsp1);
1469: get_eg(&tnfm0);
1470: get_eg(&tnf0);
1471: if ( !Multiple )
1472: _dp_nf(gall,h,ps,!Top,&nf);
1473: else
1474: _dp_nf_ptozp(gall,h,ps,!Top,Multiple,&nf);
1475: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1476: } else
1477: nf = 0;
1478: if ( nf ) {
1479: NZR++;
1480: get_eg(&tpz0);
1481: prim_part(nf,0,&h);
1482: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1483: get_eg(&tnp0);
1484: nh = newps(h,m,subst);
1485: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1486: if ( nh < 0 )
1487: return 0;
1488: g = updbase(g,nh);
1489: gall = append_one(gall,nh);
1490: if ( Print ) {
1491: if ( !prev )
1492: fprintf(asir_out,"\n");
1493: print_split_eg(&tnf0,&tnf1); fflush(asir_out);
1494: 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]);
1495: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1496: }
1497: prev = 1;
1498: } else {
1499: if ( m )
1500: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1501: ZR++;
1502: if ( Print ) {
1503: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1504: }
1505: }
1506: }
1507: if ( Print )
1508: fprintf(asir_out,"gb done\n");
1509: return g;
1510: }
1511:
1512: NODE gb_mod(f,m)
1513: NODE f;
1514: int m;
1515: {
1516: int i,nh,prev;
1517: NODE r,g,gall;
1518: DP_pairs d,d1;
1519: DP_pairs l;
1520: DP h,nf;
1521: struct oEGT tnfm0,tnfm1,tspm0,tspm1,tmp0,tmp1,tpz0,tpz1;
1522:
1523: prev = 1;
1524: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1525: i = (int)BDY(r);
1526: d = updpairs(d,g,i);
1527: g = updbase(g,i);
1528: gall = append_one(gall,i);
1529: }
1530: while ( d ) {
1531: get_eg(&tmp0);
1532: l = minp(d,&d1); d = d1;
1533: get_eg(&tmp1); add_eg(&eg_mp,&tmp0,&tmp1); get_eg(&tspm0);
1534: if ( PCoeffs ) {
1535: dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1536: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1); get_eg(&tnfm0);
1537: dp_nf_mod(gall,h,ps,m,!Top,&nf);
1538: } else {
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: }
1543: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1544: if ( nf ) {
1545: NZR++;
1546: get_eg(&tpz0);
1547: prim_part(nf,m,&h);
1548: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1549: nh = newps_mod(h,m);
1550: if ( nh < 0 )
1551: return 0;
1552: d = updpairs(d,g,nh);
1553: g = updbase(g,nh);
1554: if ( InterReduce )
1555: interreduce_mod(g,nh,m);
1556: gall = append_one(gall,nh);
1557: if ( Print ) {
1558: if ( !prev )
1559: fprintf(asir_out,"\n");
1560: print_split_eg(&tnfm0,&tnfm1); fflush(asir_out);
1561: 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]);
1562: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1563: }
1564: prev = 1;
1565: } else {
1566: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1567: ZR++;
1568: if ( Print ) {
1569: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1570: }
1571: }
1572: }
1573: if ( Print )
1574: fprintf(asir_out,"gb_mod done\n");
1575: return g;
1576: }
1577:
1578: void interreduce_mod(g,n,m)
1579: NODE g;
1580: int n,m;
1581: {
1582: DP nf;
1583: NODE gn,t;
1584: int i;
1585:
1586: MKNODE(gn,(pointer)n,0);
1587: if ( PCoeffs )
1588: for ( t = g; t; t = NEXT(t) ) {
1589: i = (int)BDY(t);
1590: if ( i != n ) {
1591: dp_nf_mod(gn,ps[i],ps,m,1,&nf); ps[i] = nf;
1592: }
1593: }
1594: else
1595: for ( t = g; t; t = NEXT(t) ) {
1596: i = (int)BDY(t);
1597: if ( i != n ) {
1598: _dp_nf_mod(gn,ps[i],ps,m,1,&nf); ps[i] = nf;
1599: }
1600: }
1601: }
1602:
1603: void _dp_nf_mod(b,g,ps,mod,full,rp)
1604: NODE b;
1605: DP g;
1606: DP *ps;
1607: int mod,full;
1608: DP *rp;
1609: {
1610: DP u,p,d,s,t;
1611: NODE l;
1612: MP m,mr;
1613: int sugar,psugar;
1614:
1615: if ( !g ) {
1616: *rp = 0; return;
1617: }
1618: sugar = g->sugar;
1619: for ( d = 0; g; ) {
1620: for ( u = 0, l = b; l; l = NEXT(l) ) {
1621: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1622: _dp_red_mod(g,p,mod,&u);
1623: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1624: sugar = MAX(sugar,psugar);
1625: if ( !u ) {
1626: if ( d )
1627: d->sugar = sugar;
1628: *rp = d; return;
1629: }
1630: break;
1631: }
1632: }
1633: if ( u )
1634: g = u;
1635: else if ( !full ) {
1636: if ( g ) {
1637: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1638: }
1639: *rp = g; return;
1640: } else {
1641: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1642: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1643: _addmd(CO,mod,d,t,&s); d = s;
1644: dp_rest(g,&t); g = t;
1645: }
1646: }
1647: if ( d )
1648: d->sugar = sugar;
1649: *rp = d;
1650: }
1651:
1652: void _dp_nf(b,g,ps,full,rp)
1653: NODE b;
1654: DP g;
1655: DP *ps;
1656: int full;
1657: DP *rp;
1658: {
1659: DP u,p,d,s,t;
1660: P dmy;
1661: NODE l;
1662: MP m,mr;
1663: int sugar,psugar;
1664:
1665: if ( !g ) {
1666: *rp = 0; return;
1667: }
1668: sugar = g->sugar;
1669: for ( d = 0; g; ) {
1670: for ( u = 0, l = b; l; l = NEXT(l) ) {
1671: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
1672: dp_load((int)BDY(l),&p);
1673: dp_red(d,g,p,&t,&u,&dmy);
1674: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1675: sugar = MAX(sugar,psugar);
1676: if ( !u ) {
1677: if ( d )
1678: d->sugar = sugar;
1679: *rp = d; return;
1680: }
1681: d = t;
1682: break;
1683: }
1684: }
1685: if ( u )
1686: g = u;
1687: else if ( !full ) {
1688: if ( g ) {
1689: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1690: }
1691: *rp = g; return;
1692: } else {
1693: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1694: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1695: addd(CO,d,t,&s); d = s;
1696: dp_rest(g,&t); g = t;
1697: }
1698: }
1699: if ( d )
1700: d->sugar = sugar;
1701: *rp = d;
1702: }
1703:
1704: #define SAFENM(q) ((q)?NM(q):0)
1705:
1706: double pz_t_e, pz_t_d, pz_t_d1, pz_t_c, im_t_s, im_t_r;
1707:
1708: void _dp_nf_ptozp(b,g,ps,full,multiple,r)
1709: NODE b;
1710: DP g;
1711: DP *ps;
1712: int full,multiple;
1713: DP *r;
1714: {
1715: DP u,dp,rp,t,t1,t2,red,shift;
1716: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
1717: N gn,tn,cn;
1718: NODE l;
1719: MP m,mr;
1720: int hmag,denom;
1721: int sugar,psugar;
1722: NODE dist;
1723: STRING imul;
1724: int ndist;
1725: int kara_bit;
1726: extern int kara_mag;
1727: double get_rtime();
1728: double t_0,t_00,tt,ttt,t_p,t_m,t_m1,t_m2,t_s,t_g,t_a;
1729:
1730: if ( !g ) {
1731: *r = 0; return;
1732: }
1733: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1734: t_p = t_m = t_m1 = t_m2 = t_s = t_g = t_a = 0;
1735:
1736: denom = Denominator?Denominator:1;
1737: hmag = multiple*HMAG(g)/denom;
1738: kara_bit = kara_mag*27; /* XXX */
1739: if ( Dist ) {
1740: dist = BDY(Dist);
1741: ndist = length(dist);
1742: }
1743: sugar = g->sugar;
1744:
1745: dc = 0; dp = 0; rc = ONE; rp = g;
1746: MKSTR(imul,"dp_imul_index");
1747:
1748: for ( ; rp; ) {
1749: for ( u = 0, l = b; l; l = NEXT(l) ) {
1750: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
1751: t_0 = get_rtime();
1752: dp_load((int)BDY(l),&red);
1753: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
1754: gcdn(NM(hr),NM(hred),&gn);
1755: divsn(NM(hred),gn,&tn); NTOQ(tn,SGN(hred),cr);
1756: divsn(NM(hr),gn,&tn); NTOQ(tn,SGN(hr),cred);
1757: tt = get_rtime()-t_0; t_p += tt;
1758:
1759: t_0 = get_rtime();
1760: dp_subd(rp,red,&shift);
1761: #if INET
1762: if ( Dist && ndist && HMAG(red) > kara_bit ) {
1763: NODE n0,n1,n2,n3;
1764: int i,s,id;
1765: Obj dmy;
1766: Q ind;
1767:
1768: if ( Print )
1769: fprintf(asir_out,"d");
1770: i = (int)BDY(l); STOQ(i,ind);
1771: chsgnp((P)cred,(P *)&mcred);
1772:
1773: MKNODE(n3,ind,0); MKNODE(n2,mcred,n3);
1774: MKNODE(n1,imul,n2); MKNODE(n0,BDY(dist),n1);
1775: Pox_rpc(n0,&dmy);
1776: muldc(CO,rp,(P)cr,&t);
1777: NEXT(n0)=0;
1778: Pox_pop_local(n0,&t1);
1779: } else
1780: #endif
1781: {
1782: /*
1783: if ( Print )
1784: fprintf(asir_out,"l");
1785: */
1786: t_00 = get_rtime();
1787: muldc(CO,rp,(P)cr,&t);
1788: ttt = get_rtime()-t_00; t_m1 += ttt/dp_nt(rp);
1789: t_00 = get_rtime();
1790: chsgnp((P)cred,(P *)&mcred);
1791: muldc(CO,red,(P)mcred,&t1);
1792: ttt = get_rtime()-t_00; t_m2 += ttt/dp_nt(red);
1793: }
1794: t_00 = get_rtime();
1795: muld(CO,shift,t1,&t2);
1796: addd(CO,t,t2,&u);
1797: tt = get_rtime(); t_m += tt-t_0;
1798: ttt = get_rtime(); t_s += ttt-t_00;
1799:
1800: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
1801: sugar = MAX(sugar,psugar);
1802: if ( !u ) {
1803: if ( dp )
1804: dp->sugar = sugar;
1805: *r = dp;
1806: goto final;
1807: }
1808: break;
1809: }
1810: }
1811: if ( u ) {
1812: if ( HMAG(u) > hmag ) {
1813: t_0 = get_rtime();
1814: #if INET
1815: if ( Dist && HMAG(u) > kara_bit ) {
1816: if ( Print )
1817: fprintf(asir_out,"D");
1818: dp_ptozp_d(dist,ndist,u,&t);
1819: } else {
1820: if ( Print )
1821: fprintf(asir_out,"L");
1822: dp_ptozp_d(0,0,u,&t);
1823: }
1824: #else
1825: dp_ptozp(u,&t); /* XXX */
1826: #endif
1827: tt = get_rtime()-t_0; t_g += tt;
1828: t_0 = get_rtime();
1829: divsn(NM((Q)BDY(u)->c),NM((Q)BDY(t)->c),&cn); NTOQ(cn,1,cont);
1830: if ( !dp_fcoeffs && Print ) {
1831: fprintf(asir_out,"(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1832: fflush(asir_out);
1833: }
1834: mulq(cr,dc,&dcq);
1835: mulq(cont,rc,&rcq);
1836: gcdn(SAFENM(dcq),SAFENM(rcq),&gn);
1837: divsn(SAFENM(dcq),gn,&tn); NTOQ(tn,SGN(dcq),dc);
1838: divsn(SAFENM(rcq),gn,&tn); NTOQ(tn,SGN(rcq),rc);
1839: tt = get_rtime()-t_0; t_a += tt;
1840: rp = t;
1841: hmag = multiple*HMAG(rp)/denom;
1842: } else {
1843: t_0 = get_rtime();
1844: mulq(cr,dc,&dcq); dc = dcq;
1845: tt = get_rtime()-t_0; t_a += tt;
1846: rp = u;
1847: }
1848: } else if ( !full ) {
1849: if ( rp ) {
1850: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
1851: }
1852: *r = rp;
1853: goto final;
1854: } else {
1855: t_0 = get_rtime();
1856: mulq((Q)BDY(rp)->c,rc,&c);
1857: gcdn(SAFENM(dc),SAFENM(c),&gn); NTOQ(gn,1,gcd);
1858: divsn(SAFENM(dc),gn,&tn); NTOQ(tn,SGN(dc),dcq);
1859: divsn(SAFENM(c),gn,&tn); NTOQ(tn,SGN(c),cq);
1860: muldc(CO,dp,(P)dcq,&t1);
1861: m = BDY(rp); NEWMP(mr); mr->dl = m->dl; mr->c = (P)cq;
1862: NEXT(mr) = 0; MKDP(rp->nv,mr,t); t->sugar = mr->dl->td;
1863: addd(CO,t,t1,&dp);
1864: dc = gcd;
1865:
1866: dp_rest(rp,&t); rp = t;
1867: tt = get_rtime()-t_0; t_a += tt;
1868: }
1869: }
1870: if ( dp )
1871: dp->sugar = sugar;
1872: *r = dp;
1873: final:
1874: if ( Print )
1875: fprintf(asir_out,"(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
1876: t_p,t_m,t_m1,t_m2,t_s,
1877: t_g,t_a,
1878: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
1879: }
1880:
1881: void dp_imul_d();
1882: void imulv();
1883:
1884: void _dp_nf_ptozp_mpi(b,g,ps,full,multiple,r)
1885: NODE b;
1886: DP g;
1887: DP *ps;
1888: int full,multiple;
1889: DP *r;
1890: {
1891: Obj dmy;
1892: DP u,dp,rp,t,t1,t2,red,shift;
1893: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
1894: N gn,tn,cn;
1895: NODE l,n0,n1,n2,n3;
1896: MP m,mr;
1897: int i,n;
1898: int hmag,denom;
1899: int sugar,psugar;
1900: NODE dist;
1901: STRING imul;
1902: int ndist;
1903: int kara_bit;
1904: extern int kara_mag;
1905: extern int mpi_mag;
1906: Q ind;
1907: int id;
1908: double get_rtime();
1909: double t_0,t_00,tt,ttt,t_p,t_m,t_m1,t_m2,t_s,t_g,t_a;
1910:
1911: if ( !g ) {
1912: *r = 0; return;
1913: }
1914: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1915: t_p = t_m = t_m1 = t_m2 = t_s = t_g = t_a = 0;
1916:
1917: denom = Denominator?Denominator:1;
1918: hmag = multiple*HMAG(g)/denom;
1919: kara_bit = kara_mag*27; /* XXX */
1920: if ( Dist ) {
1921: dist = BDY(Dist);
1922: ndist = length(dist);
1923: }
1924: sugar = g->sugar;
1925:
1926: dc = 0; dp = 0; rc = ONE; rp = g;
1927: MKSTR(imul,"dp_imul_index");
1928:
1929: for ( ; rp; ) {
1930: for ( u = 0, l = b; l; l = NEXT(l) ) {
1931: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
1932: t_0 = get_rtime();
1933: dp_load((int)BDY(l),&red);
1934: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
1935: gcdn(NM(hr),NM(hred),&gn);
1936: divsn(NM(hred),gn,&tn); NTOQ(tn,SGN(hred),cr);
1937: divsn(NM(hr),gn,&tn); NTOQ(tn,SGN(hr),cred);
1938: tt = get_rtime()-t_0; t_p += tt;
1939:
1940: t_0 = get_rtime();
1941: dp_subd(rp,red,&shift);
1942: t_00 = get_rtime();
1943: if ( Dist && ndist
1944: && HMAG(rp) > mpi_mag && p_mag((P)cr) > mpi_mag ) {
1945: if ( Print ) fprintf(asir_out,"~");
1946: dp_imul_d(rp,cr,&t);
1947: } else {
1948: if ( Print ) fprintf(asir_out,"_");
1949: muldc(CO,rp,(P)cr,&t);
1950: }
1951: ttt = get_rtime()-t_00; t_m1 += ttt/dp_nt(rp);
1952:
1953: t_00 = get_rtime();
1954: chsgnp((P)cred,(P *)&mcred);
1955: if ( Dist && ndist
1956: && HMAG(red) > mpi_mag && p_mag((P)mcred) > mpi_mag ) {
1957: if ( Print ) fprintf(asir_out,"=");
1958: dp_imul_d(red,mcred,&t1);
1959: } else {
1960: if ( Print ) fprintf(asir_out,"_");
1961: muldc(CO,red,(P)mcred,&t1);
1962: }
1963: ttt = get_rtime()-t_00; t_m2 += ttt/dp_nt(red);
1964:
1965: t_00 = get_rtime();
1966: muld(CO,shift,t1,&t2);
1967: addd(CO,t,t2,&u);
1968: tt = get_rtime(); t_m += tt-t_0;
1969: ttt = get_rtime(); t_s += ttt-t_00;
1970:
1971: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
1972: sugar = MAX(sugar,psugar);
1973: if ( !u ) {
1974: if ( dp )
1975: dp->sugar = sugar;
1976: *r = dp;
1977: goto final;
1978: }
1979: break;
1980: }
1981: }
1982: if ( u ) {
1983: if ( HMAG(u) > hmag ) {
1984: t_0 = get_rtime();
1985: if ( Dist && HMAG(u) > mpi_mag ) {
1986: if ( Print )
1987: fprintf(asir_out,"D");
1988: dp_ptozp_d(dist,ndist,u,&t);
1989: } else {
1990: if ( Print )
1991: fprintf(asir_out,"L");
1992: dp_ptozp_d(0,0,u,&t);
1993: }
1994: tt = get_rtime()-t_0; t_g += tt;
1995: t_0 = get_rtime();
1996: divsn(NM((Q)BDY(u)->c),NM((Q)BDY(t)->c),&cn); NTOQ(cn,1,cont);
1997: if ( Print ) {
1998: fprintf(asir_out,"(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1999: fflush(asir_out);
2000: }
2001: mulq(cr,dc,&dcq);
2002: mulq(cont,rc,&rcq);
2003: gcdn(SAFENM(dcq),SAFENM(rcq),&gn);
2004: divsn(SAFENM(dcq),gn,&tn); NTOQ(tn,SGN(dcq),dc);
2005: divsn(SAFENM(rcq),gn,&tn); NTOQ(tn,SGN(rcq),rc);
2006: tt = get_rtime()-t_0; t_a += tt;
2007: rp = t;
2008: hmag = multiple*HMAG(rp)/denom;
2009: } else {
2010: t_0 = get_rtime();
2011: mulq(cr,dc,&dcq); dc = dcq;
2012: tt = get_rtime()-t_0; t_a += tt;
2013: rp = u;
2014: }
2015: } else if ( !full ) {
2016: if ( rp ) {
2017: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
2018: }
2019: *r = rp;
2020: goto final;
2021: } else {
2022: t_0 = get_rtime();
2023: mulq((Q)BDY(rp)->c,(Q)rc,(Q *)&c);
2024: gcdn(SAFENM(dc),SAFENM(c),&gn); NTOQ(gn,1,gcd);
2025: divsn(SAFENM(dc),gn,&tn); NTOQ(tn,SGN(dc),dcq);
2026: divsn(SAFENM(c),gn,&tn); NTOQ(tn,SGN(c),cq);
2027: muldc(CO,dp,(P)dcq,&t1);
2028: m = BDY(rp); NEWMP(mr); mr->dl = m->dl; mr->c = (P)cq;
2029: NEXT(mr) = 0; MKDP(rp->nv,mr,t); t->sugar = mr->dl->td;
2030: addd(CO,t,t1,&dp);
2031: dc = gcd;
2032:
2033: dp_rest(rp,&t); rp = t;
2034: tt = get_rtime()-t_0; t_a += tt;
2035: }
2036: }
2037: if ( dp )
2038: dp->sugar = sugar;
2039: *r = dp;
2040: final:
2041: if ( Print )
2042: fprintf(asir_out,"(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
2043: t_p,t_m,t_m1,t_m2,t_s,
2044: t_g,t_a,
2045: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
2046: }
2047:
2048: void dp_imul_d(p,q,rp)
2049: DP p;
2050: Q q;
2051: DP *rp;
2052: {
2053: int nsep,ndist,i,j,k,l,n;
2054: double t0,t1,t2;
2055: Q *s;
2056: pointer *b;
2057: VECT c,cs,ri;
2058: VECT *r;
2059: MP m;
2060: NODE tn,dist,n0,n1,n2;
2061: Obj dmy;
2062: STRING imul;
2063:
2064: extern LIST Dist;
2065:
2066: if ( !p || !q ) {
2067: *rp = 0; return;
2068: }
2069: dist = BDY(Dist);
2070: for ( tn = dist, ndist = 0; tn; tn = NEXT(tn), ndist++ );
2071: nsep = ndist + 1;
2072: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
2073: if ( n <= nsep ) {
2074: muldc(CO,p,(P)q,rp); return;
2075: }
2076: MKSTR(imul,"imulv");
2077: t0 = get_rtime();
2078: dp_dtov(p,&c);
2079: sepvect(c,nsep,&cs);
2080: r = (VECT *)CALLOC(nsep,sizeof(VECT *));
2081: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
2082: n0 = mknode(4,BDY(tn),imul,b[i],q);
2083: Pox_rpc(n0,&dmy);
2084: }
2085: t1 = get_rtime();
2086: im_t_s += t1 - t0;
2087: imulv(b[i],q,&r[i]);
2088: t1 = get_rtime();
2089: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
2090: MKNODE(n0,BDY(tn),0);
2091: Pox_pop_local(n0,&r[i]);
2092: if ( OID(r[i]) == O_ERR ) {
2093: printexpr(CO,(Obj)r[i]);
2094: error("dp_imul_d : aborted");
2095: }
2096: }
2097: t2 = get_rtime();
2098: im_t_r += t2 - t1;
2099: s = (Q *)CALLOC(n,sizeof(Q));
2100: for ( i = j = 0; i < nsep; i++ ) {
2101: for ( k = 0, ri = r[i], l = ri->len; k < l; k++, j++ ) {
2102: s[j] = (Q)BDY(ri)[k];
2103: }
2104: }
2105: dp_vtod(s,p,rp);
2106: }
2107:
2108: void imulv(w,c,rp)
2109: VECT w;
2110: Q c;
2111: VECT *rp;
2112: {
2113: int n,i;
2114: VECT r;
2115:
2116: n = w->len;
2117: MKVECT(r,n); *rp = r;
2118: for ( i = 0; i < n; i++ )
2119: mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
2120: }
2121:
2122: DP_pairs updpairs( d, g, t)
2123: DP_pairs d;
2124: NODE /* of index */ g;
2125: int t;
2126: {
2127: register DP_pairs d1, dd, nd;
2128: struct oEGT tup0,tup1;
2129: int dl,dl1;
2130:
2131: if ( !g ) return d;
2132: get_eg(&tup0);
2133: if ( !NoCriB && d ) {
2134: dl = DPPlength(d);
2135: d = criterion_B( d, t );
2136: dl -= DPPlength(d); NBP += dl;
2137: }
2138: d1 = newpairs( g, t );
2139: if ( NEXT(d1) ) {
2140: dl = DPPlength(d1); TP += dl;
2141: d1 = criterion_M( d1 );
2142: dl1 = DPPlength(d1); NMP += (dl-dl1); dl = dl1;
2143: d1 = criterion_F( d1 );
2144: dl1 = DPPlength(d1); NFP += (dl-dl1); dl = dl1;
2145: } else
2146: dl = 1;
2147: for ( dd = 0; d1; d1 = nd ) {
2148: nd = NEXT(d1);
2149: if ( !criterion_2( d1->dp1, d1->dp2 ) ) {
2150: NEXT(d1) = dd;
2151: dd = d1;
2152: }
2153: }
2154: dl1 = DPPlength(dd); NDP += (dl-dl1);
2155: get_eg(&tup1);
2156: add_eg(&eg_up,&tup0,&tup1);
2157: if ( !(nd = d) ) return dd;
2158: while ( nd = NEXT(d1 = nd) ) ;
2159: NEXT(d1) = dd;
2160: return d;
2161: }
2162:
2163: DP_pairs newpairs( g, t )
2164: NODE /* of index */ g;
2165: register int t;
2166: {
2167: register NODE r;
2168: register DL tdl = psh[t];
2169: register int ts;
2170: register DP_pairs p, last;
2171: int dp;
2172: register DL dl;
2173: register int s;
2174:
2175: ts = pss[t] - tdl->td;
2176: for ( last = 0, r = g; r; r = NEXT(r) ) {
2177: NEXT(p = NEWDP_pairs) = last;
2178: last = p;
2179: dp = p->dp1 = (int)BDY(r); p->dp2 = t;
2180: p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
2181: p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
2182: }
2183: return last;
2184: }
2185:
2186: DP_pairs criterion_B( d, s )
2187: register DP_pairs d;
2188: int s;
2189: {
2190: register DP_pairs dd, p;
2191: register DL tij, t = psh[s], dltmp;
2192:
2193: if ( !d ) return 0;
2194: NEWDL( dltmp, CNVars );
2195: for ( dd = 0; d; d = p ) {
2196: p = NEXT(d),
2197: tij = d->lcm;
2198: if ( tij->td != lcm_of_DL(CNVars, tij, t, dltmp )->td
2199: || !dl_equal(CNVars, tij, dltmp )
2200: || (tij->td == lcm_of_DL(CNVars, psh[d->dp1], t, dltmp )->td
2201: && dl_equal(CNVars, dltmp, tij ))
2202: || (tij->td == lcm_of_DL(CNVars, psh[d->dp2], t, dltmp )->td
2203: && dl_equal(CNVars, dltmp, tij )) ) {
2204: NEXT(d) = dd;
2205: dd = d;
2206: }
2207: }
2208: return dd;
2209: }
2210:
2211: DP_pairs criterion_M( d1 )
2212: DP_pairs d1;
2213: {
2214: register DP_pairs dd, e, d3, d2, p;
2215: register DL itdl, jtdl;
2216: register int itdltd, jtdltd;
2217:
2218: for ( dd = 0, e = d1; e; e = d3 ) {
2219: if ( !(d2 = NEXT(e)) ) {
2220: NEXT(e) = dd;
2221: return e;
2222: }
2223: itdltd = (itdl = e->lcm)->td;
2224: for ( d3 = 0; d2; d2 = p ) {
2225: p = NEXT(d2),
2226: jtdltd = (jtdl = d2->lcm)->td;
2227: if ( jtdltd == itdltd )
2228: if ( dl_equal(CNVars, itdl, jtdl ) ) ;
2229: else if ( dl_redble( jtdl, itdl ) ) continue;
2230: else if ( dl_redble( itdl, jtdl ) ) goto delit;
2231: else ;
2232: else if ( jtdltd > itdltd )
2233: if ( dl_redble( jtdl, itdl ) ) continue;
2234: else ;
2235: else if ( dl_redble( itdl, jtdl ) ) goto delit;
2236: NEXT(d2) = d3;
2237: d3 = d2;
2238: }
2239: NEXT(e) = dd;
2240: dd = e;
2241: continue;
2242: /**/
2243: delit: NEXT(d2) = d3;
2244: d3 = d2;
2245: for ( ; p; p = d2 ) {
2246: d2 = NEXT(p);
2247: NEXT(p) = d3;
2248: d3 = p;
2249: }
2250: }
2251: return dd;
2252: }
2253:
2254: static DP_pairs collect_pairs_of_hdlcm( d1, prest )
2255: DP_pairs d1, *prest;
2256: {
2257: register DP_pairs w, p, r, s;
2258: register DL ti;
2259: register int td;
2260:
2261: td = (ti = (w = d1)->lcm)->td;
2262: s = NEXT(w);
2263: NEXT(w) = 0;
2264: for ( r = 0; s; s = p ) {
2265: p = NEXT(s);
2266: if ( td == s->lcm->td && dl_equal(CNVars, ti, s->lcm ) )
2267: {
2268: NEXT(s) = w;
2269: w = s;
2270: } else {
2271: NEXT(s) = r;
2272: r = s;
2273: }
2274: }
2275: *prest = r;
2276: return w;
2277: }
2278:
2279: int criterion_2( dp1, dp2 )
2280: int dp1, dp2;
2281: {
2282: register int i, *d1, *d2;
2283:
2284: d1 = psh[dp1]->d, d2 = psh[dp2]->d;
2285: for ( i = CNVars; --i >= 0; d1++, d2++ )
2286: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
2287: return 1;
2288: }
2289:
2290: DP_pairs criterion_F( d1 )
2291: DP_pairs d1;
2292: {
2293: DP_pairs rest, head;
2294: register DP_pairs last, p, r, w;
2295: register int s;
2296:
2297: for ( head = last = 0, p = d1; NEXT(p); ) {
2298: s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
2299: while ( w = NEXT(w) )
2300: if ( criterion_2( w->dp1, w->dp2 ) ) {
2301: r = w;
2302: break;
2303: } else if ( w->sugar < s ) s = (r = w)->sugar;
2304: if ( last ) NEXT(last) = r;
2305: else head = r;
2306: NEXT(last = r) = 0;
2307: if ( !(p = rest) ) return head;
2308: }
2309: if ( !last ) return p;
2310: NEXT(last) = p;
2311: return head;
2312: }
2313:
2314: NODE updbase(g,t)
2315: NODE g;
2316: int t;
2317: {
2318: g = remove_reducibles(g,t);
2319: g = append_one(g,t);
2320: return g;
2321: }
2322:
2323: NODE /* of index */ remove_reducibles( nd, newdp )
2324: NODE /* of index */ nd;
2325: int newdp;
2326: {
2327: register DL dl, dln;
2328: register NODE last, p, head;
2329: register int td;
2330:
2331: dl = psh[newdp];
2332: td = dl->td;
2333: for ( head = last = 0, p = nd; p; ) {
2334: dln = psh[(int)BDY(p)];
2335: if ( dln->td >= td && dl_redble( dln, dl ) ) {
2336: p = NEXT(p);
2337: if ( last ) NEXT(last) = p;
2338: } else {
2339: if ( !last ) head = p;
2340: p = NEXT(last = p);
2341: }
2342: }
2343: return head;
2344: }
2345:
2346: int dl_redble( dl1, dl2 )
2347: DL dl1, dl2;
2348: {
2349: register int n, *d1, *d2;
2350:
2351: for ( d1 = dl1->d, d2 = dl2->d, n = CNVars; --n >= 0; d1++, d2++ )
2352: if ( *d1 < *d2 ) return 0;
2353: return 1;
2354: }
2355:
2356: int gbcheck(f)
2357: NODE f;
2358: {
2359: int i;
2360: NODE r,g,gall;
2361: DP_pairs d,l;
2362: DP h,nf,dp1,dp2;
2363: struct oEGT tmp0,tmp1;
2364:
2365: if ( NoGC )
2366: return 1;
2367: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
2368: i = (int)BDY(r);
2369: d = updpairs(d,g,i);
2370: g = updbase(g,i);
2371: gall = append_one(gall,i);
2372: }
2373: if ( Print || PrintShort ) {
2374: fprintf(asir_out,"gbcheck total %d pairs\n",DPPlength(d)); fflush(asir_out);
2375: }
2376: while ( d ) {
2377: l = d; d = NEXT(d);
2378: get_eg(&tmp0);
2379: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2); dp_sp(dp1,dp2,&h);
2380: _dp_nf(gall,h,ps,1,&nf);
2381: get_eg(&tmp1); add_eg(&eg_gc,&tmp0,&tmp1);
2382: if ( Print || PrintShort ) {
2383: fprintf(asir_out,"."); fflush(asir_out);
2384: }
2385: if ( nf )
2386: return 0;
2387: }
2388: if ( Print || PrintShort )
2389: fprintf(asir_out,"\n");
2390: return 1;
2391: }
2392:
2393: int membercheck(f,x)
2394: NODE f,x;
2395: {
2396: DP g;
2397: struct oEGT tmp0,tmp1;
2398:
2399: if ( NoMC )
2400: return 1;
2401: if ( Print || PrintShort ) {
2402: fprintf(asir_out,"membercheck\n"); fflush(asir_out);
2403: }
2404: for ( ; f; f = NEXT(f) ) {
2405: get_eg(&tmp0);
2406: _dp_nf(x,(DP)BDY(f),ps,1,&g);
2407: get_eg(&tmp1); add_eg(&eg_mc,&tmp0,&tmp1);
2408: if ( Print ) {
2409: print_split_eg(&tmp0,&tmp1); fflush(asir_out);
2410: } else if ( PrintShort ) {
2411: fprintf(asir_out,"."); fflush(asir_out);
2412: }
2413: if ( g )
2414: return 0;
2415: }
2416: if ( Print || PrintShort )
2417: fprintf(asir_out,"\n");
2418: return 1;
2419: }
2420:
2421: void dp_set_flag(name,value)
2422: Obj name,value;
2423: {
2424: char *n;
2425: int v;
2426:
2427: if ( OID(name) != O_STR )
2428: return;
2429: n = BDY((STRING)name);
2430: if ( !strcmp(n,"Demand") ) {
2431: Demand = value ? BDY((STRING)value) : 0; return;
2432: }
2433: if ( !strcmp(n,"Dist") ) {
2434: Dist = (LIST)value; return;
2435: }
2436: if ( value && OID(value) != O_N )
2437: return;
2438: v = QTOS((Q)value);
2439: if ( !strcmp(n,"NoSugar") )
2440: NoSugar = v;
2441: else if ( !strcmp(n,"NoCriB") )
2442: NoCriB = v;
2443: else if ( !strcmp(n,"NoGC") )
2444: NoGC = v;
2445: else if ( !strcmp(n,"NoMC") )
2446: NoMC = v;
2447: else if ( !strcmp(n,"NoRA") )
2448: NoRA = v;
2449: else if ( !strcmp(n,"NoGCD") )
2450: NoGCD = v;
2451: else if ( !strcmp(n,"Top") )
2452: Top = v;
2453: else if ( !strcmp(n,"ShowMag") )
2454: ShowMag = v;
2455: else if ( !strcmp(n,"PrintShort") )
2456: PrintShort = v;
2457: else if ( !strcmp(n,"Print") )
2458: Print = v;
2459: else if ( !strcmp(n,"Stat") )
2460: Stat = v;
2461: else if ( !strcmp(n,"Reverse") )
2462: Reverse = v;
2463: else if ( !strcmp(n,"InterReduce") )
2464: InterReduce = v;
2465: else if ( !strcmp(n,"Multiple") )
2466: Multiple = v;
2467: else if ( !strcmp(n,"Denominator") )
2468: Denominator = v;
2469: else if ( !strcmp(n,"PtozpRA") )
2470: PtozpRA = v;
2471: }
2472:
2473: void dp_make_flaglist(list)
2474: LIST *list;
2475: {
2476: Q v;
2477: STRING name,path;
2478: NODE n,n1;
2479:
2480: STOQ(Multiple,v); MKNODE(n,v,0); MKSTR(name,"Multiple"); MKNODE(n1,name,n); n = n1;
2481: STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
2482: MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
2483: STOQ(InterReduce,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"InterReduce"); MKNODE(n1,name,n); n = n1;
2484: STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
2485: STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
2486: STOQ(Print,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Print"); MKNODE(n1,name,n); n = n1;
2487: STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
2488: STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
2489: STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
2490: STOQ(NoGCD,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGCD"); MKNODE(n1,name,n); n = n1;
2491: STOQ(NoRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoRA"); MKNODE(n1,name,n); n = n1;
2492: STOQ(NoMC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoMC"); MKNODE(n1,name,n); n = n1;
2493: STOQ(NoGC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGC"); MKNODE(n1,name,n); n = n1;
2494: STOQ(NoCriB,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoCriB"); MKNODE(n1,name,n); n = n1;
2495: STOQ(NoSugar,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoSugar"); MKNODE(n1,name,n); n = n1;
2496: if ( Demand )
2497: MKSTR(path,Demand);
2498: else
2499: path = 0;
2500: MKNODE(n1,path,n); n = n1; MKSTR(name,"Demand"); MKNODE(n1,name,n); n = n1;
2501: MKLIST(*list,n);
2502: }
2503:
2504: #define DELIM '/'
2505:
2506: #if defined(THINK_C)
2507: #undef DELIM
2508: #define DELIM ':'
2509: #endif
2510:
2511: void dp_save(index,p)
2512: int index;
2513: DP p;
2514: {
2515: FILE *fp;
2516: char path[BUFSIZ];
2517:
2518: sprintf(path,"%s%c%d",Demand,DELIM,index);
2519: #if defined(VISUAL) || defined(THINK_C)
2520: if ( !(fp = fopen(path,"wb") ) )
2521: #else
2522: if ( !(fp = fopen(path,"w") ) )
2523: #endif
2524: error("dp_save : cannot open a file");
2525: savevl(fp,VC); saveobj(fp,(Obj)p); fclose(fp);
2526: }
2527:
2528: void dp_load(index,p)
2529: int index;
2530: DP *p;
2531: {
2532: FILE *fp;
2533: char path[BUFSIZ];
2534:
2535: if ( !Demand || ps[index] )
2536: *p = ps[index];
2537: else {
2538: sprintf(path,"%s%c%d",Demand,DELIM,index);
2539: #if defined(VISUAL) || defined(THINK_C)
2540: if ( !(fp = fopen(path,"rb") ) )
2541: #else
2542: if ( !(fp = fopen(path,"r") ) )
2543: #endif
2544: error("dp_load : cannot open a file");
2545: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp);
2546: }
2547: }
2548:
2549: int dp_load_t(index,p)
2550: int index;
2551: DP *p;
2552: {
2553: FILE *fp;
2554: char path[BUFSIZ];
2555:
2556: sprintf(path,"%s%c%d",Demand,DELIM,index);
2557: #if defined(VISUAL) || defined(THINK_C)
2558: if ( !(fp = fopen(path,"rb") ) )
2559: #else
2560: if ( !(fp = fopen(path,"r") ) )
2561: #endif
2562: return 0;
2563: else {
2564: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp); return 1;
2565: }
2566: }
2567:
2568: void init_stat() {
2569: init_eg(&eg_up); init_eg(&eg_sp); init_eg(&eg_spm);
2570: init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
2571: init_eg(&eg_pz); init_eg(&eg_np); init_eg(&eg_mp);
2572: init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
2573: ZR = NZR = TP = NBP = NMP = NFP = NDP = 0;
2574: }
2575:
2576: void print_stat() {
2577: if ( !Print && !Stat )
2578: return;
2579: print_eg("UP",&eg_up); print_eg("SP",&eg_sp); print_eg("SPM",&eg_spm);
2580: print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
2581: print_eg("PZ",&eg_pz); print_eg("NP",&eg_np); print_eg("MP",&eg_mp);
2582: print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
2583: 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);
2584: }
2585:
2586: /*
2587: * Old codes
2588: */
2589:
2590: #if 0
2591: void _dp_nf_ptozp(b,g,ps,full,multiple,rp)
2592: NODE b;
2593: DP g;
2594: DP *ps;
2595: int full,multiple;
2596: DP *rp;
2597: {
2598: DP u,p,d,s,t,dmy;
2599: NODE l;
2600: MP m,mr;
2601: int i,n;
2602: int hmag,denom;
2603: int sugar,psugar;
2604: NODE dist;
2605: int ndist;
2606: int kara_bit;
2607: extern int kara_mag;
2608:
2609: if ( !g ) {
2610: *rp = 0; return;
2611: }
2612: denom = Denominator?Denominator:1;
2613: hmag = multiple*HMAG(g)/denom;
2614: kara_bit = kara_mag*27; /* XXX */
2615: if ( Dist ) {
2616: dist = BDY(Dist);
2617: ndist = length(dist);
2618: }
2619: sugar = g->sugar;
2620: for ( d = 0; g; ) {
2621: for ( u = 0, l = b; l; l = NEXT(l) ) {
2622: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
2623: dp_load((int)BDY(l),&p);
2624: dp_red(d,g,p,&t,&u,&dmy);
2625: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
2626: sugar = MAX(sugar,psugar);
2627: if ( !u ) {
2628: if ( d )
2629: d->sugar = sugar;
2630: *rp = d; return;
2631: }
2632: d = t;
2633: break;
2634: }
2635: }
2636: if ( u ) {
2637: g = u;
2638: if ( d ) {
2639: if ( HMAG(d) > hmag ) {
2640: if ( Dist && HMAG(g) > kara_bit )
2641: dp_ptozp2_d(dist,ndist,d,g,&t,&u);
2642: else
2643: dp_ptozp2(d,g,&t,&u);
2644: d = t; g = u;
2645: hmag = multiple*HMAG(d)/denom;
2646: }
2647: } else {
2648: if ( HMAG(g) > hmag ) {
2649: #if INET
2650: if ( Dist && HMAG(g) > kara_bit ) {
2651: dp_ptozp_d(dist,ndist,g,&t);
2652: } else
2653: #endif
2654: dp_ptozp(g,&t);
2655: g = t;
2656: hmag = multiple*HMAG(g)/denom;
2657: }
2658: }
2659: }
2660: else if ( !full ) {
2661: if ( g ) {
2662: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
2663: }
2664: *rp = g; return;
2665: } else {
2666: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
2667: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
2668: addd(CO,d,t,&s); d = s;
2669: dp_rest(g,&t); g = t;
2670: }
2671: }
2672: if ( d )
2673: d->sugar = sugar;
2674: *rp = d;
2675: }
2676:
2677: int criterion_2(dp1,dp2)
2678: {
2679: DP g1,g2,g,s1,s2;
2680:
2681: monomial_gcd(ps[dp1],&g1); monomial_gcd(ps[dp2],&g2);
2682: dp_gcd(g1,g2,&g);
2683: dp_subd(ps[dp1],g,&s1); dp_subd(ps[dp2],g,&s2);
2684: return _criterion_2(s1,s2);
2685: }
2686:
2687: int _criterion_2( dp1, dp2 )
2688: DP dp1, dp2;
2689: {
2690: register int i, *d1, *d2;
2691:
2692: d1 = BDY(dp1)->dl->d; d2 = BDY(dp2)->dl->d;
2693: for ( i = CNVars; --i >= 0; d1++, d2++ )
2694: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
2695: return 1;
2696: }
2697:
2698: void dp_gcd(p1,p2,rp)
2699: DP p1,p2;
2700: DP *rp;
2701: {
2702: int i,n,td;
2703: DL d1,d2,d;
2704: MP m;
2705: DP s;
2706:
2707: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
2708: NEWDL(d,n);
2709: for ( i = 0, td = 0; i < n; i++ ) {
2710: d->d[i] = MIN(d1->d[i],d2->d[i]);
2711: td += d->d[i];
2712: }
2713: d->td = td;
2714: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
2715: *rp = s;
2716: }
2717:
2718: void monomial_gcd(p,rp)
2719: DP p;
2720: DP *rp;
2721: {
2722: int n,i,td;
2723: DL d,d1;
2724: MP m;
2725: DP s;
2726:
2727: n = p->nv; m = BDY(p); d = m->dl;
2728: NEWDL(d1,n);
2729: for ( i = 0; i < n; i++ )
2730: d1->d[i] = d->d[i];
2731: for ( m = NEXT(m); m; m = NEXT(m) ) {
2732: d = m->dl;
2733: for ( i = 0; i < n; i++ )
2734: d1->d[i] = MIN(d1->d[i],d->d[i]);
2735: }
2736: for ( i = 0, td = 0; i < n; i++ )
2737: td += d1->d[i];
2738: NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
2739: *rp = s;
2740: }
2741: #endif