![[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.15 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.14 2000/12/08 02:39:05 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:
1.14 noro 77: struct oEGT eg_nf,eg_nfm;
78: struct oEGT eg_znfm,eg_pz,eg_np,eg_ra,eg_mc,eg_gc;
1.1 noro 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.14 noro 88: void dp_imul_d(DP,Q,DP *);
1.1 noro 89: void print_stat(void);
90: void init_stat(void);
91: int dp_load_t(int,DP *);
92: void dp_load(int,DP *);
1.7 noro 93: void dp_save(int,Obj,char *);
1.1 noro 94: void dp_make_flaglist(LIST *);
95: void dp_set_flag(Obj,Obj);
96: int membercheck(NODE,NODE);
97: int gbcheck(NODE);
98: int dl_redble(DL,DL);
99: NODE remove_reducibles(NODE,int);
100: NODE updbase(NODE,int);
101: DP_pairs criterion_F(DP_pairs);
102: int criterion_2(int,int);
103: static DP_pairs collect_pairs_of_hdlcm(DP_pairs,DP_pairs *);
104: DP_pairs criterion_M(DP_pairs);
105: DP_pairs criterion_B(DP_pairs,int);
106: DP_pairs newpairs(NODE,int);
107: DP_pairs updpairs(DP_pairs,NODE,int);
1.14 noro 108: void _dp_nf(NODE,DP,DP *,int,int,DP *);
1.1 noro 109: NODE gb_mod(NODE,int);
110: NODE gbd(NODE,int,NODE,NODE);
111: NODE gb(NODE,int,NODE);
112: NODE gb_f4(NODE);
113: NODE gb_f4_mod(NODE,int);
114: DP_pairs minp(DP_pairs, DP_pairs *);
115: void minsugar(DP_pairs,DP_pairs *,DP_pairs *);
116: NODE append_one(NODE,int);
117: void reducebase_dehomo(NODE,NODE *);
118: int newps_mod(DP,int);
119: int newps_nosave(DP,int,NODE);
120: int newps(DP,int,NODE);
121: void reduceall_mod(NODE,int,NODE *);
122: void reduceall(NODE,NODE *);
123: NODE NODE_sortbi(NODE,int);
124: NODE NODE_sortbi_insert(int, NODE,int);
125: NODE NODE_sortb(NODE,int);
126: NODE NODE_sortb_insert(DP,NODE,int);
127: void prim_part(DP,int,DP *);
128: void setup_arrays(NODE,int,NODE *);
129: int validhc(P,int,NODE);
130: void vlminus(VL,VL,VL *);
131: void printsubst(NODE);
132: void makesubst(VL,NODE *);
133: void pltovl(LIST,VL *);
134: void printdl(DL);
135: int DPPlength(DP_pairs);
136: void dp_gr_mod_main(LIST,LIST,Num,int,struct order_spec *,LIST *);
137: void dp_gr_main(LIST,LIST,Num,int,struct order_spec *,LIST *);
138: void dp_f4_main(LIST,LIST,struct order_spec *,LIST *);
139: void dp_f4_mod_main(LIST,LIST,int,struct order_spec *,LIST *);
140: double get_rtime();
141: void _dpmod_to_vect(DP,DL *,int *);
142: void dp_to_vect(DP,DL *,Q *);
1.12 noro 143: NODE dp_dllist(DP f);
144: NODE symb_merge(NODE,NODE,int),_symb_merge(NODE,NODE,int);
1.1 noro 145: extern int dp_nelim;
146: extern int dp_fcoeffs;
147: static DP *ps,*psm;
148: static DL *psh;
149: static P *psc;
150:
151: static int *pss;
152: static int psn,pslen;
153: static int NVars,CNVars,PCoeffs;
154: static VL VC;
1.14 noro 155:
156: int DP_Print = 0;
157: int DP_Multiple = 0;
158: LIST Dist = 0;
159: int NoGCD = 0;
160: int GenTrace = 0;
161: int OXCheck = -1;
162:
163: static DP_NFStat = 0;
1.1 noro 164: static int NoSugar = 0;
165: static int NoCriB = 0;
166: static int NoGC = 0;
167: static int NoMC = 0;
168: static int NoRA = 0;
1.12 noro 169: static int DP_PrintShort = 0;
1.1 noro 170: static int ShowMag = 0;
171: static int Stat = 0;
172: static int Denominator = 1;
173: static int Top = 0;
174: static int Reverse = 0;
175: static int Max_mag = 0;
176: static char *Demand = 0;
177: static int PtozpRA = 0;
1.14 noro 178:
1.1 noro 179: int doing_f4;
1.7 noro 180: NODE TraceList;
1.1 noro 181:
182: int eqdl(nv,dl1,dl2)
183: int nv;
184: DL dl1,dl2;
185: {
186: int i;
187: int *b1,*b2;
188:
189: if ( dl1->td != dl2->td )
190: return 0;
191: for ( i = 0, b1 = dl1->d, b2 = dl2->d; i < nv; i++ )
192: if ( b1[i] != b2[i] )
193: break;
194: if ( i == nv )
195: return 1;
196: else
197: return 0;
198: }
199:
1.5 noro 200: /* b[] should be cleared */
201:
1.1 noro 202: void _dpmod_to_vect(f,at,b)
203: DP f;
204: DL *at;
205: int *b;
206: {
207: int i,nv;
208: MP m;
209:
210: nv = f->nv;
211: for ( m = BDY(f), i = 0; m; m = NEXT(m), i++ ) {
212: for ( ; !eqdl(nv,m->dl,at[i]); i++ );
213: b[i] = ITOS(m->c);
214: }
215: }
216:
217: void dp_to_vect(f,at,b)
218: DP f;
219: DL *at;
220: Q *b;
221: {
222: int i,nv;
223: MP m;
224:
225: nv = f->nv;
226: for ( m = BDY(f), i = 0; m; m = NEXT(m), i++ ) {
227: for ( ; !eqdl(nv,m->dl,at[i]); i++ );
228: b[i] =(Q)m->c;
229: }
230: }
231:
1.2 noro 232: NODE dp_dllist(f)
1.1 noro 233: DP f;
234: {
1.2 noro 235: MP m;
236: NODE mp,mp0;
1.1 noro 237:
1.2 noro 238: if ( !f )
239: return 0;
1.1 noro 240: mp0 = 0;
241: for ( m = BDY(f); m; m = NEXT(m) ) {
1.2 noro 242: NEXTNODE(mp0,mp); BDY(mp) = (pointer)m->dl;
1.1 noro 243: }
244: NEXT(mp) = 0;
1.2 noro 245: return mp0;
246: }
247:
248: void pdl(f)
249: NODE f;
250: {
251: while ( f ) {
252: printdl(BDY(f)); f = NEXT(f);
253: }
254: fflush(stdout);
255: printf("\n");
1.1 noro 256: }
257:
258: void dp_gr_main(f,v,homo,modular,ord,rp)
259: LIST f,v;
260: Num homo;
261: int modular;
262: struct order_spec *ord;
263: LIST *rp;
264: {
265: int i,mindex,m,nochk;
266: struct order_spec ord1;
267: VL fv,vv,vc;
268: NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
269:
270: mindex = 0; nochk = 0; dp_fcoeffs = 0;
271: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
272: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
273: CNVars = homo ? NVars+1 : NVars;
274: if ( ord->id && NVars != ord->nv )
275: error("dp_gr_main : invalid order specification");
276: initd(ord);
277: if ( homo ) {
278: homogenize_order(ord,NVars,&ord1);
279: for ( fd0 = fi0 = 0, t = BDY(f); t; t = NEXT(t) ) {
280: NEXTNODE(fd0,fd); NEXTNODE(fi0,fi);
281: ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fi)); dp_homo((DP)BDY(fi),(DP *)&BDY(fd));
282: }
283: if ( fd0 ) NEXT(fd) = 0;
284: if ( fi0 ) NEXT(fi) = 0;
285: initd(&ord1);
286: } else {
287: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
288: NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
289: }
290: if ( fd0 ) NEXT(fd) = 0;
291: fi0 = fd0;
292: }
293: if ( modular < 0 ) {
294: modular = -modular; nochk = 1;
295: }
296: if ( modular )
297: m = modular > 1 ? modular : lprime[mindex];
298: else
299: m = 0;
300: makesubst(vc,&subst);
301: setup_arrays(fd0,0,&s);
302: init_stat();
303: while ( 1 ) {
304: if ( homo ) {
305: initd(&ord1); CNVars = NVars+1;
306: }
1.12 noro 307: if ( DP_Print && modular ) {
1.1 noro 308: fprintf(asir_out,"mod= %d, eval = ",m); printsubst(subst);
309: }
310: x = gb(s,m,subst);
311: if ( x ) {
312: if ( homo ) {
313: reducebase_dehomo(x,&xx); x = xx;
314: initd(ord); CNVars = NVars;
315: }
316: reduceall(x,&xx); x = xx;
317: if ( modular ) {
318: if ( nochk || (membercheck(fi0,x) && gbcheck(x)) )
319: break;
320: } else
321: break;
322: }
323: if ( modular )
324: if ( modular > 1 ) {
325: *rp = 0; return;
326: } else
327: m = lprime[++mindex];
328: makesubst(vc,&subst);
329: psn = length(s);
330: for ( i = psn; i < pslen; i++ ) {
331: pss[i] = 0; psh[i] = 0; psc[i] = 0; ps[i] = 0;
332: }
333: }
334: for ( r0 = 0; x; x = NEXT(x) ) {
335: NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
336: dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
337: }
338: if ( r0 ) NEXT(r) = 0;
339: MKLIST(*rp,r0);
340: print_stat();
341: if ( ShowMag )
342: fprintf(asir_out,"\nMax_mag=%d\n",Max_mag);
343: }
344:
345: void dp_gr_mod_main(f,v,homo,m,ord,rp)
346: LIST f,v;
347: Num homo;
348: int m;
349: struct order_spec *ord;
350: LIST *rp;
351: {
352: struct order_spec ord1;
353: VL fv,vv,vc;
354: NODE fd,fd0,r,r0,t,x,s,xx;
355: DP a,b,c;
356:
357: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
358: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
359: CNVars = homo ? NVars+1 : NVars;
360: if ( ord->id && NVars != ord->nv )
361: error("dp_gr_mod_main : invalid order specification");
362: initd(ord);
363: if ( homo ) {
364: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
365: ptod(CO,vv,(P)BDY(t),&a); dp_homo(a,&b);
366: if ( PCoeffs )
367: dp_mod(b,m,0,&c);
368: else
369: _dp_mod(b,m,(NODE)0,&c);
370: if ( c ) {
371: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
372: }
373: }
374: homogenize_order(ord,NVars,&ord1); initd(&ord1);
375: } else {
376: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
377: ptod(CO,vv,(P)BDY(t),&b);
378: if ( PCoeffs )
379: dp_mod(b,m,0,&c);
380: else
381: _dp_mod(b,m,0,&c);
382: if ( c ) {
383: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
384: }
385: }
386: }
387: if ( fd0 ) NEXT(fd) = 0;
388: setup_arrays(fd0,m,&s);
389: init_stat();
390: if ( homo ) {
391: initd(&ord1); CNVars = NVars+1;
392: }
393: x = gb_mod(s,m);
394: if ( homo ) {
395: reducebase_dehomo(x,&xx); x = xx;
396: initd(ord); CNVars = NVars;
397: }
398: reduceall_mod(x,m,&xx); x = xx;
399: if ( PCoeffs )
400: for ( r0 = 0; x; x = NEXT(x) ) {
401: NEXTNODE(r0,r); mdtop(CO,m,vv,ps[(int)BDY(x)],(P *)&BDY(r));
402: }
403: else
404: for ( r0 = 0; x; x = NEXT(x) ) {
405: NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
406: }
407: print_stat();
408: if ( r0 ) NEXT(r) = 0;
409: MKLIST(*rp,r0);
410: }
411:
412: void dp_f4_main(f,v,ord,rp)
413: LIST f,v;
414: struct order_spec *ord;
415: LIST *rp;
416: {
417: int i,mindex,m,nochk;
418: struct order_spec ord1;
419: VL fv,vv,vc;
420: NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
421:
422: dp_fcoeffs = 0;
423: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
424: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
425: CNVars = NVars;
426: if ( ord->id && NVars != ord->nv )
427: error("dp_f4_main : invalid order specification");
428: initd(ord);
429: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
430: NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
431: }
432: if ( fd0 ) NEXT(fd) = 0;
433: setup_arrays(fd0,0,&s);
434: x = gb_f4(s);
435: reduceall(x,&xx); x = xx;
436: for ( r0 = 0; x; x = NEXT(x) ) {
437: NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
438: dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
439: }
440: if ( r0 ) NEXT(r) = 0;
441: MKLIST(*rp,r0);
442: }
443:
444: void dp_f4_mod_main(f,v,m,ord,rp)
445: LIST f,v;
446: int m;
447: struct order_spec *ord;
448: LIST *rp;
449: {
450: int i;
451: struct order_spec ord1;
452: VL fv,vv,vc;
1.5 noro 453: DP b,c,c1;
1.1 noro 454: NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
455:
456: dp_fcoeffs = 0;
457: get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
458: NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
459: CNVars = NVars;
460: if ( ord->id && NVars != ord->nv )
461: error("dp_f4_mod_main : invalid order specification");
462: initd(ord);
463: for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
464: ptod(CO,vv,(P)BDY(t),&b);
465: _dp_mod(b,m,0,&c);
1.5 noro 466: _dp_monic(c,m,&c1);
1.1 noro 467: if ( c ) {
1.5 noro 468: NEXTNODE(fd0,fd); BDY(fd) = (pointer)c1;
1.1 noro 469: }
470: }
471: if ( fd0 ) NEXT(fd) = 0;
472: setup_arrays(fd0,m,&s);
473: x = gb_f4_mod(s,m);
474: reduceall_mod(x,m,&xx); x = xx;
475: for ( r0 = 0; x; x = NEXT(x) ) {
476: NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
477: }
478: if ( r0 ) NEXT(r) = 0;
479: MKLIST(*rp,r0);
480: }
481:
482: NODE gb_f4(f)
483: NODE f;
484: {
485: int i,j,k,nh,row,col,nv;
486: NODE r,g,gall;
1.2 noro 487: NODE s,s0;
1.1 noro 488: DP_pairs d,dm,dr,t;
1.2 noro 489: DP h,nf,nf1,f1,f2,f21,f21r,sp,sp1,sd,sdm,tdp;
1.1 noro 490: MP mp,mp0;
491: NODE blist,bt,nt;
492: DL *ht,*at;
493: MAT mat,nm;
494: int *colstat;
495: int *rind,*cind;
496: int rank,nred;
497: Q dn;
498: struct oEGT tmp0,tmp1,tmp2,eg_split_symb,eg_split_elim;
499: extern struct oEGT eg_mod,eg_elim,eg_chrem,eg_gschk,eg_intrat,eg_symb;
500:
501: init_eg(&eg_mod); init_eg(&eg_elim); init_eg(&eg_chrem);
502: init_eg(&eg_gschk); init_eg(&eg_intrat); init_eg(&eg_symb);
503:
504: doing_f4 = 1;
505: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
506: i = (int)BDY(r);
507: d = updpairs(d,g,i);
508: g = updbase(g,i);
509: gall = append_one(gall,i);
510: }
511: if ( gall )
512: nv = ((DP)ps[(int)BDY(gall)])->nv;
513: while ( d ) {
514: get_eg(&tmp0);
515: minsugar(d,&dm,&dr); d = dr;
1.12 noro 516: if ( DP_Print )
1.1 noro 517: fprintf(asir_out,"sugar=%d\n",dm->sugar);
518: blist = 0; s0 = 0;
519: /* asph : sum of all head terms of spoly */
520: for ( t = dm; t; t = NEXT(t) ) {
521: dp_sp(ps[t->dp1],ps[t->dp2],&sp);
522: if ( sp ) {
523: MKNODE(bt,sp,blist); blist = bt;
1.2 noro 524: s0 = symb_merge(s0,dp_dllist(sp),nv);
1.1 noro 525: }
526: }
527: /* s0 : all the terms appeared in symbolic redunction */
1.2 noro 528: for ( s = s0, nred = 0; s; s = NEXT(s) ) {
1.1 noro 529: for ( r = gall; r; r = NEXT(r) )
1.2 noro 530: if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
1.1 noro 531: break;
532: if ( r ) {
1.2 noro 533: dltod(BDY(s),nv,&tdp);
534: dp_subd(tdp,ps[(int)BDY(r)],&sd);
1.5 noro 535: muld(CO,sd,ps[(int)BDY(r)],&f2);
1.1 noro 536: MKNODE(bt,f2,blist); blist = bt;
1.2 noro 537: s = symb_merge(s,dp_dllist(f2),nv);
1.1 noro 538: nred++;
539: }
540: }
541:
542: /* the first nred polys in blist are reducers */
543: /* row = the number of all the polys */
544: for ( r = blist, row = 0; r; r = NEXT(r), row++ );
545: ht = (DL *)MALLOC(nred*sizeof(DL));
546: for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
547: ht[i] = BDY((DP)BDY(r))->dl;
1.2 noro 548: for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.1 noro 549: at = (DL *)MALLOC(col*sizeof(DL));
1.2 noro 550: for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
551: at[i] = (DL)BDY(s);
1.1 noro 552: MKMAT(mat,row,col);
553: for ( i = 0, r = blist; i < row; r = NEXT(r), i++ )
554: dp_to_vect(BDY(r),at,(Q *)mat->body[i]);
555: get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
556: init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
1.12 noro 557: if ( DP_Print ) {
1.1 noro 558: print_eg("Symb",&eg_split_symb);
559: fprintf(asir_out,"mat : %d x %d",row,col);
560: fflush(asir_out);
561: }
1.3 noro 562: #if 0
563: rank = generic_gauss_elim_hensel(mat,&nm,&dn,&rind,&cind);
564: #else
1.1 noro 565: rank = generic_gauss_elim(mat,&nm,&dn,&rind,&cind);
1.3 noro 566: #endif
1.12 noro 567: if ( DP_Print )
1.1 noro 568: fprintf(asir_out,"done rank = %d\n",rank,row,col);
569: for ( i = 0; i < rank; i++ ) {
570: for ( k = 0; k < nred; k++ )
571: if ( !cmpdl(nv,at[rind[i]],ht[k]) )
572: break;
573: if ( k == nred ) {
574: /* this is a new base */
575: mp0 = 0;
576: NEXTMP(mp0,mp); mp->dl = at[rind[i]]; mp->c = (P)dn;
577: for ( k = 0; k < col-rank; k++ )
578: if ( nm->body[i][k] ) {
579: NEXTMP(mp0,mp); mp->dl = at[cind[k]];
580: mp->c = (P)nm->body[i][k];
581: }
582: NEXT(mp) = 0;
583: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
584: dp_ptozp(nf,&nf1);
585: nh = newps(nf1,0,0);
586: d = updpairs(d,g,nh);
587: g = updbase(g,nh);
588: gall = append_one(gall,nh);
589: }
590: }
591: }
1.12 noro 592: if ( DP_Print ) {
1.1 noro 593: print_eg("Symb",&eg_symb);
594: print_eg("Mod",&eg_mod); print_eg("GaussElim",&eg_elim);
595: print_eg("ChRem",&eg_chrem); print_eg("IntToRat",&eg_intrat);
596: print_eg("Check",&eg_gschk);
597: }
598: return g;
599: }
600:
1.5 noro 601: /* initial bases are monic */
602:
1.1 noro 603: NODE gb_f4_mod(f,m)
604: NODE f;
605: int m;
606: {
607: int i,j,k,nh,row,col,nv;
608: NODE r,g,gall;
1.2 noro 609: NODE s,s0;
1.1 noro 610: DP_pairs d,dm,dr,t;
1.2 noro 611: DP h,nf,f1,f2,f21,f21r,sp,sp1,sd,sdm,tdp;
1.1 noro 612: MP mp,mp0;
613: NODE blist,bt,nt;
1.5 noro 614: DL *ht,*at,*st;
615: int **spmat,**redmat;
616: int *colstat,*w;
617: int rank,nred,nsp,nonzero,spcol;
618: int *indred,*isred,*ri;
619: struct oEGT tmp0,tmp1,tmp2,eg_split_symb,eg_split_elim1,eg_split_elim2;
620: extern struct oEGT eg_symb,eg_elim1,eg_elim2;
1.1 noro 621:
1.5 noro 622: init_eg(&eg_symb); init_eg(&eg_elim1); init_eg(&eg_elim2);
1.1 noro 623: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
624: i = (int)BDY(r);
625: d = updpairs(d,g,i);
626: g = updbase(g,i);
627: gall = append_one(gall,i);
628: }
629: if ( gall )
630: nv = ((DP)ps[(int)BDY(gall)])->nv;
631: while ( d ) {
632: get_eg(&tmp0);
633: minsugar(d,&dm,&dr); d = dr;
1.12 noro 634: if ( DP_Print )
1.1 noro 635: fprintf(asir_out,"sugar=%d\n",dm->sugar);
636: blist = 0; s0 = 0;
637: /* asph : sum of all head terms of spoly */
638: for ( t = dm; t; t = NEXT(t) ) {
1.13 noro 639: _dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp);
1.1 noro 640: if ( sp ) {
641: MKNODE(bt,sp,blist); blist = bt;
1.2 noro 642: s0 = symb_merge(s0,dp_dllist(sp),nv);
1.1 noro 643: }
644: }
645: /* s0 : all the terms appeared in symbolic redunction */
1.2 noro 646: for ( s = s0, nred = 0; s; s = NEXT(s) ) {
1.1 noro 647: for ( r = gall; r; r = NEXT(r) )
1.2 noro 648: if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
1.1 noro 649: break;
650: if ( r ) {
1.2 noro 651: dltod(BDY(s),nv,&tdp);
652: dp_subd(tdp,ps[(int)BDY(r)],&sd);
1.1 noro 653: _dp_mod(sd,m,0,&sdm);
1.13 noro 654: mulmd_dup(m,sdm,ps[(int)BDY(r)],&f2);
1.1 noro 655: MKNODE(bt,f2,blist); blist = bt;
1.2 noro 656: s = symb_merge(s,dp_dllist(f2),nv);
1.1 noro 657: nred++;
658: }
659: }
1.5 noro 660:
661: get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
662: init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
1.1 noro 663:
664: /* the first nred polys in blist are reducers */
665: /* row = the number of all the polys */
666: for ( r = blist, row = 0; r; r = NEXT(r), row++ );
1.5 noro 667:
668: /* head terms of reducers */
1.1 noro 669: ht = (DL *)MALLOC(nred*sizeof(DL));
670: for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
671: ht[i] = BDY((DP)BDY(r))->dl;
1.5 noro 672:
673: /* col = number of all terms */
1.2 noro 674: for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.5 noro 675:
676: /* head terms of all terms */
1.1 noro 677: at = (DL *)MALLOC(col*sizeof(DL));
1.2 noro 678: for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
679: at[i] = (DL)BDY(s);
1.5 noro 680:
681: /* store coefficients separately in spmat and redmat */
682: nsp = row-nred;
683:
684: /* reducer matrix */
685: redmat = (int **)almat(nred,col);
686: for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ )
687: _dpmod_to_vect(BDY(r),at,redmat[i]);
688: /* XXX */
689: /* reduce_reducers_mod(redmat,nred,col,m); */
690: /* register the position of the head term */
691: indred = (int *)MALLOC(nred*sizeof(int));
692: bzero(indred,nred*sizeof(int));
693: isred = (int *)MALLOC(col*sizeof(int));
694: bzero(isred,col*sizeof(int));
695: for ( i = 0; i < nred; i++ ) {
696: ri = redmat[i];
697: for ( j = 0; j < col && !ri[j]; j++ );
698: indred[i] = j;
699: isred[j] = 1;
700: }
701:
702: spcol = col-nred;
703: /* head terms not in ht */
704: st = (DL *)MALLOC(spcol*sizeof(DL));
705: for ( j = 0, k = 0; j < col; j++ )
706: if ( !isred[j] )
707: st[k++] = at[j];
708:
709: /* spoly matrix; stored in reduced form; terms in ht[] are omitted */
710: spmat = almat(nsp,spcol);
711: w = (int *)MALLOC(col*sizeof(int));
712: for ( ; i < row; r = NEXT(r), i++ ) {
713: bzero(w,col*sizeof(int));
714: _dpmod_to_vect(BDY(r),at,w);
715: reduce_sp_by_red_mod(w,redmat,indred,nred,col,m);
716: for ( j = 0, k = 0; j < col; j++ )
717: if ( !isred[j] )
718: spmat[i-nred][k++] = w[j];
719: }
720:
721: get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);
722: init_eg(&eg_split_elim1); add_eg(&eg_split_elim1,&tmp1,&tmp0);
723:
724: colstat = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
725: for ( i = 0, nonzero=0; i < nsp; i++ )
726: for ( j = 0; j < spcol; j++ )
727: if ( spmat[i][j] )
1.1 noro 728: nonzero++;
1.12 noro 729: if ( DP_Print )
1.5 noro 730: fprintf(asir_out,"spmat : %d x %d (nonzero=%f%%)...",
731: nsp,spcol,((double)nonzero*100)/(nsp*spcol));
732: rank = generic_gauss_elim_mod(spmat,nsp,spcol,m,colstat);
733:
734: get_eg(&tmp1); add_eg(&eg_elim2,&tmp0,&tmp1);
735: init_eg(&eg_split_elim2); add_eg(&eg_split_elim2,&tmp0,&tmp1);
736:
1.12 noro 737: if ( DP_Print ) {
1.1 noro 738: fprintf(asir_out,"done rank = %d\n",rank,row,col);
739: print_eg("Symb",&eg_split_symb);
1.5 noro 740: print_eg("Elim1",&eg_split_elim1);
741: print_eg("Elim2",&eg_split_elim2);
1.1 noro 742: fprintf(asir_out,"\n");
743: }
1.5 noro 744: for ( j = 0, i = 0; j < spcol; j++ )
1.1 noro 745: if ( colstat[j] ) {
1.5 noro 746: mp0 = 0;
1.13 noro 747: NEXTMP(mp0,mp); mp->dl = st[j]; mp->c = STOI(1);
1.5 noro 748: for ( k = j+1; k < spcol; k++ )
749: if ( !colstat[k] && spmat[i][k] ) {
1.13 noro 750: NEXTMP(mp0,mp); mp->dl = st[k];
1.5 noro 751: mp->c = STOI(spmat[i][k]);
1.1 noro 752: }
1.5 noro 753: NEXT(mp) = 0;
754: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
755: nh = newps_mod(nf,m);
756: d = updpairs(d,g,nh);
757: g = updbase(g,nh);
758: gall = append_one(gall,nh);
1.1 noro 759: i++;
760: }
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.14 noro 1063: _dp_nf(top,ps[w[i]],ps,1,PtozpRA?DP_Multiple:0,&g);
1.1 noro 1064: prim_part(g,0,&g1);
1065: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12 noro 1066: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1067: fprintf(asir_out,"."); fflush(asir_out);
1068: }
1069: w[i] = newps(g1,0,(NODE)0);
1070: }
1071: for ( top = 0, j = n-1; j >= 0; j-- ) {
1072: MKNODE(t,(pointer)w[j],top); top = t;
1073: }
1074: *h = top;
1.12 noro 1075: if ( DP_Print || DP_PrintShort )
1.1 noro 1076: fprintf(asir_out,"\n");
1077: }
1078:
1079: void reduceall_mod(in,m,h)
1080: NODE in;
1081: int m;
1082: NODE *h;
1083: {
1084: NODE r,t,top;
1085: int n,i,j;
1086: int *w;
1.12 noro 1087: DP g,p;
1.1 noro 1088: struct oEGT tmp0,tmp1;
1089:
1090: if ( NoRA ) {
1091: *h = in; return;
1092: }
1.12 noro 1093: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1094: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1095: }
1096: r = NODE_sortbi(in,0);
1097: n = length(r);
1098: w = (int *)ALLOCA(n*sizeof(int));
1099: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1100: w[i] = (int)BDY(t);
1101: for ( i = 0; i < n; i++ ) {
1102: for ( top = 0, j = n-1; j >= 0; j-- )
1103: if ( j != i ) {
1104: MKNODE(t,(pointer)w[j],top); top = t;
1105: }
1106: get_eg(&tmp0);
1107: if ( PCoeffs )
1108: dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1.12 noro 1109: else {
1110: dpto_dp(ps[w[i]],&p);
1111: _dp_nf_mod_destructive(top,p,ps,m,1,&g);
1112: }
1.1 noro 1113: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12 noro 1114: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1115: fprintf(asir_out,"."); fflush(asir_out);
1116: }
1117: w[i] = newps_mod(g,m);
1118: }
1119: for ( top = 0, j = n-1; j >= 0; j-- ) {
1120: MKNODE(t,(pointer)w[j],top); top = t;
1121: }
1122: *h = top;
1.12 noro 1123: if ( DP_Print || DP_PrintShort )
1.1 noro 1124: fprintf(asir_out,"\n");
1125: }
1126:
1127: int newps(a,m,subst)
1128: DP a;
1129: int m;
1130: NODE subst;
1131: {
1132: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1133: return -1;
1134: if ( psn == pslen ) {
1135: pslen *= 2;
1136: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1137: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1138: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1139: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1140: if ( m )
1141: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1142: }
1143: if ( Demand ) {
1144: if ( doing_f4 )
1145: ps[psn] = a;
1146: else
1147: ps[psn] = 0;
1.7 noro 1148: dp_save(psn,(Obj)a,0);
1.1 noro 1149: } else
1150: ps[psn] = a;
1151: psh[psn] = BDY(a)->dl;
1152: pss[psn] = a->sugar;
1153: psc[psn] = BDY(a)->c;
1154: if ( m )
1155: _dp_mod(a,m,subst,&psm[psn]);
1.7 noro 1156: if ( GenTrace ) {
1157: NODE tn,tr,tr1;
1158: LIST trace;
1159:
1160: /* reverse the TraceList */
1161: tn = TraceList;
1162: for ( tr = 0; tn; tn = NEXT(tn) ) {
1163: MKNODE(tr1,BDY(tn),tr); tr = tr1;
1164: }
1165: MKLIST(trace,tr);
1166: if ( OXCheck >= 0 ) {
1167: NODE arg;
1168: Q q1,q2;
1169: STRING fname;
1170: Obj dmy;
1171:
1172: STOQ(OXCheck,q1);
1173: MKSTR(fname,"check_trace");
1174: STOQ(psn,q2);
1175: arg = mknode(5,q1,fname,a,q2,trace);
1176: Pox_cmo_rpc(arg,&dmy);
1177: } else
1178: dp_save(psn,(Obj)trace,"t");
1179: TraceList = 0;
1180: }
1.1 noro 1181: return psn++;
1182: }
1183:
1184: int newps_nosave(a,m,subst)
1185: DP a;
1186: int m;
1187: NODE subst;
1188: {
1189: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1190: return -1;
1191: if ( psn == pslen ) {
1192: pslen *= 2;
1193: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1194: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1195: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1196: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1197: if ( m )
1198: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1199: }
1200: ps[psn] = 0;
1201: psh[psn] = BDY(a)->dl;
1202: pss[psn] = a->sugar;
1203: psc[psn] = BDY(a)->c;
1204: if ( m )
1205: _dp_mod(a,m,subst,&psm[psn]);
1206: return psn++;
1207: }
1208:
1209: int newps_mod(a,m)
1210: DP a;
1211: int m;
1212: {
1213: if ( psn == pslen ) {
1214: pslen *= 2;
1215: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1216: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1217: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1218: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P)); /* XXX */
1219: }
1220: ps[psn] = a;
1221: psh[psn] = BDY(ps[psn])->dl;
1222: pss[psn] = ps[psn]->sugar;
1223: return psn++;
1224: }
1225:
1226: void reducebase_dehomo(f,g)
1227: NODE f,*g;
1228: {
1229: int n,i,j,k;
1230: int *r;
1231: DL *w,d;
1232: DP u;
1233: NODE t,top;
1234:
1235: n = length(f);
1236: w = (DL *)ALLOCA(n*sizeof(DL));
1237: r = (int *)ALLOCA(n*sizeof(int));
1238: for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {
1239: r[i] = (int)BDY(t); w[i] = psh[r[i]];
1240: }
1241: for ( i = 0; i < n; i++ ) {
1242: for ( j = 0, d = w[i]; j < n; j++ ) {
1243: if ( j != i ) {
1244: for ( k = 0; k < NVars; k++ )
1245: if ( d->d[k] < w[j]->d[k] )
1246: break;
1247: if ( k == NVars )
1248: break;
1249: }
1250: }
1251: if ( j != n )
1252: r[i] = -1;
1253: }
1254: for ( top = 0, i = n-1; i >= 0; i-- )
1255: if ( r[i] >= 0 ) {
1.7 noro 1256: dp_load(r[i],&ps[r[i]]); dp_dehomo(ps[r[i]],&u);
1257: if ( GenTrace ) {
1258: Q q;
1259: LIST hist;
1260: NODE node;
1261:
1262: STOQ(r[i],q);
1263: node = mknode(4,0,q,0,0);
1264: MKLIST(hist,node);
1265: MKNODE(TraceList,hist,0);
1266: }
1267: j = newps(u,0,0);
1.1 noro 1268: MKNODE(t,j,top); top = t;
1269: }
1270: *g = top;
1271: }
1272:
1273: NODE append_one(f,n)
1274: NODE f;
1275: int n;
1276: {
1277: NODE t;
1278:
1279: if ( Reverse || !f ) {
1280: MKNODE(t,(pointer)n,f); return t;
1281: } else {
1282: for ( t = f; NEXT(t); t = NEXT(t) );
1283: MKNODE(NEXT(t),(pointer)n,0);
1284: return f;
1285: }
1286: }
1287:
1288: DP_pairs minp( d, prest )
1289: DP_pairs d, *prest;
1290: {
1291: register DP_pairs m, ml, p, l;
1292: register DL lcm;
1293: register int s, nv = CNVars;
1294: register int (*cmpfun)() = cmpdl;
1295:
1296: if ( !(p = NEXT(m = d)) ) {
1297: *prest = p;
1298: NEXT(m) = 0;
1299: return m;
1300: }
1301: for ( lcm = m->lcm, s = m->sugar, ml = 0, l = m; p; p = NEXT(l = p) )
1302: if ( NoSugar ? (*cmpfun)( nv, lcm, p->lcm ) >= 0 :
1303: (s > p->sugar || s == p->sugar && (*cmpfun)( nv, lcm, p->lcm ) >= 0) )
1304: ml = l, lcm = (m = p)->lcm, s = p->sugar;
1305: if ( !ml ) *prest = NEXT(m);
1306: else {
1307: NEXT(ml) = NEXT(m);
1308: *prest = d;
1309: }
1310: NEXT(m) = 0;
1311: return m;
1312: }
1313:
1314: void minsugar(d,dm,dr)
1315: DP_pairs d;
1316: DP_pairs *dm,*dr;
1317: {
1318: int msugar;
1319: DP_pairs t,dm0,dr0,dmt,drt;
1320:
1321: for ( msugar = d->sugar, t = NEXT(d); t; t = NEXT(t) )
1322: if ( t->sugar < msugar )
1323: msugar = t->sugar;
1324: dm0 = 0; dr0 = 0;
1325: for ( t = d; t; t = NEXT(t) ) {
1326: if ( t->sugar == msugar ) {
1327: NEXTDPP(dm0,dmt);
1328: dmt->dp1 = t->dp1; dmt->dp2 = t->dp2;
1329: dmt->lcm = t->lcm; dmt->sugar = t->sugar;
1330: } else {
1331: NEXTDPP(dr0,drt);
1332: drt->dp1 = t->dp1; drt->dp2 = t->dp2;
1333: drt->lcm = t->lcm; drt->sugar = t->sugar;
1334: }
1335: }
1336: if ( dm0 ) NEXT(dmt) = 0;
1337: if ( dr0 ) NEXT(drt) = 0;
1338: *dm = dm0; *dr = dr0;
1339: }
1340:
1341: NODE gb(f,m,subst)
1342: NODE f;
1343: int m;
1344: NODE subst;
1345: {
1346: int i,nh,prev,mag;
1347: NODE r,g,gall;
1348: DP_pairs d,d1;
1349: DP_pairs l;
1350: DP h,nf,nfm,dp1,dp2;
1351: MP mp;
1352: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tsp0,tsp1,tspm0,tspm1,tnp0,tnp1,tmp0,tmp1;
1353: int skip_nf_flag;
1354: double t_0;
1.7 noro 1355: Q q;
1.6 noro 1356: int new_sugar;
1.1 noro 1357: static prev_sugar = -1;
1358:
1359: Max_mag = 0;
1360: prev = 1;
1361: doing_f4 = 0;
1362: if ( m ) {
1363: psm = (DP *)MALLOC(pslen*sizeof(DP));
1364: for ( i = 0; i < psn; i++ )
1365: if ( psh[i] && !validhc(psc[i],m,subst) )
1366: return 0;
1367: else
1368: _dp_mod(ps[i],m,subst,&psm[i]);
1369: }
1370: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1371: i = (int)BDY(r);
1372: d = updpairs(d,g,i);
1373: g = updbase(g,i);
1374: gall = append_one(gall,i);
1375: }
1376: while ( d ) {
1.14 noro 1377: l = minp(d,&d);
1.1 noro 1378: if ( m ) {
1.6 noro 1379: _dp_sp_mod_dup(psm[l->dp1],psm[l->dp2],m,&h);
1.10 noro 1380: if ( h )
1381: new_sugar = h->sugar;
1.1 noro 1382: get_eg(&tnfm0);
1.6 noro 1383: _dp_nf_mod_destructive(gall,h,psm,m,0,&nfm);
1.1 noro 1384: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1385: } else
1386: nfm = (DP)1;
1387: if ( nfm ) {
1388: if ( Demand ) {
1389: if ( dp_load_t(psn,&nf) ) {
1390: skip_nf_flag = 1;
1.14 noro 1391: tnf1=tsp1;
1.1 noro 1392: goto skip_nf;
1393: } else {
1394: skip_nf_flag = 0;
1395: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1396: dp_sp(dp1,dp2,&h);
1397: }
1398: } else
1399: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1.7 noro 1400: if ( GenTrace ) {
1401: STOQ(l->dp1,q); ARG1(BDY((LIST)BDY(NEXT(TraceList)))) = q;
1402: STOQ(l->dp2,q); ARG1(BDY((LIST)BDY(TraceList))) = q;
1403: }
1.10 noro 1404: if ( h )
1405: new_sugar = h->sugar;
1.1 noro 1406: get_eg(&tnf0);
1407: t_0 = get_rtime();
1.14 noro 1408: _dp_nf(gall,h,ps,!Top,DP_Multiple,&nf);
1.12 noro 1409: if ( DP_Print )
1.1 noro 1410: fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
1411: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1412: } else
1413: nf = 0;
1414: skip_nf:
1415: if ( nf ) {
1416: NZR++;
1417: get_eg(&tpz0);
1418: prim_part(nf,0,&h);
1419: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1420: get_eg(&tnp0);
1421: if ( Demand && skip_nf_flag )
1422: nh = newps_nosave(h,m,subst);
1423: else
1424: nh = newps(h,m,subst);
1425: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1426: if ( nh < 0 )
1427: return 0;
1428: d = updpairs(d,g,nh);
1429: g = updbase(g,nh);
1430: gall = append_one(gall,nh);
1431: if ( !dp_fcoeffs && ShowMag ) {
1432: for ( mag = 0, mp = BDY(h); mp; mp = NEXT(mp) )
1433: mag += p_mag((P)mp->c);
1434: Max_mag = MAX(Max_mag,mag);
1435: }
1.12 noro 1436: if ( DP_Print ) {
1.1 noro 1437: if ( !prev )
1438: fprintf(asir_out,"\n");
1439: print_split_e(&tnf0,&tnf1); print_split_e(&tpz0,&tpz1);
1440: printdl(psh[nh]);
1441: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",
1442: l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
1443: pss[nh]);
1444: if ( ShowMag )
1445: fprintf(asir_out,",mag=%d",mag);
1446: fprintf(asir_out,"\n"); fflush(asir_out);
1.12 noro 1447: } else if ( DP_PrintShort ) {
1.1 noro 1448: fprintf(asir_out,"+"); fflush(asir_out);
1449: }
1450: prev = 1;
1451: } else {
1452: if ( m )
1453: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1454: ZR++;
1.12 noro 1455: if ( DP_Print || DP_PrintShort ) {
1.6 noro 1456: if ( new_sugar != prev_sugar ) {
1457: fprintf(asir_out,"[%d]",new_sugar);
1458: prev_sugar = new_sugar;
1.1 noro 1459: }
1460: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1461: }
1462: }
1463: }
1.12 noro 1464: if ( DP_Print || DP_PrintShort )
1.1 noro 1465: fprintf(asir_out,"gb done\n");
1466: return g;
1467: }
1468:
1469: NODE gb_mod(f,m)
1470: NODE f;
1471: int m;
1472: {
1473: int i,nh,prev;
1474: NODE r,g,gall;
1475: DP_pairs d,d1;
1476: DP_pairs l;
1477: DP h,nf;
1478: struct oEGT tnfm0,tnfm1,tspm0,tspm1,tmp0,tmp1,tpz0,tpz1;
1479:
1480: prev = 1;
1481: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1482: i = (int)BDY(r);
1483: d = updpairs(d,g,i);
1484: g = updbase(g,i);
1485: gall = append_one(gall,i);
1486: }
1487: while ( d ) {
1.14 noro 1488: l = minp(d,&d);
1.1 noro 1489: if ( PCoeffs ) {
1490: dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1.14 noro 1491: get_eg(&tnfm0);
1.1 noro 1492: dp_nf_mod(gall,h,ps,m,!Top,&nf);
1493: } else {
1.6 noro 1494: _dp_sp_mod_dup(ps[l->dp1],ps[l->dp2],m,&h);
1.14 noro 1495: get_eg(&tnfm0);
1.6 noro 1496: _dp_nf_mod_destructive(gall,h,ps,m,!Top,&nf);
1.1 noro 1497: }
1498: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1499: if ( nf ) {
1500: NZR++;
1501: get_eg(&tpz0);
1502: prim_part(nf,m,&h);
1503: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1504: nh = newps_mod(h,m);
1505: if ( nh < 0 )
1506: return 0;
1507: d = updpairs(d,g,nh);
1508: g = updbase(g,nh);
1509: gall = append_one(gall,nh);
1.12 noro 1510: if ( DP_Print ) {
1.1 noro 1511: if ( !prev )
1512: fprintf(asir_out,"\n");
1513: print_split_eg(&tnfm0,&tnfm1); fflush(asir_out);
1514: 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]);
1515: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1.12 noro 1516: } else if ( DP_PrintShort ) {
1.6 noro 1517: fprintf(asir_out,"+"); fflush(asir_out);
1.1 noro 1518: }
1519: prev = 1;
1520: } else {
1521: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1522: ZR++;
1.12 noro 1523: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1524: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1525: }
1526: }
1527: }
1.12 noro 1528: if ( DP_Print || DP_PrintShort )
1.1 noro 1529: fprintf(asir_out,"gb_mod done\n");
1530: return g;
1531: }
1532:
1.14 noro 1533: DP_pairs updpairs( d, g, t)
1534: DP_pairs d;
1535: NODE /* of index */ g;
1536: int t;
1.1 noro 1537: {
1.14 noro 1538: register DP_pairs d1, dd, nd;
1539: struct oEGT tup0,tup1;
1540: int dl,dl1;
1.1 noro 1541:
1.14 noro 1542: if ( !g ) return d;
1543: if ( !NoCriB && d ) {
1544: dl = DPPlength(d);
1545: d = criterion_B( d, t );
1546: dl -= DPPlength(d); NBP += dl;
1.1 noro 1547: }
1.14 noro 1548: d1 = newpairs( g, t );
1549: if ( NEXT(d1) ) {
1550: dl = DPPlength(d1); TP += dl;
1551: d1 = criterion_M( d1 );
1552: dl1 = DPPlength(d1); NMP += (dl-dl1); dl = dl1;
1553: d1 = criterion_F( d1 );
1554: dl1 = DPPlength(d1); NFP += (dl-dl1); dl = dl1;
1555: } else
1556: dl = 1;
1557: if ( !do_weyl )
1558: for ( dd = 0; d1; d1 = nd ) {
1559: nd = NEXT(d1);
1560: if ( !criterion_2( d1->dp1, d1->dp2 ) ) {
1561: NEXT(d1) = dd;
1562: dd = d1;
1.1 noro 1563: }
1564: }
1.14 noro 1565: else
1566: dd = d1;
1567: dl1 = DPPlength(dd); NDP += (dl-dl1);
1568: if ( !(nd = d) ) return dd;
1569: while ( nd = NEXT(d1 = nd) ) ;
1570: NEXT(d1) = dd;
1571: return d;
1.1 noro 1572: }
1573:
1.14 noro 1574: DP_pairs newpairs( g, t )
1575: NODE /* of index */ g;
1576: register int t;
1577: {
1578: register NODE r;
1579: register DL tdl = psh[t];
1580: register int ts;
1581: register DP_pairs p, last;
1582: int dp;
1583: register DL dl;
1584: register int s;
1.1 noro 1585:
1.14 noro 1586: ts = pss[t] - tdl->td;
1587: for ( last = 0, r = g; r; r = NEXT(r) ) {
1588: NEXT(p = NEWDP_pairs) = last;
1589: last = p;
1590: dp = p->dp1 = (int)BDY(r); p->dp2 = t;
1591: p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
1592: #if 0
1593: if ( do_weyl )
1594: p->sugar = dl_weight(p->lcm);
1595: else
1596: #endif
1597: p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
1598: }
1599: return last;
1600: }
1.1 noro 1601:
1.14 noro 1602: DP_pairs criterion_B( d, s )
1603: register DP_pairs d;
1604: int s;
1.1 noro 1605: {
1.14 noro 1606: register DP_pairs dd, p;
1607: register DL tij, t = psh[s], dltmp;
1.1 noro 1608:
1.14 noro 1609: if ( !d ) return 0;
1610: NEWDL( dltmp, CNVars );
1611: for ( dd = 0; d; d = p ) {
1612: p = NEXT(d),
1613: tij = d->lcm;
1614: if ( tij->td != lcm_of_DL(CNVars, tij, t, dltmp )->td
1615: || !dl_equal(CNVars, tij, dltmp )
1616: || (tij->td == lcm_of_DL(CNVars, psh[d->dp1], t, dltmp )->td
1617: && dl_equal(CNVars, dltmp, tij ))
1618: || (tij->td == lcm_of_DL(CNVars, psh[d->dp2], t, dltmp )->td
1619: && dl_equal(CNVars, dltmp, tij )) ) {
1620: NEXT(d) = dd;
1621: dd = d;
1622: }
1.1 noro 1623: }
1.14 noro 1624: return dd;
1625: }
1.1 noro 1626:
1.14 noro 1627: DP_pairs criterion_M( d1 )
1628: DP_pairs d1;
1629: {
1630: register DP_pairs dd, e, d3, d2, p;
1631: register DL itdl, jtdl;
1632: register int itdltd, jtdltd;
1.1 noro 1633:
1.14 noro 1634: for ( dd = 0, e = d1; e; e = d3 ) {
1635: if ( !(d2 = NEXT(e)) ) {
1636: NEXT(e) = dd;
1637: return e;
1638: }
1639: itdltd = (itdl = e->lcm)->td;
1640: for ( d3 = 0; d2; d2 = p ) {
1641: p = NEXT(d2),
1642: jtdltd = (jtdl = d2->lcm)->td;
1643: if ( jtdltd == itdltd )
1644: if ( dl_equal(CNVars, itdl, jtdl ) ) ;
1645: else if ( dl_redble( jtdl, itdl ) ) continue;
1646: else if ( dl_redble( itdl, jtdl ) ) goto delit;
1647: else ;
1648: else if ( jtdltd > itdltd )
1649: if ( dl_redble( jtdl, itdl ) ) continue;
1.1 noro 1650: else ;
1651: else if ( dl_redble( itdl, jtdl ) ) goto delit;
1652: NEXT(d2) = d3;
1653: d3 = d2;
1654: }
1655: NEXT(e) = dd;
1656: dd = e;
1657: continue;
1658: /**/
1659: delit: NEXT(d2) = d3;
1660: d3 = d2;
1661: for ( ; p; p = d2 ) {
1662: d2 = NEXT(p);
1663: NEXT(p) = d3;
1664: d3 = p;
1665: }
1666: }
1667: return dd;
1668: }
1669:
1670: static DP_pairs collect_pairs_of_hdlcm( d1, prest )
1671: DP_pairs d1, *prest;
1672: {
1673: register DP_pairs w, p, r, s;
1674: register DL ti;
1675: register int td;
1676:
1677: td = (ti = (w = d1)->lcm)->td;
1678: s = NEXT(w);
1679: NEXT(w) = 0;
1680: for ( r = 0; s; s = p ) {
1681: p = NEXT(s);
1682: if ( td == s->lcm->td && dl_equal(CNVars, ti, s->lcm ) )
1683: {
1684: NEXT(s) = w;
1685: w = s;
1686: } else {
1687: NEXT(s) = r;
1688: r = s;
1689: }
1690: }
1691: *prest = r;
1692: return w;
1693: }
1694:
1695: int criterion_2( dp1, dp2 )
1696: int dp1, dp2;
1697: {
1698: register int i, *d1, *d2;
1699:
1700: d1 = psh[dp1]->d, d2 = psh[dp2]->d;
1701: for ( i = CNVars; --i >= 0; d1++, d2++ )
1702: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
1703: return 1;
1704: }
1705:
1706: DP_pairs criterion_F( d1 )
1707: DP_pairs d1;
1708: {
1709: DP_pairs rest, head;
1710: register DP_pairs last, p, r, w;
1711: register int s;
1712:
1713: for ( head = last = 0, p = d1; NEXT(p); ) {
1714: s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
1715: while ( w = NEXT(w) )
1716: if ( criterion_2( w->dp1, w->dp2 ) ) {
1717: r = w;
1718: break;
1719: } else if ( w->sugar < s ) s = (r = w)->sugar;
1720: if ( last ) NEXT(last) = r;
1721: else head = r;
1722: NEXT(last = r) = 0;
1723: if ( !(p = rest) ) return head;
1724: }
1725: if ( !last ) return p;
1726: NEXT(last) = p;
1727: return head;
1728: }
1729:
1730: NODE updbase(g,t)
1731: NODE g;
1732: int t;
1733: {
1734: g = remove_reducibles(g,t);
1735: g = append_one(g,t);
1736: return g;
1737: }
1738:
1739: NODE /* of index */ remove_reducibles( nd, newdp )
1740: NODE /* of index */ nd;
1741: int newdp;
1742: {
1743: register DL dl, dln;
1744: register NODE last, p, head;
1745: register int td;
1746:
1747: dl = psh[newdp];
1748: td = dl->td;
1749: for ( head = last = 0, p = nd; p; ) {
1750: dln = psh[(int)BDY(p)];
1751: if ( dln->td >= td && dl_redble( dln, dl ) ) {
1752: p = NEXT(p);
1753: if ( last ) NEXT(last) = p;
1754: } else {
1755: if ( !last ) head = p;
1756: p = NEXT(last = p);
1757: }
1758: }
1759: return head;
1760: }
1761:
1762: int dl_redble( dl1, dl2 )
1763: DL dl1, dl2;
1764: {
1765: register int n, *d1, *d2;
1766:
1767: for ( d1 = dl1->d, d2 = dl2->d, n = CNVars; --n >= 0; d1++, d2++ )
1768: if ( *d1 < *d2 ) return 0;
1769: return 1;
1.5 noro 1770: }
1771:
1772: int dl_weight(dl)
1773: DL dl;
1774: {
1775: int n,w,i;
1776:
1777: n = CNVars/2;
1778: for ( i = 0, w = 0; i < n; i++ )
1779: w += (-dl->d[i]+dl->d[n+i]);
1780: return w;
1.1 noro 1781: }
1782:
1783: int gbcheck(f)
1784: NODE f;
1785: {
1786: int i;
1787: NODE r,g,gall;
1788: DP_pairs d,l;
1789: DP h,nf,dp1,dp2;
1790: struct oEGT tmp0,tmp1;
1791:
1792: if ( NoGC )
1793: return 1;
1794: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1795: i = (int)BDY(r);
1796: d = updpairs(d,g,i);
1797: g = updbase(g,i);
1798: gall = append_one(gall,i);
1799: }
1.12 noro 1800: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1801: fprintf(asir_out,"gbcheck total %d pairs\n",DPPlength(d)); fflush(asir_out);
1802: }
1803: while ( d ) {
1804: l = d; d = NEXT(d);
1805: get_eg(&tmp0);
1806: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2); dp_sp(dp1,dp2,&h);
1.14 noro 1807: _dp_nf(gall,h,ps,1,0,&nf);
1.1 noro 1808: get_eg(&tmp1); add_eg(&eg_gc,&tmp0,&tmp1);
1.12 noro 1809: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1810: fprintf(asir_out,"."); fflush(asir_out);
1811: }
1812: if ( nf )
1813: return 0;
1814: }
1.12 noro 1815: if ( DP_Print || DP_PrintShort )
1.1 noro 1816: fprintf(asir_out,"\n");
1817: return 1;
1818: }
1819:
1820: int membercheck(f,x)
1821: NODE f,x;
1822: {
1823: DP g;
1824: struct oEGT tmp0,tmp1;
1825:
1826: if ( NoMC )
1827: return 1;
1.12 noro 1828: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1829: fprintf(asir_out,"membercheck\n"); fflush(asir_out);
1830: }
1831: for ( ; f; f = NEXT(f) ) {
1832: get_eg(&tmp0);
1.14 noro 1833: _dp_nf(x,(DP)BDY(f),ps,1,0,&g);
1.1 noro 1834: get_eg(&tmp1); add_eg(&eg_mc,&tmp0,&tmp1);
1.12 noro 1835: if ( DP_Print ) {
1.1 noro 1836: print_split_eg(&tmp0,&tmp1); fflush(asir_out);
1.12 noro 1837: } else if ( DP_PrintShort ) {
1.1 noro 1838: fprintf(asir_out,"."); fflush(asir_out);
1839: }
1840: if ( g )
1841: return 0;
1842: }
1.12 noro 1843: if ( DP_Print || DP_PrintShort )
1.1 noro 1844: fprintf(asir_out,"\n");
1845: return 1;
1846: }
1847:
1848: void dp_set_flag(name,value)
1849: Obj name,value;
1850: {
1851: char *n;
1852: int v;
1853:
1854: if ( OID(name) != O_STR )
1855: return;
1856: n = BDY((STRING)name);
1857: if ( !strcmp(n,"Demand") ) {
1858: Demand = value ? BDY((STRING)value) : 0; return;
1859: }
1860: if ( !strcmp(n,"Dist") ) {
1861: Dist = (LIST)value; return;
1862: }
1863: if ( value && OID(value) != O_N )
1864: return;
1865: v = QTOS((Q)value);
1866: if ( !strcmp(n,"NoSugar") )
1867: NoSugar = v;
1868: else if ( !strcmp(n,"NoCriB") )
1869: NoCriB = v;
1870: else if ( !strcmp(n,"NoGC") )
1871: NoGC = v;
1872: else if ( !strcmp(n,"NoMC") )
1873: NoMC = v;
1874: else if ( !strcmp(n,"NoRA") )
1875: NoRA = v;
1876: else if ( !strcmp(n,"NoGCD") )
1877: NoGCD = v;
1878: else if ( !strcmp(n,"Top") )
1879: Top = v;
1880: else if ( !strcmp(n,"ShowMag") )
1881: ShowMag = v;
1.14 noro 1882: else if ( !strcmp(n,"PrintShort") )
1.12 noro 1883: DP_PrintShort = v;
1.14 noro 1884: else if ( !strcmp(n,"Print") )
1.12 noro 1885: DP_Print = v;
1.14 noro 1886: else if ( !strcmp(n,"NFStat") )
1887: DP_NFStat = v;
1.1 noro 1888: else if ( !strcmp(n,"Stat") )
1889: Stat = v;
1890: else if ( !strcmp(n,"Reverse") )
1891: Reverse = v;
1892: else if ( !strcmp(n,"Multiple") )
1.14 noro 1893: DP_Multiple = v;
1.1 noro 1894: else if ( !strcmp(n,"Denominator") )
1895: Denominator = v;
1896: else if ( !strcmp(n,"PtozpRA") )
1897: PtozpRA = v;
1.7 noro 1898: else if ( !strcmp(n,"GenTrace") )
1899: GenTrace = v;
1900: else if ( !strcmp(n,"OXCheck") )
1901: OXCheck = v;
1.1 noro 1902: }
1903:
1904: void dp_make_flaglist(list)
1905: LIST *list;
1906: {
1907: Q v;
1908: STRING name,path;
1909: NODE n,n1;
1910:
1.14 noro 1911: STOQ(DP_Multiple,v); MKNODE(n,v,0); MKSTR(name,"DP_Multiple"); MKNODE(n1,name,n); n = n1;
1.1 noro 1912: STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
1913: MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
1914: STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
1915: STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
1.14 noro 1916: STOQ(DP_Print,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Print"); MKNODE(n1,name,n); n = n1;
1917: STOQ(DP_NFStat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NFStat"); MKNODE(n1,name,n); n = n1;
1.7 noro 1918: STOQ(OXCheck,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"OXCheck"); MKNODE(n1,name,n); n = n1;
1919: STOQ(GenTrace,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenTrace"); MKNODE(n1,name,n); n = n1;
1.1 noro 1920: STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
1921: STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
1922: STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
1923: STOQ(NoGCD,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGCD"); MKNODE(n1,name,n); n = n1;
1924: STOQ(NoRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoRA"); MKNODE(n1,name,n); n = n1;
1925: STOQ(NoMC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoMC"); MKNODE(n1,name,n); n = n1;
1926: STOQ(NoGC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGC"); MKNODE(n1,name,n); n = n1;
1927: STOQ(NoCriB,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoCriB"); MKNODE(n1,name,n); n = n1;
1928: STOQ(NoSugar,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoSugar"); MKNODE(n1,name,n); n = n1;
1929: if ( Demand )
1930: MKSTR(path,Demand);
1931: else
1932: path = 0;
1933: MKNODE(n1,path,n); n = n1; MKSTR(name,"Demand"); MKNODE(n1,name,n); n = n1;
1934: MKLIST(*list,n);
1935: }
1936:
1937: #define DELIM '/'
1938:
1.7 noro 1939: void dp_save(index,p,prefix)
1.1 noro 1940: int index;
1.7 noro 1941: Obj p;
1942: char *prefix;
1.1 noro 1943: {
1944: FILE *fp;
1945: char path[BUFSIZ];
1946:
1.7 noro 1947: if ( prefix )
1948: sprintf(path,"%s%c%s%d",Demand,DELIM,prefix,index);
1949: else
1950: sprintf(path,"%s%c%d",Demand,DELIM,index);
1.1 noro 1951: if ( !(fp = fopen(path,"wb") ) )
1952: error("dp_save : cannot open a file");
1.7 noro 1953: savevl(fp,VC); saveobj(fp,p); fclose(fp);
1.1 noro 1954: }
1955:
1956: void dp_load(index,p)
1957: int index;
1958: DP *p;
1959: {
1960: FILE *fp;
1961: char path[BUFSIZ];
1962:
1963: if ( !Demand || ps[index] )
1964: *p = ps[index];
1965: else {
1966: sprintf(path,"%s%c%d",Demand,DELIM,index);
1967: if ( !(fp = fopen(path,"rb") ) )
1968: error("dp_load : cannot open a file");
1969: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp);
1970: }
1971: }
1972:
1973: int dp_load_t(index,p)
1974: int index;
1975: DP *p;
1976: {
1977: FILE *fp;
1978: char path[BUFSIZ];
1979:
1980: sprintf(path,"%s%c%d",Demand,DELIM,index);
1981: if ( !(fp = fopen(path,"rb") ) )
1982: return 0;
1983: else {
1984: skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp); return 1;
1985: }
1986: }
1987:
1988: void init_stat() {
1989: init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
1.14 noro 1990: init_eg(&eg_pz); init_eg(&eg_np);
1.1 noro 1991: init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
1.14 noro 1992: ZR = NZR = TP = NBP = NFP = NDP = 0;
1.1 noro 1993: }
1994:
1995: void print_stat() {
1.12 noro 1996: if ( !DP_Print && !Stat )
1.1 noro 1997: return;
1998: print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
1.14 noro 1999: print_eg("PZ",&eg_pz); print_eg("NP",&eg_np);
1.1 noro 2000: print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
2001: 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);
2002: }
1.14 noro 2003:
2004: /*
2005: * dp_nf used in gb()
2006: *
2007: */
2008:
2009: double pz_t_e, pz_t_d, pz_t_d1, pz_t_c, im_t_s, im_t_r;
2010:
2011: extern int GenTrace;
2012: extern NODE TraceList;
1.15 ! noro 2013: extern int mpi_mag;
! 2014:
! 2015: void dp_mulc_d(p,c,r)
! 2016: DP p;
! 2017: P c;
! 2018: DP *r;
! 2019: {
! 2020: if ( Dist && BDY(Dist)
! 2021: && HMAG(p) > mpi_mag
! 2022: && p_mag((P)c) > mpi_mag ) {
! 2023: if ( DP_NFStat ) fprintf(asir_out,"~");
! 2024: dp_imul_d(p,(Q)c,r);
! 2025: } else {
! 2026: if ( DP_NFStat ) fprintf(asir_out,"_");
! 2027: muldc(CO,p,c,r);
! 2028: }
! 2029: }
1.14 noro 2030:
2031: void _dp_nf(b,g,ps,full,multiple,r)
2032: NODE b;
2033: DP g;
2034: DP *ps;
2035: int full,multiple;
2036: DP *r;
2037: {
2038: DP u,dp,rp,t,t1,t2,red,shift;
2039: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
2040: N gn,tn,cn;
2041: NODE l;
2042: MP m,mr;
2043: int hmag,denom;
2044: int sugar,psugar;
2045: NODE dist;
2046: STRING imul;
2047: int ndist;
2048: int kara_bit;
2049: double get_rtime();
1.15 ! noro 2050: double t_0,t_00,tt,ttt,t_p,t_m,t_g,t_a;
1.14 noro 2051: LIST hist;
2052: NODE node;
2053: Q rcred,mrcred;
2054:
2055: if ( !g ) {
2056: *r = 0; return;
2057: }
2058: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1.15 ! noro 2059: t_p = t_m = t_g = t_a = 0;
1.14 noro 2060:
2061: denom = Denominator?Denominator:1;
2062: hmag = multiple*HMAG(g)/denom;
2063: if ( Dist ) {
2064: dist = BDY(Dist);
2065: ndist = length(dist);
2066: }
2067: sugar = g->sugar;
2068:
2069: dc = 0; dp = 0; rc = ONE; rp = g;
2070: MKSTR(imul,"dp_imul_index");
2071:
2072: /* g = dc*dp+rc*rp */
2073: for ( ; rp; ) {
2074: for ( u = 0, l = b; l; l = NEXT(l) ) {
2075: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
2076: t_0 = get_rtime();
2077: dp_load((int)BDY(l),&red);
2078: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
2079: igcd_cofactor((Q)BDY(rp)->c,(Q)BDY(red)->c,&gcd,&cred,&cr);
1.15 ! noro 2080: tt = get_rtime(); t_p += tt-t_0;
1.14 noro 2081:
2082: dp_subd(rp,red,&shift);
1.15 ! noro 2083: dp_mulc_d(rp,cr,&t);
1.14 noro 2084: chsgnp((P)cred,(P *)&mcred);
1.15 ! noro 2085: dp_mulc_d(red,mcred,&t1);
! 2086: muld(CO,shift,t1,&t1);
! 2087: addd(CO,t,t1,&u);
! 2088: t_m += get_rtime()-tt;
1.14 noro 2089:
2090: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
2091: sugar = MAX(sugar,psugar);
1.15 ! noro 2092:
1.14 noro 2093: if ( GenTrace ) {
2094: /* u = cr*rp + (-cred)*shift*red */
2095: STOQ((int)BDY(l),cq);
2096: node = mknode(4,cr,cq,0,0);
2097: mulq(cred,rc,&rcred);
2098: chsgnnum((Num)rcred,(Num *)&mrcred);
2099: muldc(CO,shift,(P)mrcred,(DP *)&ARG2(node));
2100: MKLIST(hist,node);
2101: }
1.15 ! noro 2102:
1.14 noro 2103: if ( !u ) {
2104: if ( dp )
2105: dp->sugar = sugar;
2106: *r = dp;
2107: if ( GenTrace ) {
2108: ARG3(BDY(hist)) = ONE;
2109: MKNODE(node,hist,TraceList); TraceList = node;
2110: }
2111: goto final;
2112: }
2113: break;
2114: }
2115: }
2116: if ( u ) {
2117: if ( multiple && HMAG(u) > hmag ) {
2118: t_0 = get_rtime();
2119: if ( Dist && HMAG(u) > mpi_mag ) {
2120: if ( DP_NFStat )
2121: fprintf(asir_out,"D");
1.15 ! noro 2122: dp_ptozp_d(dist,ndist,u,&rp);
1.14 noro 2123: } else {
2124: if ( DP_NFStat )
2125: fprintf(asir_out,"L");
1.15 ! noro 2126: dp_ptozp_d(0,0,u,&rp);
1.14 noro 2127: }
1.15 ! noro 2128: tt = get_rtime(); t_g += tt-t_0;
! 2129:
! 2130: divsq((Q)BDY(u)->c,(Q)BDY(rp)->c,&cont);
1.14 noro 2131: if ( !dp_fcoeffs && DP_NFStat ) {
1.15 ! noro 2132: fprintf(asir_out,
! 2133: "(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1.14 noro 2134: fflush(asir_out);
2135: }
1.15 ! noro 2136: mulq(cr,dc,&dcq); mulq(cont,rc,&rcq);
1.14 noro 2137: igcd_cofactor(dcq,rcq,&gcd,&dc,&rc);
1.15 ! noro 2138: t_a = get_rtime()-tt;
! 2139:
1.14 noro 2140: hmag = multiple*HMAG(rp)/denom;
2141: if ( GenTrace ) {
2142: ARG3(BDY(hist)) = (pointer)gcd;
2143: MKNODE(node,hist,TraceList); TraceList = node;
2144: }
2145: } else {
1.15 ! noro 2146: rp = u;
1.14 noro 2147: t_0 = get_rtime();
1.15 ! noro 2148: mulq(cr,dc,&dc);
! 2149: t_a += get_rtime()-t_0;
1.14 noro 2150: if ( GenTrace ) {
2151: ARG3(BDY(hist)) = (pointer)ONE;
2152: MKNODE(node,hist,TraceList); TraceList = node;
2153: }
2154: }
2155: } else if ( !full ) {
2156: if ( rp ) {
2157: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
2158: }
2159: *r = rp;
2160: goto final;
2161: } else {
2162: t_0 = get_rtime();
2163: mulq((Q)BDY(rp)->c,rc,&c);
1.15 ! noro 2164: igcd_cofactor(dc,c,&dc,&dcq,&cq);
! 2165: muldc(CO,dp,(P)dcq,&t);
! 2166: dp_hm(rp,&t1); BDY(t1)->c = (P)cq; addd(CO,t,t1,&dp);
! 2167: dp_rest(rp,&rp);
! 2168: t_a += get_rtime()-t_0;
1.14 noro 2169: }
2170: }
2171: if ( GenTrace ) {
2172: mulq(ARG3(BDY((LIST)BDY(TraceList))),dc,&cq);
2173: ARG3(BDY((LIST)BDY(TraceList))) = (pointer)cq;
2174: }
2175: if ( dp )
2176: dp->sugar = sugar;
2177: *r = dp;
2178: final:
2179: if ( DP_NFStat )
1.15 ! noro 2180: fprintf(asir_out,
! 2181: "(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
! 2182: t_p,t_m,t_g,t_a,
1.14 noro 2183: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
2184: }
2185:
2186: void imulv();
2187:
2188: void dp_imul_d(p,q,rp)
2189: DP p;
2190: Q q;
2191: DP *rp;
2192: {
2193: int nsep,ndist,i,j,k,l,n;
2194: double t0,t1,t2;
2195: Q *s;
2196: pointer *b;
2197: VECT c,cs,ri;
2198: VECT *r;
2199: MP m;
2200: NODE tn,dist,n0,n1,n2;
2201: Obj dmy;
2202: STRING imul;
2203:
2204: extern LIST Dist;
2205:
2206: if ( !p || !q ) {
2207: *rp = 0; return;
2208: }
2209: dist = BDY(Dist);
2210: for ( tn = dist, ndist = 0; tn; tn = NEXT(tn), ndist++ );
2211: nsep = ndist + 1;
2212: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
2213: if ( n <= nsep ) {
2214: muldc(CO,p,(P)q,rp); return;
2215: }
2216: MKSTR(imul,"imulv");
2217: t0 = get_rtime();
2218: dp_dtov(p,&c);
2219: sepvect(c,nsep,&cs);
2220: r = (VECT *)CALLOC(nsep,sizeof(VECT *));
2221: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
2222: n0 = mknode(4,BDY(tn),imul,b[i],q);
2223: Pox_rpc(n0,&dmy);
2224: }
2225: t1 = get_rtime();
2226: im_t_s += t1 - t0;
2227: imulv(b[i],q,&r[i]);
2228: t1 = get_rtime();
2229: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
2230: MKNODE(n0,BDY(tn),0);
2231: Pox_pop_local(n0,&r[i]);
2232: if ( OID(r[i]) == O_ERR ) {
2233: printexpr(CO,(Obj)r[i]);
2234: error("dp_imul_d : aborted");
2235: }
2236: }
2237: t2 = get_rtime();
2238: im_t_r += t2 - t1;
2239: s = (Q *)CALLOC(n,sizeof(Q));
2240: for ( i = j = 0; i < nsep; i++ ) {
2241: for ( k = 0, ri = r[i], l = ri->len; k < l; k++, j++ ) {
2242: s[j] = (Q)BDY(ri)[k];
2243: }
2244: }
2245: dp_vtod(s,p,rp);
2246: }
2247:
2248: void imulv(w,c,rp)
2249: VECT w;
2250: Q c;
2251: VECT *rp;
2252: {
2253: int n,i;
2254: VECT r;
2255:
2256: n = w->len;
2257: MKVECT(r,n); *rp = r;
2258: for ( i = 0; i < n; i++ )
2259: mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
2260: }
2261: