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