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