Annotation of OpenXM_contrib2/asir2000/builtin/gr.c, Revision 1.6
1.6 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.5 2000/05/29 08:54:45 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.6 ! noro 717: _mulmd(CO,m,sdm,ps[(int)BDY(r)],&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;
1.6 ! noro 1358: int new_sugar;
1.1 noro 1359: static prev_sugar = -1;
1360:
1361: Max_mag = 0;
1362: prev = 1;
1363: doing_f4 = 0;
1364: if ( m ) {
1365: psm = (DP *)MALLOC(pslen*sizeof(DP));
1366: for ( i = 0; i < psn; i++ )
1367: if ( psh[i] && !validhc(psc[i],m,subst) )
1368: return 0;
1369: else
1370: _dp_mod(ps[i],m,subst,&psm[i]);
1371: }
1372: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1373: i = (int)BDY(r);
1374: d = updpairs(d,g,i);
1375: g = updbase(g,i);
1376: gall = append_one(gall,i);
1377: }
1378: while ( d ) {
1379: get_eg(&tmp0);
1380: l = minp(d,&d1); d = d1;
1381: get_eg(&tmp1); add_eg(&eg_mp,&tmp0,&tmp1);
1382: if ( m ) {
1383: get_eg(&tspm0);
1.6 ! noro 1384: _dp_sp_mod_dup(psm[l->dp1],psm[l->dp2],m,&h);
! 1385: new_sugar = h->sugar;
1.1 noro 1386: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1);
1387: get_eg(&tnfm0);
1.6 ! noro 1388: _dp_nf_mod_destructive(gall,h,psm,m,0,&nfm);
1.1 noro 1389: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1390: } else
1391: nfm = (DP)1;
1392: if ( nfm ) {
1393: get_eg(&tsp0);
1394: if ( Demand ) {
1395: if ( dp_load_t(psn,&nf) ) {
1396: skip_nf_flag = 1;
1397: get_eg(&tsp1); add_eg(&eg_nf,&tsp0,&tsp1);
1398: tnf0=tsp0; tnf1=tsp1;
1399: goto skip_nf;
1400: } else {
1401: skip_nf_flag = 0;
1402: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1403: dp_sp(dp1,dp2,&h);
1404: }
1405: } else
1406: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1.6 ! noro 1407: new_sugar = h->sugar;
1.1 noro 1408: get_eg(&tsp1); add_eg(&eg_sp,&tsp0,&tsp1);
1409: get_eg(&tnf0);
1410: t_0 = get_rtime();
1411: if ( !Multiple )
1412: _dp_nf(gall,h,ps,!Top,&nf);
1413: else
1414: #if MPI
1415: _dp_nf_ptozp_mpi(gall,h,ps,!Top,Multiple,&nf);
1416: #else
1417: _dp_nf_ptozp(gall,h,ps,!Top,Multiple,&nf);
1418: #endif
1419: if ( Print )
1420: fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
1421: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1422: } else
1423: nf = 0;
1424: skip_nf:
1425: if ( nf ) {
1426: NZR++;
1427: get_eg(&tpz0);
1428: prim_part(nf,0,&h);
1429: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1430: get_eg(&tnp0);
1431: if ( Demand && skip_nf_flag )
1432: nh = newps_nosave(h,m,subst);
1433: else
1434: nh = newps(h,m,subst);
1435: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1436: if ( nh < 0 )
1437: return 0;
1438: d = updpairs(d,g,nh);
1439: g = updbase(g,nh);
1440: gall = append_one(gall,nh);
1441: if ( !dp_fcoeffs && ShowMag ) {
1442: for ( mag = 0, mp = BDY(h); mp; mp = NEXT(mp) )
1443: mag += p_mag((P)mp->c);
1444: Max_mag = MAX(Max_mag,mag);
1445: }
1446: if ( Print ) {
1447: if ( !prev )
1448: fprintf(asir_out,"\n");
1449: print_split_e(&tnf0,&tnf1); print_split_e(&tpz0,&tpz1);
1450: printdl(psh[nh]);
1451: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",
1452: l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
1453: pss[nh]);
1454: if ( ShowMag )
1455: fprintf(asir_out,",mag=%d",mag);
1456: fprintf(asir_out,"\n"); fflush(asir_out);
1457: } else if ( PrintShort ) {
1458: fprintf(asir_out,"+"); fflush(asir_out);
1459: }
1460: prev = 1;
1461: } else {
1462: if ( m )
1463: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1464: ZR++;
1465: if ( Print || PrintShort ) {
1.6 ! noro 1466: if ( new_sugar != prev_sugar ) {
! 1467: fprintf(asir_out,"[%d]",new_sugar);
! 1468: prev_sugar = new_sugar;
1.1 noro 1469: }
1470: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1471: }
1472: }
1473: }
1474: if ( Print || PrintShort )
1475: fprintf(asir_out,"gb done\n");
1476: return g;
1477: }
1478:
1479: NODE gbd(f,m,subst,dlist)
1480: NODE f;
1481: int m;
1482: NODE subst;
1483: NODE dlist;
1484: {
1485: int i,nh,prev;
1486: NODE r,g,gall;
1487: struct dp_pairs ol;
1488: DP_pairs l;
1489: NODE pair;
1490: DP h,nf,nfm,dp1,dp2;
1491: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tsp0,tsp1,tspm0,tspm1,tnp0,tnp1;
1492:
1493: prev = 1;
1494: l = &ol;
1495: if ( m ) {
1496: psm = (DP *)MALLOC(pslen*sizeof(DP));
1497: for ( i = 0; i < psn; i++ )
1498: if ( psh[i] && !validhc(psc[i],m,subst) )
1499: return 0;
1500: else
1501: _dp_mod(ps[i],m,subst,&psm[i]);
1502: }
1503: for ( gall = g = 0, r = f; r; r = NEXT(r) ) {
1504: i = (int)BDY(r);
1505: g = updbase(g,i);
1506: gall = append_one(gall,i);
1507: }
1508: while ( dlist ) {
1509: pair = BDY((LIST)BDY(dlist)); dlist = NEXT(dlist);
1510: l->dp1 = QTOS((Q)BDY(pair)); pair = NEXT(pair);
1511: l->dp2 = QTOS((Q)BDY(pair));
1512: if ( m ) {
1513: get_eg(&tspm0);
1.6 ! noro 1514: _dp_sp_mod_dup(ps[l->dp1],ps[l->dp2],m,&h);
1.1 noro 1515: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1);
1516: get_eg(&tnfm0);
1.6 ! noro 1517: _dp_nf_mod_destructive(gall,h,ps,m,!Top,&nf);
1.1 noro 1518: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1519: } else
1520: nfm = (DP)1;
1521: if ( nfm ) {
1522: get_eg(&tsp0);
1523: if ( Demand ) {
1524: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1525: dp_sp(dp1,dp2,&h);
1526: } else
1527: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1528: get_eg(&tsp1); add_eg(&eg_sp,&tsp0,&tsp1);
1529: get_eg(&tnfm0);
1530: get_eg(&tnf0);
1531: if ( !Multiple )
1532: _dp_nf(gall,h,ps,!Top,&nf);
1533: else
1534: _dp_nf_ptozp(gall,h,ps,!Top,Multiple,&nf);
1535: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1536: } else
1537: nf = 0;
1538: if ( nf ) {
1539: NZR++;
1540: get_eg(&tpz0);
1541: prim_part(nf,0,&h);
1542: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1543: get_eg(&tnp0);
1544: nh = newps(h,m,subst);
1545: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1546: if ( nh < 0 )
1547: return 0;
1548: g = updbase(g,nh);
1549: gall = append_one(gall,nh);
1550: if ( Print ) {
1551: if ( !prev )
1552: fprintf(asir_out,"\n");
1553: print_split_eg(&tnf0,&tnf1); fflush(asir_out);
1554: 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]);
1555: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1556: }
1557: prev = 1;
1558: } else {
1559: if ( m )
1560: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1561: ZR++;
1562: if ( Print ) {
1563: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1564: }
1565: }
1566: }
1567: if ( Print )
1568: fprintf(asir_out,"gb done\n");
1569: return g;
1570: }
1571:
1572: NODE gb_mod(f,m)
1573: NODE f;
1574: int m;
1575: {
1576: int i,nh,prev;
1577: NODE r,g,gall;
1578: DP_pairs d,d1;
1579: DP_pairs l;
1580: DP h,nf;
1581: struct oEGT tnfm0,tnfm1,tspm0,tspm1,tmp0,tmp1,tpz0,tpz1;
1582:
1583: prev = 1;
1584: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1585: i = (int)BDY(r);
1586: d = updpairs(d,g,i);
1587: g = updbase(g,i);
1588: gall = append_one(gall,i);
1589: }
1590: while ( d ) {
1591: get_eg(&tmp0);
1592: l = minp(d,&d1); d = d1;
1593: get_eg(&tmp1); add_eg(&eg_mp,&tmp0,&tmp1); get_eg(&tspm0);
1594: if ( PCoeffs ) {
1595: dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1596: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1); get_eg(&tnfm0);
1597: dp_nf_mod(gall,h,ps,m,!Top,&nf);
1598: } else {
1.6 ! noro 1599: _dp_sp_mod_dup(ps[l->dp1],ps[l->dp2],m,&h);
1.1 noro 1600: get_eg(&tspm1); add_eg(&eg_spm,&tspm0,&tspm1); get_eg(&tnfm0);
1.6 ! noro 1601: _dp_nf_mod_destructive(gall,h,ps,m,!Top,&nf);
1.1 noro 1602: }
1603: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1604: if ( nf ) {
1605: NZR++;
1606: get_eg(&tpz0);
1607: prim_part(nf,m,&h);
1608: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1609: nh = newps_mod(h,m);
1610: if ( nh < 0 )
1611: return 0;
1612: d = updpairs(d,g,nh);
1613: g = updbase(g,nh);
1614: if ( InterReduce )
1615: interreduce_mod(g,nh,m);
1616: gall = append_one(gall,nh);
1617: if ( Print ) {
1618: if ( !prev )
1619: fprintf(asir_out,"\n");
1620: print_split_eg(&tnfm0,&tnfm1); fflush(asir_out);
1621: 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]);
1622: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1.6 ! noro 1623: } else if ( PrintShort ) {
! 1624: fprintf(asir_out,"+"); fflush(asir_out);
1.1 noro 1625: }
1626: prev = 1;
1627: } else {
1628: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1629: ZR++;
1.6 ! noro 1630: if ( Print || PrintShort ) {
1.1 noro 1631: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1632: }
1633: }
1634: }
1.6 ! noro 1635: if ( Print || PrintShort )
1.1 noro 1636: fprintf(asir_out,"gb_mod done\n");
1637: return g;
1638: }
1639:
1640: void interreduce_mod(g,n,m)
1641: NODE g;
1642: int n,m;
1643: {
1644: DP nf;
1645: NODE gn,t;
1646: int i;
1647:
1648: MKNODE(gn,(pointer)n,0);
1649: if ( PCoeffs )
1650: for ( t = g; t; t = NEXT(t) ) {
1651: i = (int)BDY(t);
1652: if ( i != n ) {
1653: dp_nf_mod(gn,ps[i],ps,m,1,&nf); ps[i] = nf;
1654: }
1655: }
1656: else
1657: for ( t = g; t; t = NEXT(t) ) {
1658: i = (int)BDY(t);
1659: if ( i != n ) {
1660: _dp_nf_mod(gn,ps[i],ps,m,1,&nf); ps[i] = nf;
1661: }
1662: }
1663: }
1664:
1665: void _dp_nf_mod(b,g,ps,mod,full,rp)
1666: NODE b;
1667: DP g;
1668: DP *ps;
1669: int mod,full;
1670: DP *rp;
1671: {
1672: DP u,p,d,s,t;
1673: NODE l;
1674: MP m,mr;
1675: int sugar,psugar;
1676:
1677: if ( !g ) {
1678: *rp = 0; return;
1679: }
1680: sugar = g->sugar;
1681: for ( d = 0; g; ) {
1682: for ( u = 0, l = b; l; l = NEXT(l) ) {
1683: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1684: _dp_red_mod(g,p,mod,&u);
1685: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1686: sugar = MAX(sugar,psugar);
1687: if ( !u ) {
1688: if ( d )
1689: d->sugar = sugar;
1690: *rp = d; return;
1691: }
1692: break;
1693: }
1694: }
1695: if ( u )
1696: g = u;
1697: else if ( !full ) {
1698: if ( g ) {
1699: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1700: }
1701: *rp = g; return;
1702: } else {
1703: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1704: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1705: _addmd(CO,mod,d,t,&s); d = s;
1706: dp_rest(g,&t); g = t;
1707: }
1708: }
1709: if ( d )
1710: d->sugar = sugar;
1711: *rp = d;
1712: }
1713:
1714: void _dp_nf(b,g,ps,full,rp)
1715: NODE b;
1716: DP g;
1717: DP *ps;
1718: int full;
1719: DP *rp;
1720: {
1721: DP u,p,d,s,t;
1722: P dmy;
1723: NODE l;
1724: MP m,mr;
1725: int sugar,psugar;
1726:
1727: if ( !g ) {
1728: *rp = 0; return;
1729: }
1730: sugar = g->sugar;
1731: for ( d = 0; g; ) {
1732: for ( u = 0, l = b; l; l = NEXT(l) ) {
1733: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
1734: dp_load((int)BDY(l),&p);
1735: dp_red(d,g,p,&t,&u,&dmy);
1736: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1737: sugar = MAX(sugar,psugar);
1738: if ( !u ) {
1739: if ( d )
1740: d->sugar = sugar;
1741: *rp = d; return;
1742: }
1743: d = t;
1744: break;
1745: }
1746: }
1747: if ( u )
1748: g = u;
1749: else if ( !full ) {
1750: if ( g ) {
1751: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1752: }
1753: *rp = g; return;
1754: } else {
1755: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1756: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1757: addd(CO,d,t,&s); d = s;
1758: dp_rest(g,&t); g = t;
1759: }
1760: }
1761: if ( d )
1762: d->sugar = sugar;
1763: *rp = d;
1764: }
1765:
1766: #define SAFENM(q) ((q)?NM(q):0)
1767:
1768: double pz_t_e, pz_t_d, pz_t_d1, pz_t_c, im_t_s, im_t_r;
1769:
1770: void _dp_nf_ptozp(b,g,ps,full,multiple,r)
1771: NODE b;
1772: DP g;
1773: DP *ps;
1774: int full,multiple;
1775: DP *r;
1776: {
1777: DP u,dp,rp,t,t1,t2,red,shift;
1778: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
1779: N gn,tn,cn;
1780: NODE l;
1781: MP m,mr;
1782: int hmag,denom;
1783: int sugar,psugar;
1784: NODE dist;
1785: STRING imul;
1786: int ndist;
1787: int kara_bit;
1788: extern int kara_mag;
1789: double get_rtime();
1790: double t_0,t_00,tt,ttt,t_p,t_m,t_m1,t_m2,t_s,t_g,t_a;
1791:
1792: if ( !g ) {
1793: *r = 0; return;
1794: }
1795: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1796: t_p = t_m = t_m1 = t_m2 = t_s = t_g = t_a = 0;
1797:
1798: denom = Denominator?Denominator:1;
1799: hmag = multiple*HMAG(g)/denom;
1800: kara_bit = kara_mag*27; /* XXX */
1801: if ( Dist ) {
1802: dist = BDY(Dist);
1803: ndist = length(dist);
1804: }
1805: sugar = g->sugar;
1806:
1807: dc = 0; dp = 0; rc = ONE; rp = g;
1808: MKSTR(imul,"dp_imul_index");
1809:
1810: for ( ; rp; ) {
1811: for ( u = 0, l = b; l; l = NEXT(l) ) {
1812: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
1813: t_0 = get_rtime();
1814: dp_load((int)BDY(l),&red);
1815: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
1816: gcdn(NM(hr),NM(hred),&gn);
1817: divsn(NM(hred),gn,&tn); NTOQ(tn,SGN(hred),cr);
1818: divsn(NM(hr),gn,&tn); NTOQ(tn,SGN(hr),cred);
1819: tt = get_rtime()-t_0; t_p += tt;
1820:
1821: t_0 = get_rtime();
1822: dp_subd(rp,red,&shift);
1823: #if INET
1824: if ( Dist && ndist && HMAG(red) > kara_bit ) {
1825: NODE n0,n1,n2,n3;
1826: int i,s,id;
1827: Obj dmy;
1828: Q ind;
1829:
1830: if ( Print )
1831: fprintf(asir_out,"d");
1832: i = (int)BDY(l); STOQ(i,ind);
1833: chsgnp((P)cred,(P *)&mcred);
1834:
1835: MKNODE(n3,ind,0); MKNODE(n2,mcred,n3);
1836: MKNODE(n1,imul,n2); MKNODE(n0,BDY(dist),n1);
1837: Pox_rpc(n0,&dmy);
1838: muldc(CO,rp,(P)cr,&t);
1839: NEXT(n0)=0;
1840: Pox_pop_local(n0,&t1);
1841: } else
1842: #endif
1843: {
1844: /*
1845: if ( Print )
1846: fprintf(asir_out,"l");
1847: */
1848: t_00 = get_rtime();
1849: muldc(CO,rp,(P)cr,&t);
1850: ttt = get_rtime()-t_00; t_m1 += ttt/dp_nt(rp);
1851: t_00 = get_rtime();
1852: chsgnp((P)cred,(P *)&mcred);
1853: muldc(CO,red,(P)mcred,&t1);
1854: ttt = get_rtime()-t_00; t_m2 += ttt/dp_nt(red);
1855: }
1856: t_00 = get_rtime();
1.5 noro 1857: muld(CO,shift,t1,&t2);
1.1 noro 1858: addd(CO,t,t2,&u);
1859: tt = get_rtime(); t_m += tt-t_0;
1860: ttt = get_rtime(); t_s += ttt-t_00;
1861:
1862: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
1863: sugar = MAX(sugar,psugar);
1864: if ( !u ) {
1865: if ( dp )
1866: dp->sugar = sugar;
1867: *r = dp;
1868: goto final;
1869: }
1870: break;
1871: }
1872: }
1873: if ( u ) {
1874: if ( HMAG(u) > hmag ) {
1875: t_0 = get_rtime();
1876: #if INET
1877: if ( Dist && HMAG(u) > kara_bit ) {
1878: if ( Print )
1879: fprintf(asir_out,"D");
1880: dp_ptozp_d(dist,ndist,u,&t);
1881: } else {
1882: if ( Print )
1883: fprintf(asir_out,"L");
1884: dp_ptozp_d(0,0,u,&t);
1885: }
1886: #else
1887: dp_ptozp(u,&t); /* XXX */
1888: #endif
1889: tt = get_rtime()-t_0; t_g += tt;
1890: t_0 = get_rtime();
1891: divsn(NM((Q)BDY(u)->c),NM((Q)BDY(t)->c),&cn); NTOQ(cn,1,cont);
1892: if ( !dp_fcoeffs && Print ) {
1893: fprintf(asir_out,"(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1894: fflush(asir_out);
1895: }
1896: mulq(cr,dc,&dcq);
1897: mulq(cont,rc,&rcq);
1898: gcdn(SAFENM(dcq),SAFENM(rcq),&gn);
1899: divsn(SAFENM(dcq),gn,&tn); NTOQ(tn,SGN(dcq),dc);
1900: divsn(SAFENM(rcq),gn,&tn); NTOQ(tn,SGN(rcq),rc);
1901: tt = get_rtime()-t_0; t_a += tt;
1902: rp = t;
1903: hmag = multiple*HMAG(rp)/denom;
1904: } else {
1905: t_0 = get_rtime();
1906: mulq(cr,dc,&dcq); dc = dcq;
1907: tt = get_rtime()-t_0; t_a += tt;
1908: rp = u;
1909: }
1910: } else if ( !full ) {
1911: if ( rp ) {
1912: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
1913: }
1914: *r = rp;
1915: goto final;
1916: } else {
1917: t_0 = get_rtime();
1918: mulq((Q)BDY(rp)->c,rc,&c);
1919: gcdn(SAFENM(dc),SAFENM(c),&gn); NTOQ(gn,1,gcd);
1920: divsn(SAFENM(dc),gn,&tn); NTOQ(tn,SGN(dc),dcq);
1921: divsn(SAFENM(c),gn,&tn); NTOQ(tn,SGN(c),cq);
1922: muldc(CO,dp,(P)dcq,&t1);
1923: m = BDY(rp); NEWMP(mr); mr->dl = m->dl; mr->c = (P)cq;
1924: NEXT(mr) = 0; MKDP(rp->nv,mr,t); t->sugar = mr->dl->td;
1925: addd(CO,t,t1,&dp);
1926: dc = gcd;
1927:
1928: dp_rest(rp,&t); rp = t;
1929: tt = get_rtime()-t_0; t_a += tt;
1930: }
1931: }
1932: if ( dp )
1933: dp->sugar = sugar;
1934: *r = dp;
1935: final:
1936: if ( Print )
1937: fprintf(asir_out,"(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
1938: t_p,t_m,t_m1,t_m2,t_s,
1939: t_g,t_a,
1940: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
1941: }
1942:
1943: void dp_imul_d();
1944: void imulv();
1945:
1946: void _dp_nf_ptozp_mpi(b,g,ps,full,multiple,r)
1947: NODE b;
1948: DP g;
1949: DP *ps;
1950: int full,multiple;
1951: DP *r;
1952: {
1953: Obj dmy;
1954: DP u,dp,rp,t,t1,t2,red,shift;
1955: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
1956: N gn,tn,cn;
1957: NODE l,n0,n1,n2,n3;
1958: MP m,mr;
1959: int i,n;
1960: int hmag,denom;
1961: int sugar,psugar;
1962: NODE dist;
1963: STRING imul;
1964: int ndist;
1965: int kara_bit;
1966: extern int kara_mag;
1967: extern int mpi_mag;
1968: Q ind;
1969: int id;
1970: double get_rtime();
1971: double t_0,t_00,tt,ttt,t_p,t_m,t_m1,t_m2,t_s,t_g,t_a;
1972:
1973: if ( !g ) {
1974: *r = 0; return;
1975: }
1976: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1977: t_p = t_m = t_m1 = t_m2 = t_s = t_g = t_a = 0;
1978:
1979: denom = Denominator?Denominator:1;
1980: hmag = multiple*HMAG(g)/denom;
1981: kara_bit = kara_mag*27; /* XXX */
1982: if ( Dist ) {
1983: dist = BDY(Dist);
1984: ndist = length(dist);
1985: }
1986: sugar = g->sugar;
1987:
1988: dc = 0; dp = 0; rc = ONE; rp = g;
1989: MKSTR(imul,"dp_imul_index");
1990:
1991: for ( ; rp; ) {
1992: for ( u = 0, l = b; l; l = NEXT(l) ) {
1993: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
1994: t_0 = get_rtime();
1995: dp_load((int)BDY(l),&red);
1996: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
1997: gcdn(NM(hr),NM(hred),&gn);
1998: divsn(NM(hred),gn,&tn); NTOQ(tn,SGN(hred),cr);
1999: divsn(NM(hr),gn,&tn); NTOQ(tn,SGN(hr),cred);
2000: tt = get_rtime()-t_0; t_p += tt;
2001:
2002: t_0 = get_rtime();
2003: dp_subd(rp,red,&shift);
2004: t_00 = get_rtime();
2005: if ( Dist && ndist
2006: && HMAG(rp) > mpi_mag && p_mag((P)cr) > mpi_mag ) {
2007: if ( Print ) fprintf(asir_out,"~");
2008: dp_imul_d(rp,cr,&t);
2009: } else {
2010: if ( Print ) fprintf(asir_out,"_");
2011: muldc(CO,rp,(P)cr,&t);
2012: }
2013: ttt = get_rtime()-t_00; t_m1 += ttt/dp_nt(rp);
2014:
2015: t_00 = get_rtime();
2016: chsgnp((P)cred,(P *)&mcred);
2017: if ( Dist && ndist
2018: && HMAG(red) > mpi_mag && p_mag((P)mcred) > mpi_mag ) {
2019: if ( Print ) fprintf(asir_out,"=");
2020: dp_imul_d(red,mcred,&t1);
2021: } else {
2022: if ( Print ) fprintf(asir_out,"_");
2023: muldc(CO,red,(P)mcred,&t1);
2024: }
2025: ttt = get_rtime()-t_00; t_m2 += ttt/dp_nt(red);
2026:
2027: t_00 = get_rtime();
1.5 noro 2028: muld(CO,shift,t1,&t2);
1.1 noro 2029: addd(CO,t,t2,&u);
2030: tt = get_rtime(); t_m += tt-t_0;
2031: ttt = get_rtime(); t_s += ttt-t_00;
2032:
2033: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
2034: sugar = MAX(sugar,psugar);
2035: if ( !u ) {
2036: if ( dp )
2037: dp->sugar = sugar;
2038: *r = dp;
2039: goto final;
2040: }
2041: break;
2042: }
2043: }
2044: if ( u ) {
2045: if ( HMAG(u) > hmag ) {
2046: t_0 = get_rtime();
2047: if ( Dist && HMAG(u) > mpi_mag ) {
2048: if ( Print )
2049: fprintf(asir_out,"D");
2050: dp_ptozp_d(dist,ndist,u,&t);
2051: } else {
2052: if ( Print )
2053: fprintf(asir_out,"L");
2054: dp_ptozp_d(0,0,u,&t);
2055: }
2056: tt = get_rtime()-t_0; t_g += tt;
2057: t_0 = get_rtime();
2058: divsn(NM((Q)BDY(u)->c),NM((Q)BDY(t)->c),&cn); NTOQ(cn,1,cont);
2059: if ( Print ) {
2060: fprintf(asir_out,"(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
2061: fflush(asir_out);
2062: }
2063: mulq(cr,dc,&dcq);
2064: mulq(cont,rc,&rcq);
2065: gcdn(SAFENM(dcq),SAFENM(rcq),&gn);
2066: divsn(SAFENM(dcq),gn,&tn); NTOQ(tn,SGN(dcq),dc);
2067: divsn(SAFENM(rcq),gn,&tn); NTOQ(tn,SGN(rcq),rc);
2068: tt = get_rtime()-t_0; t_a += tt;
2069: rp = t;
2070: hmag = multiple*HMAG(rp)/denom;
2071: } else {
2072: t_0 = get_rtime();
2073: mulq(cr,dc,&dcq); dc = dcq;
2074: tt = get_rtime()-t_0; t_a += tt;
2075: rp = u;
2076: }
2077: } else if ( !full ) {
2078: if ( rp ) {
2079: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
2080: }
2081: *r = rp;
2082: goto final;
2083: } else {
2084: t_0 = get_rtime();
2085: mulq((Q)BDY(rp)->c,(Q)rc,(Q *)&c);
2086: gcdn(SAFENM(dc),SAFENM(c),&gn); NTOQ(gn,1,gcd);
2087: divsn(SAFENM(dc),gn,&tn); NTOQ(tn,SGN(dc),dcq);
2088: divsn(SAFENM(c),gn,&tn); NTOQ(tn,SGN(c),cq);
2089: muldc(CO,dp,(P)dcq,&t1);
2090: m = BDY(rp); NEWMP(mr); mr->dl = m->dl; mr->c = (P)cq;
2091: NEXT(mr) = 0; MKDP(rp->nv,mr,t); t->sugar = mr->dl->td;
2092: addd(CO,t,t1,&dp);
2093: dc = gcd;
2094:
2095: dp_rest(rp,&t); rp = t;
2096: tt = get_rtime()-t_0; t_a += tt;
2097: }
2098: }
2099: if ( dp )
2100: dp->sugar = sugar;
2101: *r = dp;
2102: final:
2103: if ( Print )
2104: fprintf(asir_out,"(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
2105: t_p,t_m,t_m1,t_m2,t_s,
2106: t_g,t_a,
2107: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
2108: }
2109:
2110: void dp_imul_d(p,q,rp)
2111: DP p;
2112: Q q;
2113: DP *rp;
2114: {
2115: int nsep,ndist,i,j,k,l,n;
2116: double t0,t1,t2;
2117: Q *s;
2118: pointer *b;
2119: VECT c,cs,ri;
2120: VECT *r;
2121: MP m;
2122: NODE tn,dist,n0,n1,n2;
2123: Obj dmy;
2124: STRING imul;
2125:
2126: extern LIST Dist;
2127:
2128: if ( !p || !q ) {
2129: *rp = 0; return;
2130: }
2131: dist = BDY(Dist);
2132: for ( tn = dist, ndist = 0; tn; tn = NEXT(tn), ndist++ );
2133: nsep = ndist + 1;
2134: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
2135: if ( n <= nsep ) {
2136: muldc(CO,p,(P)q,rp); return;
2137: }
2138: MKSTR(imul,"imulv");
2139: t0 = get_rtime();
2140: dp_dtov(p,&c);
2141: sepvect(c,nsep,&cs);
2142: r = (VECT *)CALLOC(nsep,sizeof(VECT *));
2143: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
2144: n0 = mknode(4,BDY(tn),imul,b[i],q);
2145: Pox_rpc(n0,&dmy);
2146: }
2147: t1 = get_rtime();
2148: im_t_s += t1 - t0;
2149: imulv(b[i],q,&r[i]);
2150: t1 = get_rtime();
2151: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
2152: MKNODE(n0,BDY(tn),0);
2153: Pox_pop_local(n0,&r[i]);
2154: if ( OID(r[i]) == O_ERR ) {
2155: printexpr(CO,(Obj)r[i]);
2156: error("dp_imul_d : aborted");
2157: }
2158: }
2159: t2 = get_rtime();
2160: im_t_r += t2 - t1;
2161: s = (Q *)CALLOC(n,sizeof(Q));
2162: for ( i = j = 0; i < nsep; i++ ) {
2163: for ( k = 0, ri = r[i], l = ri->len; k < l; k++, j++ ) {
2164: s[j] = (Q)BDY(ri)[k];
2165: }
2166: }
2167: dp_vtod(s,p,rp);
2168: }
2169:
2170: void imulv(w,c,rp)
2171: VECT w;
2172: Q c;
2173: VECT *rp;
2174: {
2175: int n,i;
2176: VECT r;
2177:
2178: n = w->len;
2179: MKVECT(r,n); *rp = r;
2180: for ( i = 0; i < n; i++ )
2181: mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
2182: }
2183:
2184: DP_pairs updpairs( d, g, t)
2185: DP_pairs d;
2186: NODE /* of index */ g;
2187: int t;
2188: {
2189: register DP_pairs d1, dd, nd;
2190: struct oEGT tup0,tup1;
2191: int dl,dl1;
2192:
2193: if ( !g ) return d;
2194: get_eg(&tup0);
2195: if ( !NoCriB && d ) {
2196: dl = DPPlength(d);
2197: d = criterion_B( d, t );
2198: dl -= DPPlength(d); NBP += dl;
2199: }
2200: d1 = newpairs( g, t );
2201: if ( NEXT(d1) ) {
2202: dl = DPPlength(d1); TP += dl;
2203: d1 = criterion_M( d1 );
2204: dl1 = DPPlength(d1); NMP += (dl-dl1); dl = dl1;
2205: d1 = criterion_F( d1 );
2206: dl1 = DPPlength(d1); NFP += (dl-dl1); dl = dl1;
2207: } else
2208: dl = 1;
1.5 noro 2209: if ( !do_weyl )
2210: for ( dd = 0; d1; d1 = nd ) {
2211: nd = NEXT(d1);
2212: if ( !criterion_2( d1->dp1, d1->dp2 ) ) {
2213: NEXT(d1) = dd;
2214: dd = d1;
2215: }
1.1 noro 2216: }
1.5 noro 2217: else
2218: dd = d1;
1.1 noro 2219: dl1 = DPPlength(dd); NDP += (dl-dl1);
2220: get_eg(&tup1);
2221: add_eg(&eg_up,&tup0,&tup1);
2222: if ( !(nd = d) ) return dd;
2223: while ( nd = NEXT(d1 = nd) ) ;
2224: NEXT(d1) = dd;
2225: return d;
2226: }
2227:
2228: DP_pairs newpairs( g, t )
2229: NODE /* of index */ g;
2230: register int t;
2231: {
2232: register NODE r;
2233: register DL tdl = psh[t];
2234: register int ts;
2235: register DP_pairs p, last;
2236: int dp;
2237: register DL dl;
2238: register int s;
2239:
2240: ts = pss[t] - tdl->td;
2241: for ( last = 0, r = g; r; r = NEXT(r) ) {
2242: NEXT(p = NEWDP_pairs) = last;
2243: last = p;
2244: dp = p->dp1 = (int)BDY(r); p->dp2 = t;
2245: p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
1.5 noro 2246: #if 0
2247: if ( do_weyl )
2248: p->sugar = dl_weight(p->lcm);
2249: else
2250: #endif
2251: p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
1.1 noro 2252: }
2253: return last;
2254: }
2255:
2256: DP_pairs criterion_B( d, s )
2257: register DP_pairs d;
2258: int s;
2259: {
2260: register DP_pairs dd, p;
2261: register DL tij, t = psh[s], dltmp;
2262:
2263: if ( !d ) return 0;
2264: NEWDL( dltmp, CNVars );
2265: for ( dd = 0; d; d = p ) {
2266: p = NEXT(d),
2267: tij = d->lcm;
2268: if ( tij->td != lcm_of_DL(CNVars, tij, t, dltmp )->td
2269: || !dl_equal(CNVars, tij, dltmp )
2270: || (tij->td == lcm_of_DL(CNVars, psh[d->dp1], t, dltmp )->td
2271: && dl_equal(CNVars, dltmp, tij ))
2272: || (tij->td == lcm_of_DL(CNVars, psh[d->dp2], t, dltmp )->td
2273: && dl_equal(CNVars, dltmp, tij )) ) {
2274: NEXT(d) = dd;
2275: dd = d;
2276: }
2277: }
2278: return dd;
2279: }
2280:
2281: DP_pairs criterion_M( d1 )
2282: DP_pairs d1;
2283: {
2284: register DP_pairs dd, e, d3, d2, p;
2285: register DL itdl, jtdl;
2286: register int itdltd, jtdltd;
2287:
2288: for ( dd = 0, e = d1; e; e = d3 ) {
2289: if ( !(d2 = NEXT(e)) ) {
2290: NEXT(e) = dd;
2291: return e;
2292: }
2293: itdltd = (itdl = e->lcm)->td;
2294: for ( d3 = 0; d2; d2 = p ) {
2295: p = NEXT(d2),
2296: jtdltd = (jtdl = d2->lcm)->td;
2297: if ( jtdltd == itdltd )
2298: if ( dl_equal(CNVars, itdl, jtdl ) ) ;
2299: else if ( dl_redble( jtdl, itdl ) ) continue;
2300: else if ( dl_redble( itdl, jtdl ) ) goto delit;
2301: else ;
2302: else if ( jtdltd > itdltd )
2303: if ( dl_redble( jtdl, itdl ) ) continue;
2304: else ;
2305: else if ( dl_redble( itdl, jtdl ) ) goto delit;
2306: NEXT(d2) = d3;
2307: d3 = d2;
2308: }
2309: NEXT(e) = dd;
2310: dd = e;
2311: continue;
2312: /**/
2313: delit: NEXT(d2) = d3;
2314: d3 = d2;
2315: for ( ; p; p = d2 ) {
2316: d2 = NEXT(p);
2317: NEXT(p) = d3;
2318: d3 = p;
2319: }
2320: }
2321: return dd;
2322: }
2323:
2324: static DP_pairs collect_pairs_of_hdlcm( d1, prest )
2325: DP_pairs d1, *prest;
2326: {
2327: register DP_pairs w, p, r, s;
2328: register DL ti;
2329: register int td;
2330:
2331: td = (ti = (w = d1)->lcm)->td;
2332: s = NEXT(w);
2333: NEXT(w) = 0;
2334: for ( r = 0; s; s = p ) {
2335: p = NEXT(s);
2336: if ( td == s->lcm->td && dl_equal(CNVars, ti, s->lcm ) )
2337: {
2338: NEXT(s) = w;
2339: w = s;
2340: } else {
2341: NEXT(s) = r;
2342: r = s;
2343: }
2344: }
2345: *prest = r;
2346: return w;
2347: }
2348:
2349: int criterion_2( dp1, dp2 )
2350: int dp1, dp2;
2351: {
2352: register int i, *d1, *d2;
2353:
2354: d1 = psh[dp1]->d, d2 = psh[dp2]->d;
2355: for ( i = CNVars; --i >= 0; d1++, d2++ )
2356: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
2357: return 1;
2358: }
2359:
2360: DP_pairs criterion_F( d1 )
2361: DP_pairs d1;
2362: {
2363: DP_pairs rest, head;
2364: register DP_pairs last, p, r, w;
2365: register int s;
2366:
2367: for ( head = last = 0, p = d1; NEXT(p); ) {
2368: s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
2369: while ( w = NEXT(w) )
2370: if ( criterion_2( w->dp1, w->dp2 ) ) {
2371: r = w;
2372: break;
2373: } else if ( w->sugar < s ) s = (r = w)->sugar;
2374: if ( last ) NEXT(last) = r;
2375: else head = r;
2376: NEXT(last = r) = 0;
2377: if ( !(p = rest) ) return head;
2378: }
2379: if ( !last ) return p;
2380: NEXT(last) = p;
2381: return head;
2382: }
2383:
2384: NODE updbase(g,t)
2385: NODE g;
2386: int t;
2387: {
2388: g = remove_reducibles(g,t);
2389: g = append_one(g,t);
2390: return g;
2391: }
2392:
2393: NODE /* of index */ remove_reducibles( nd, newdp )
2394: NODE /* of index */ nd;
2395: int newdp;
2396: {
2397: register DL dl, dln;
2398: register NODE last, p, head;
2399: register int td;
2400:
2401: dl = psh[newdp];
2402: td = dl->td;
2403: for ( head = last = 0, p = nd; p; ) {
2404: dln = psh[(int)BDY(p)];
2405: if ( dln->td >= td && dl_redble( dln, dl ) ) {
2406: p = NEXT(p);
2407: if ( last ) NEXT(last) = p;
2408: } else {
2409: if ( !last ) head = p;
2410: p = NEXT(last = p);
2411: }
2412: }
2413: return head;
2414: }
2415:
2416: int dl_redble( dl1, dl2 )
2417: DL dl1, dl2;
2418: {
2419: register int n, *d1, *d2;
2420:
2421: for ( d1 = dl1->d, d2 = dl2->d, n = CNVars; --n >= 0; d1++, d2++ )
2422: if ( *d1 < *d2 ) return 0;
2423: return 1;
1.5 noro 2424: }
2425:
2426: int dl_weight(dl)
2427: DL dl;
2428: {
2429: int n,w,i;
2430:
2431: n = CNVars/2;
2432: for ( i = 0, w = 0; i < n; i++ )
2433: w += (-dl->d[i]+dl->d[n+i]);
2434: return w;
1.1 noro 2435: }
2436:
2437: int gbcheck(f)
2438: NODE f;
2439: {
2440: int i;
2441: NODE r,g,gall;
2442: DP_pairs d,l;
2443: DP h,nf,dp1,dp2;
2444: struct oEGT tmp0,tmp1;
2445:
2446: if ( NoGC )
2447: return 1;
2448: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
2449: i = (int)BDY(r);
2450: d = updpairs(d,g,i);
2451: g = updbase(g,i);
2452: gall = append_one(gall,i);
2453: }
2454: if ( Print || PrintShort ) {
2455: fprintf(asir_out,"gbcheck total %d pairs\n",DPPlength(d)); fflush(asir_out);
2456: }
2457: while ( d ) {
2458: l = d; d = NEXT(d);
2459: get_eg(&tmp0);
2460: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2); dp_sp(dp1,dp2,&h);
2461: _dp_nf(gall,h,ps,1,&nf);
2462: get_eg(&tmp1); add_eg(&eg_gc,&tmp0,&tmp1);
2463: if ( Print || PrintShort ) {
2464: fprintf(asir_out,"."); fflush(asir_out);
2465: }
2466: if ( nf )
2467: return 0;
2468: }
2469: if ( Print || PrintShort )
2470: fprintf(asir_out,"\n");
2471: return 1;
2472: }
2473:
2474: int membercheck(f,x)
2475: NODE f,x;
2476: {
2477: DP g;
2478: struct oEGT tmp0,tmp1;
2479:
2480: if ( NoMC )
2481: return 1;
2482: if ( Print || PrintShort ) {
2483: fprintf(asir_out,"membercheck\n"); fflush(asir_out);
2484: }
2485: for ( ; f; f = NEXT(f) ) {
2486: get_eg(&tmp0);
2487: _dp_nf(x,(DP)BDY(f),ps,1,&g);
2488: get_eg(&tmp1); add_eg(&eg_mc,&tmp0,&tmp1);
2489: if ( Print ) {
2490: print_split_eg(&tmp0,&tmp1); fflush(asir_out);
2491: } else if ( PrintShort ) {
2492: fprintf(asir_out,"."); fflush(asir_out);
2493: }
2494: if ( g )
2495: return 0;
2496: }
2497: if ( Print || PrintShort )
2498: fprintf(asir_out,"\n");
2499: return 1;
2500: }
2501:
2502: void dp_set_flag(name,value)
2503: Obj name,value;
2504: {
2505: char *n;
2506: int v;
2507:
2508: if ( OID(name) != O_STR )
2509: return;
2510: n = BDY((STRING)name);
2511: if ( !strcmp(n,"Demand") ) {
2512: Demand = value ? BDY((STRING)value) : 0; return;
2513: }
2514: if ( !strcmp(n,"Dist") ) {
2515: Dist = (LIST)value; return;
2516: }
2517: if ( value && OID(value) != O_N )
2518: return;
2519: v = QTOS((Q)value);
2520: if ( !strcmp(n,"NoSugar") )
2521: NoSugar = v;
2522: else if ( !strcmp(n,"NoCriB") )
2523: NoCriB = v;
2524: else if ( !strcmp(n,"NoGC") )
2525: NoGC = v;
2526: else if ( !strcmp(n,"NoMC") )
2527: NoMC = v;
2528: else if ( !strcmp(n,"NoRA") )
2529: NoRA = v;
2530: else if ( !strcmp(n,"NoGCD") )
2531: NoGCD = v;
2532: else if ( !strcmp(n,"Top") )
2533: Top = v;
2534: else if ( !strcmp(n,"ShowMag") )
2535: ShowMag = v;
2536: else if ( !strcmp(n,"PrintShort") )
2537: PrintShort = v;
2538: else if ( !strcmp(n,"Print") )
2539: Print = v;
2540: else if ( !strcmp(n,"Stat") )
2541: Stat = v;
2542: else if ( !strcmp(n,"Reverse") )
2543: Reverse = v;
2544: else if ( !strcmp(n,"InterReduce") )
2545: InterReduce = v;
2546: else if ( !strcmp(n,"Multiple") )
2547: Multiple = v;
2548: else if ( !strcmp(n,"Denominator") )
2549: Denominator = v;
2550: else if ( !strcmp(n,"PtozpRA") )
2551: PtozpRA = v;
2552: }
2553:
2554: void dp_make_flaglist(list)
2555: LIST *list;
2556: {
2557: Q v;
2558: STRING name,path;
2559: NODE n,n1;
2560:
2561: STOQ(Multiple,v); MKNODE(n,v,0); MKSTR(name,"Multiple"); MKNODE(n1,name,n); n = n1;
2562: STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
2563: MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
2564: STOQ(InterReduce,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"InterReduce"); MKNODE(n1,name,n); n = n1;
2565: STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
2566: STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
2567: STOQ(Print,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Print"); MKNODE(n1,name,n); n = n1;
2568: STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
2569: STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
2570: STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
2571: STOQ(NoGCD,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGCD"); MKNODE(n1,name,n); n = n1;
2572: STOQ(NoRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoRA"); MKNODE(n1,name,n); n = n1;
2573: STOQ(NoMC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoMC"); MKNODE(n1,name,n); n = n1;
2574: STOQ(NoGC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGC"); MKNODE(n1,name,n); n = n1;
2575: STOQ(NoCriB,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoCriB"); MKNODE(n1,name,n); n = n1;
2576: STOQ(NoSugar,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoSugar"); MKNODE(n1,name,n); n = n1;
2577: if ( Demand )
2578: MKSTR(path,Demand);
2579: else
2580: path = 0;
2581: MKNODE(n1,path,n); n = n1; MKSTR(name,"Demand"); MKNODE(n1,name,n); n = n1;
2582: MKLIST(*list,n);
2583: }
2584:
2585: #define DELIM '/'
2586:
2587: #if defined(THINK_C)
2588: #undef DELIM
2589: #define DELIM ':'
2590: #endif
2591:
2592: void dp_save(index,p)
2593: int index;
2594: DP p;
2595: {
2596: FILE *fp;
2597: char path[BUFSIZ];
2598:
2599: sprintf(path,"%s%c%d",Demand,DELIM,index);
2600: #if defined(VISUAL) || defined(THINK_C)
2601: if ( !(fp = fopen(path,"wb") ) )
2602: #else
2603: if ( !(fp = fopen(path,"w") ) )
2604: #endif
2605: error("dp_save : cannot open a file");
2606: savevl(fp,VC); saveobj(fp,(Obj)p); fclose(fp);
2607: }
2608:
2609: void dp_load(index,p)
2610: int index;
2611: DP *p;
2612: {
2613: FILE *fp;
2614: char path[BUFSIZ];
2615:
2616: if ( !Demand || ps[index] )
2617: *p = ps[index];
2618: else {
2619: sprintf(path,"%s%c%d",Demand,DELIM,index);
2620: #if defined(VISUAL) || defined(THINK_C)
2621: if ( !(fp = fopen(path,"rb") ) )
2622: #else
2623: if ( !(fp = fopen(path,"r") ) )
2624: #endif
2625: error("dp_load : cannot open a file");
2626: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp);
2627: }
2628: }
2629:
2630: int dp_load_t(index,p)
2631: int index;
2632: DP *p;
2633: {
2634: FILE *fp;
2635: char path[BUFSIZ];
2636:
2637: sprintf(path,"%s%c%d",Demand,DELIM,index);
2638: #if defined(VISUAL) || defined(THINK_C)
2639: if ( !(fp = fopen(path,"rb") ) )
2640: #else
2641: if ( !(fp = fopen(path,"r") ) )
2642: #endif
2643: return 0;
2644: else {
2645: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp); return 1;
2646: }
2647: }
2648:
2649: void init_stat() {
2650: init_eg(&eg_up); init_eg(&eg_sp); init_eg(&eg_spm);
2651: init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
2652: init_eg(&eg_pz); init_eg(&eg_np); init_eg(&eg_mp);
2653: init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
2654: ZR = NZR = TP = NBP = NMP = NFP = NDP = 0;
2655: }
2656:
2657: void print_stat() {
2658: if ( !Print && !Stat )
2659: return;
2660: print_eg("UP",&eg_up); print_eg("SP",&eg_sp); print_eg("SPM",&eg_spm);
2661: print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
2662: print_eg("PZ",&eg_pz); print_eg("NP",&eg_np); print_eg("MP",&eg_mp);
2663: print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
2664: 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);
2665: }
2666:
2667: /*
2668: * Old codes
2669: */
2670:
2671: #if 0
2672: void _dp_nf_ptozp(b,g,ps,full,multiple,rp)
2673: NODE b;
2674: DP g;
2675: DP *ps;
2676: int full,multiple;
2677: DP *rp;
2678: {
2679: DP u,p,d,s,t,dmy;
2680: NODE l;
2681: MP m,mr;
2682: int i,n;
2683: int hmag,denom;
2684: int sugar,psugar;
2685: NODE dist;
2686: int ndist;
2687: int kara_bit;
2688: extern int kara_mag;
2689:
2690: if ( !g ) {
2691: *rp = 0; return;
2692: }
2693: denom = Denominator?Denominator:1;
2694: hmag = multiple*HMAG(g)/denom;
2695: kara_bit = kara_mag*27; /* XXX */
2696: if ( Dist ) {
2697: dist = BDY(Dist);
2698: ndist = length(dist);
2699: }
2700: sugar = g->sugar;
2701: for ( d = 0; g; ) {
2702: for ( u = 0, l = b; l; l = NEXT(l) ) {
2703: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
2704: dp_load((int)BDY(l),&p);
2705: dp_red(d,g,p,&t,&u,&dmy);
2706: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
2707: sugar = MAX(sugar,psugar);
2708: if ( !u ) {
2709: if ( d )
2710: d->sugar = sugar;
2711: *rp = d; return;
2712: }
2713: d = t;
2714: break;
2715: }
2716: }
2717: if ( u ) {
2718: g = u;
2719: if ( d ) {
2720: if ( HMAG(d) > hmag ) {
2721: if ( Dist && HMAG(g) > kara_bit )
2722: dp_ptozp2_d(dist,ndist,d,g,&t,&u);
2723: else
2724: dp_ptozp2(d,g,&t,&u);
2725: d = t; g = u;
2726: hmag = multiple*HMAG(d)/denom;
2727: }
2728: } else {
2729: if ( HMAG(g) > hmag ) {
2730: #if INET
2731: if ( Dist && HMAG(g) > kara_bit ) {
2732: dp_ptozp_d(dist,ndist,g,&t);
2733: } else
2734: #endif
2735: dp_ptozp(g,&t);
2736: g = t;
2737: hmag = multiple*HMAG(g)/denom;
2738: }
2739: }
2740: }
2741: else if ( !full ) {
2742: if ( g ) {
2743: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
2744: }
2745: *rp = g; return;
2746: } else {
2747: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
2748: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
2749: addd(CO,d,t,&s); d = s;
2750: dp_rest(g,&t); g = t;
2751: }
2752: }
2753: if ( d )
2754: d->sugar = sugar;
2755: *rp = d;
2756: }
2757:
2758: int criterion_2(dp1,dp2)
2759: {
2760: DP g1,g2,g,s1,s2;
2761:
2762: monomial_gcd(ps[dp1],&g1); monomial_gcd(ps[dp2],&g2);
2763: dp_gcd(g1,g2,&g);
2764: dp_subd(ps[dp1],g,&s1); dp_subd(ps[dp2],g,&s2);
2765: return _criterion_2(s1,s2);
2766: }
2767:
2768: int _criterion_2( dp1, dp2 )
2769: DP dp1, dp2;
2770: {
2771: register int i, *d1, *d2;
2772:
2773: d1 = BDY(dp1)->dl->d; d2 = BDY(dp2)->dl->d;
2774: for ( i = CNVars; --i >= 0; d1++, d2++ )
2775: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
2776: return 1;
2777: }
2778:
2779: void dp_gcd(p1,p2,rp)
2780: DP p1,p2;
2781: DP *rp;
2782: {
2783: int i,n,td;
2784: DL d1,d2,d;
2785: MP m;
2786: DP s;
2787:
2788: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
2789: NEWDL(d,n);
2790: for ( i = 0, td = 0; i < n; i++ ) {
2791: d->d[i] = MIN(d1->d[i],d2->d[i]);
2792: td += d->d[i];
2793: }
2794: d->td = td;
2795: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
2796: *rp = s;
2797: }
2798:
2799: void monomial_gcd(p,rp)
2800: DP p;
2801: DP *rp;
2802: {
2803: int n,i,td;
2804: DL d,d1;
2805: MP m;
2806: DP s;
2807:
2808: n = p->nv; m = BDY(p); d = m->dl;
2809: NEWDL(d1,n);
2810: for ( i = 0; i < n; i++ )
2811: d1->d[i] = d->d[i];
2812: for ( m = NEXT(m); m; m = NEXT(m) ) {
2813: d = m->dl;
2814: for ( i = 0; i < n; i++ )
2815: d1->d[i] = MIN(d1->d[i],d->d[i]);
2816: }
2817: for ( i = 0, td = 0; i < n; i++ )
2818: td += d1->d[i];
2819: NEWMP(m); m->dl = d1; m->c = (P)ONE; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
2820: *rp = s;
2821: }
2822: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>