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