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