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