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