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