Annotation of OpenXM_contrib2/asir2000/builtin/gr.c, Revision 1.48
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.48 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.47 2003/05/30 06:03:28 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;
1.46 noro 92: int DP_PrintShort = 0;
1.14 noro 93: int DP_Multiple = 0;
1.16 noro 94: int DP_NFStat = 0;
1.14 noro 95: LIST Dist = 0;
96: int NoGCD = 0;
97: int GenTrace = 0;
98: int OXCheck = -1;
99:
1.1 noro 100: static int NoSugar = 0;
101: static int NoCriB = 0;
102: static int NoGC = 0;
103: static int NoMC = 0;
104: static int NoRA = 0;
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.48 ! 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.48 ! noro 671: /* fprintf(stderr,"%d-",length(s0)); */
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 {
1.44 noro 1164: if ( dp_fcoeffs || PCoeffs )
1.1 noro 1165: dp_prim(f,r);
1166: else
1167: dp_ptozp(f,r);
1.7 noro 1168: if ( GenTrace && TraceList ) {
1.23 noro 1169: /* adust the denominator according to the final
1170: content reduction */
1.7 noro 1171: divsp(CO,BDY(f)->c,BDY(*r)->c,&d);
1172: mulp(CO,(P)ARG3(BDY((LIST)BDY(TraceList))),d,&t);
1173: ARG3(BDY((LIST)BDY(TraceList))) = t;
1174: }
1.1 noro 1175: }
1176: }
1177:
1.37 noro 1178: NODE /* of DP */ NODE_sortb_insert( DP newdp, NODE /* of DP */ nd, int dec )
1.1 noro 1179: {
1180: register NODE last, p;
1181: register DL newdl = BDY(newdp)->dl;
1182: register int (*cmpfun)() = cmpdl, nv = CNVars;
1183: NODE newnd;
1184: int sgn = dec ? 1 : -1;
1185: MKNODE( newnd, newdp, 0 );
1186: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(last))->dl ) > 0 ) {
1187: NEXT(newnd) = last;
1188: return newnd;
1189: }
1190: for ( ; p = NEXT(last); last = p )
1191: if ( sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(p))->dl ) > 0 ) break;
1192: if ( p ) NEXT(NEXT(last) = newnd) = p;
1193: else NEXT(last) = newnd;
1194: return nd;
1195: }
1196:
1.37 noro 1197: NODE NODE_sortb( NODE node, int dec )
1.1 noro 1198: {
1199: register NODE nd, ans;
1200:
1201: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
1202: ans = NODE_sortb_insert( (DP) BDY(nd), ans, dec );
1203: return ans;
1204: }
1205:
1.37 noro 1206: NODE /* of index */ NODE_sortbi_insert( int newdpi, NODE /* of index */ nd, int dec )
1.1 noro 1207: {
1208: register NODE last, p;
1209: register DL newdl = psh[newdpi];
1210: register int (*cmpfun)() = cmpdl, nv = CNVars;
1211: NODE newnd;
1212: int sgn = dec ? 1 : -1;
1213: MKNODE( newnd, newdpi, 0 );
1214: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(last)] ) > 0 ) {
1215: NEXT(newnd) = last;
1216: return newnd;
1217: }
1218: for ( ; p = NEXT(last); last = p )
1219: if ( sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(p)] ) > 0 ) break;
1220: if ( p ) NEXT(NEXT(last) = newnd) = p;
1221: else NEXT(last) = newnd;
1222: return nd;
1223: }
1224:
1.37 noro 1225: NODE NODE_sortbi( NODE node, int dec )
1.1 noro 1226: {
1227: register NODE nd, ans;
1228:
1229: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
1230: ans = NODE_sortbi_insert( (int) BDY(nd), ans, dec );
1231: return ans;
1232: }
1233:
1.37 noro 1234: void reduceall(NODE in,NODE *h)
1.1 noro 1235: {
1236: NODE r,t,top;
1237: int n,i,j;
1238: int *w;
1239: DP g,g1;
1240: struct oEGT tmp0,tmp1;
1241:
1242: if ( NoRA ) {
1243: *h = in; return;
1244: }
1.12 noro 1245: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1246: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1247: }
1248: r = NODE_sortbi(in,0);
1249: n = length(r);
1250: w = (int *)ALLOCA(n*sizeof(int));
1251: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1252: w[i] = (int)BDY(t);
1253: for ( i = 0; i < n; i++ ) {
1254: for ( top = 0, j = n-1; j >= 0; j-- )
1255: if ( j != i ) {
1256: MKNODE(t,(pointer)w[j],top); top = t;
1257: }
1258: get_eg(&tmp0);
1259: dp_load(w[i],&ps[w[i]]);
1260:
1.7 noro 1261: if ( GenTrace ) {
1262: Q q;
1263: NODE node;
1264: LIST hist;
1265:
1266: STOQ(w[i],q);
1267: node = mknode(4,ONE,q,ONE,ONE);
1268: MKLIST(hist,node);
1269: MKNODE(TraceList,hist,0);
1270: }
1.16 noro 1271: _dp_nf(top,ps[w[i]],ps,1,&g);
1.1 noro 1272: prim_part(g,0,&g1);
1273: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12 noro 1274: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1275: fprintf(asir_out,"."); fflush(asir_out);
1276: }
1277: w[i] = newps(g1,0,(NODE)0);
1278: }
1279: for ( top = 0, j = n-1; j >= 0; j-- ) {
1280: MKNODE(t,(pointer)w[j],top); top = t;
1281: }
1282: *h = top;
1.12 noro 1283: if ( DP_Print || DP_PrintShort )
1.1 noro 1284: fprintf(asir_out,"\n");
1285: }
1286:
1.37 noro 1287: void reduceall_mod(NODE in,int m,NODE *h)
1.1 noro 1288: {
1289: NODE r,t,top;
1290: int n,i,j;
1291: int *w;
1.12 noro 1292: DP g,p;
1.1 noro 1293: struct oEGT tmp0,tmp1;
1294:
1295: if ( NoRA ) {
1296: *h = in; return;
1297: }
1.12 noro 1298: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1299: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1300: }
1301: r = NODE_sortbi(in,0);
1302: n = length(r);
1303: w = (int *)ALLOCA(n*sizeof(int));
1304: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1305: w[i] = (int)BDY(t);
1306: for ( i = 0; i < n; i++ ) {
1307: for ( top = 0, j = n-1; j >= 0; j-- )
1308: if ( j != i ) {
1309: MKNODE(t,(pointer)w[j],top); top = t;
1310: }
1311: get_eg(&tmp0);
1312: if ( PCoeffs )
1313: dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1.12 noro 1314: else {
1315: dpto_dp(ps[w[i]],&p);
1316: _dp_nf_mod_destructive(top,p,ps,m,1,&g);
1317: }
1.1 noro 1318: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12 noro 1319: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1320: fprintf(asir_out,"."); fflush(asir_out);
1321: }
1322: w[i] = newps_mod(g,m);
1323: }
1324: for ( top = 0, j = n-1; j >= 0; j-- ) {
1325: MKNODE(t,(pointer)w[j],top); top = t;
1326: }
1327: *h = top;
1.12 noro 1328: if ( DP_Print || DP_PrintShort )
1.1 noro 1329: fprintf(asir_out,"\n");
1330: }
1331:
1.37 noro 1332: int newps(DP a,int m,NODE subst)
1.1 noro 1333: {
1334: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1335: return -1;
1336: if ( psn == pslen ) {
1337: pslen *= 2;
1338: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1339: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1340: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1341: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1342: if ( m )
1343: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1344: }
1345: if ( Demand ) {
1346: if ( doing_f4 )
1347: ps[psn] = a;
1348: else
1349: ps[psn] = 0;
1.7 noro 1350: dp_save(psn,(Obj)a,0);
1.1 noro 1351: } else
1352: ps[psn] = a;
1353: psh[psn] = BDY(a)->dl;
1354: pss[psn] = a->sugar;
1355: psc[psn] = BDY(a)->c;
1356: if ( m )
1357: _dp_mod(a,m,subst,&psm[psn]);
1.7 noro 1358: if ( GenTrace ) {
1359: NODE tn,tr,tr1;
1.23 noro 1360: LIST trace,trace1;
1361: NODE arg;
1362: Q q1,q2;
1363: STRING fname;
1.37 noro 1364: Obj obj;
1.7 noro 1365:
1366: /* reverse the TraceList */
1367: tn = TraceList;
1368: for ( tr = 0; tn; tn = NEXT(tn) ) {
1369: MKNODE(tr1,BDY(tn),tr); tr = tr1;
1370: }
1371: MKLIST(trace,tr);
1372: if ( OXCheck >= 0 ) {
1373: STOQ(OXCheck,q1);
1374: MKSTR(fname,"check_trace");
1375: STOQ(psn,q2);
1376: arg = mknode(5,q1,fname,a,q2,trace);
1.37 noro 1377: Pox_cmo_rpc(arg,&obj);
1.23 noro 1378: } else if ( OXCheck < 0 ) {
1379: STOQ(psn,q1);
1380: tn = mknode(2,q1,trace);
1381: MKLIST(trace1,tn);
1382: MKNODE(tr,trace1,AllTraceList);
1383: AllTraceList = tr;
1.7 noro 1384: } else
1385: dp_save(psn,(Obj)trace,"t");
1386: TraceList = 0;
1387: }
1.1 noro 1388: return psn++;
1389: }
1390:
1.37 noro 1391: int newps_nosave(DP a,int m,NODE subst)
1.1 noro 1392: {
1393: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1394: return -1;
1395: if ( psn == pslen ) {
1396: pslen *= 2;
1397: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1398: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1399: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1400: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1401: if ( m )
1402: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1403: }
1404: ps[psn] = 0;
1405: psh[psn] = BDY(a)->dl;
1406: pss[psn] = a->sugar;
1407: psc[psn] = BDY(a)->c;
1408: if ( m )
1409: _dp_mod(a,m,subst,&psm[psn]);
1410: return psn++;
1411: }
1412:
1.37 noro 1413: int newps_mod(DP a,int m)
1.1 noro 1414: {
1415: if ( psn == pslen ) {
1416: pslen *= 2;
1417: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1418: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1419: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1420: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P)); /* XXX */
1421: }
1422: ps[psn] = a;
1423: psh[psn] = BDY(ps[psn])->dl;
1424: pss[psn] = ps[psn]->sugar;
1425: return psn++;
1426: }
1427:
1.37 noro 1428: void reducebase_dehomo(NODE f,NODE *g)
1.1 noro 1429: {
1430: int n,i,j,k;
1431: int *r;
1432: DL *w,d;
1433: DP u;
1434: NODE t,top;
1435:
1436: n = length(f);
1437: w = (DL *)ALLOCA(n*sizeof(DL));
1438: r = (int *)ALLOCA(n*sizeof(int));
1439: for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {
1440: r[i] = (int)BDY(t); w[i] = psh[r[i]];
1441: }
1442: for ( i = 0; i < n; i++ ) {
1443: for ( j = 0, d = w[i]; j < n; j++ ) {
1444: if ( j != i ) {
1445: for ( k = 0; k < NVars; k++ )
1446: if ( d->d[k] < w[j]->d[k] )
1447: break;
1448: if ( k == NVars )
1449: break;
1450: }
1451: }
1452: if ( j != n )
1453: r[i] = -1;
1454: }
1455: for ( top = 0, i = n-1; i >= 0; i-- )
1456: if ( r[i] >= 0 ) {
1.7 noro 1457: dp_load(r[i],&ps[r[i]]); dp_dehomo(ps[r[i]],&u);
1458: if ( GenTrace ) {
1459: Q q;
1460: LIST hist;
1461: NODE node;
1462:
1463: STOQ(r[i],q);
1464: node = mknode(4,0,q,0,0);
1465: MKLIST(hist,node);
1466: MKNODE(TraceList,hist,0);
1467: }
1468: j = newps(u,0,0);
1.1 noro 1469: MKNODE(t,j,top); top = t;
1470: }
1471: *g = top;
1472: }
1473:
1.37 noro 1474: NODE append_one(NODE f,int n)
1.1 noro 1475: {
1476: NODE t;
1477:
1478: if ( Reverse || !f ) {
1479: MKNODE(t,(pointer)n,f); return t;
1480: } else {
1481: for ( t = f; NEXT(t); t = NEXT(t) );
1482: MKNODE(NEXT(t),(pointer)n,0);
1483: return f;
1484: }
1485: }
1486:
1.37 noro 1487: DP_pairs minp( DP_pairs d, DP_pairs *prest )
1.1 noro 1488: {
1489: register DP_pairs m, ml, p, l;
1490: register DL lcm;
1491: register int s, nv = CNVars;
1492: register int (*cmpfun)() = cmpdl;
1493:
1494: if ( !(p = NEXT(m = d)) ) {
1495: *prest = p;
1496: NEXT(m) = 0;
1497: return m;
1498: }
1499: for ( lcm = m->lcm, s = m->sugar, ml = 0, l = m; p; p = NEXT(l = p) )
1500: if ( NoSugar ? (*cmpfun)( nv, lcm, p->lcm ) >= 0 :
1501: (s > p->sugar || s == p->sugar && (*cmpfun)( nv, lcm, p->lcm ) >= 0) )
1502: ml = l, lcm = (m = p)->lcm, s = p->sugar;
1503: if ( !ml ) *prest = NEXT(m);
1504: else {
1505: NEXT(ml) = NEXT(m);
1506: *prest = d;
1507: }
1508: NEXT(m) = 0;
1509: return m;
1510: }
1511:
1.37 noro 1512: void minsugar(DP_pairs d,DP_pairs *dm,DP_pairs *dr)
1.1 noro 1513: {
1514: int msugar;
1515: DP_pairs t,dm0,dr0,dmt,drt;
1516:
1517: for ( msugar = d->sugar, t = NEXT(d); t; t = NEXT(t) )
1518: if ( t->sugar < msugar )
1519: msugar = t->sugar;
1520: dm0 = 0; dr0 = 0;
1521: for ( t = d; t; t = NEXT(t) ) {
1522: if ( t->sugar == msugar ) {
1523: NEXTDPP(dm0,dmt);
1524: dmt->dp1 = t->dp1; dmt->dp2 = t->dp2;
1525: dmt->lcm = t->lcm; dmt->sugar = t->sugar;
1526: } else {
1527: NEXTDPP(dr0,drt);
1528: drt->dp1 = t->dp1; drt->dp2 = t->dp2;
1529: drt->lcm = t->lcm; drt->sugar = t->sugar;
1530: }
1531: }
1532: if ( dm0 ) NEXT(dmt) = 0;
1533: if ( dr0 ) NEXT(drt) = 0;
1534: *dm = dm0; *dr = dr0;
1535: }
1536:
1.37 noro 1537: NODE gb(NODE f,int m,NODE subst)
1.1 noro 1538: {
1.42 noro 1539: int i,nh,prev,mag,mag0,magt;
1.1 noro 1540: NODE r,g,gall;
1.37 noro 1541: DP_pairs d;
1.1 noro 1542: DP_pairs l;
1543: DP h,nf,nfm,dp1,dp2;
1544: MP mp;
1.37 noro 1545: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tnp0,tnp1;
1.1 noro 1546: int skip_nf_flag;
1547: double t_0;
1.7 noro 1548: Q q;
1.6 noro 1549: int new_sugar;
1.1 noro 1550: static prev_sugar = -1;
1551:
1552: Max_mag = 0;
1.42 noro 1553: Max_coef = 0;
1.1 noro 1554: prev = 1;
1555: doing_f4 = 0;
1556: if ( m ) {
1557: psm = (DP *)MALLOC(pslen*sizeof(DP));
1558: for ( i = 0; i < psn; i++ )
1559: if ( psh[i] && !validhc(psc[i],m,subst) )
1560: return 0;
1561: else
1562: _dp_mod(ps[i],m,subst,&psm[i]);
1563: }
1564: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1565: i = (int)BDY(r);
1566: d = updpairs(d,g,i);
1567: g = updbase(g,i);
1568: gall = append_one(gall,i);
1569: }
1570: while ( d ) {
1.14 noro 1571: l = minp(d,&d);
1.1 noro 1572: if ( m ) {
1.6 noro 1573: _dp_sp_mod_dup(psm[l->dp1],psm[l->dp2],m,&h);
1.10 noro 1574: if ( h )
1575: new_sugar = h->sugar;
1.1 noro 1576: get_eg(&tnfm0);
1.6 noro 1577: _dp_nf_mod_destructive(gall,h,psm,m,0,&nfm);
1.1 noro 1578: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1579: } else
1580: nfm = (DP)1;
1581: if ( nfm ) {
1582: if ( Demand ) {
1583: if ( dp_load_t(psn,&nf) ) {
1584: skip_nf_flag = 1;
1585: goto skip_nf;
1586: } else {
1587: skip_nf_flag = 0;
1588: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1589: dp_sp(dp1,dp2,&h);
1590: }
1591: } else
1592: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1.7 noro 1593: if ( GenTrace ) {
1594: STOQ(l->dp1,q); ARG1(BDY((LIST)BDY(NEXT(TraceList)))) = q;
1595: STOQ(l->dp2,q); ARG1(BDY((LIST)BDY(TraceList))) = q;
1596: }
1.10 noro 1597: if ( h )
1598: new_sugar = h->sugar;
1.1 noro 1599: get_eg(&tnf0);
1600: t_0 = get_rtime();
1.18 noro 1601: if ( PCoeffs || dp_fcoeffs )
1.16 noro 1602: _dp_nf(gall,h,ps,!Top,&nf);
1603: else
1.21 noro 1604: _dp_nf_z(gall,h,ps,!Top,DP_Multiple,&nf);
1.12 noro 1605: if ( DP_Print )
1.1 noro 1606: fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
1607: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1608: } else
1609: nf = 0;
1610: skip_nf:
1611: if ( nf ) {
1612: NZR++;
1613: get_eg(&tpz0);
1614: prim_part(nf,0,&h);
1615: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1616: get_eg(&tnp0);
1617: if ( Demand && skip_nf_flag )
1618: nh = newps_nosave(h,m,subst);
1619: else
1620: nh = newps(h,m,subst);
1621: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1622: if ( nh < 0 )
1623: return 0;
1624: d = updpairs(d,g,nh);
1625: g = updbase(g,nh);
1626: gall = append_one(gall,nh);
1627: if ( !dp_fcoeffs && ShowMag ) {
1.42 noro 1628: for ( mag = 0, mag0 = 0, mp = BDY(h); mp; mp = NEXT(mp) ) {
1629: magt = p_mag((P)mp->c);
1630: mag0 = MAX(mag0,magt);
1631: mag += magt;
1632: }
1633: Max_coef = MAX(Max_coef,mag0);
1.1 noro 1634: Max_mag = MAX(Max_mag,mag);
1635: }
1.12 noro 1636: if ( DP_Print ) {
1.1 noro 1637: if ( !prev )
1638: fprintf(asir_out,"\n");
1639: print_split_e(&tnf0,&tnf1); print_split_e(&tpz0,&tpz1);
1640: printdl(psh[nh]);
1641: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",
1642: l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
1643: pss[nh]);
1644: if ( ShowMag )
1.42 noro 1645: fprintf(asir_out,",mag=(%d,%d)",mag,mag0);
1.1 noro 1646: fprintf(asir_out,"\n"); fflush(asir_out);
1.12 noro 1647: } else if ( DP_PrintShort ) {
1.1 noro 1648: fprintf(asir_out,"+"); fflush(asir_out);
1649: }
1650: prev = 1;
1651: } else {
1652: if ( m )
1653: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1654: ZR++;
1.12 noro 1655: if ( DP_Print || DP_PrintShort ) {
1.6 noro 1656: if ( new_sugar != prev_sugar ) {
1657: fprintf(asir_out,"[%d]",new_sugar);
1658: prev_sugar = new_sugar;
1.1 noro 1659: }
1660: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1661: }
1662: }
1663: }
1.12 noro 1664: if ( DP_Print || DP_PrintShort )
1.1 noro 1665: fprintf(asir_out,"gb done\n");
1666: return g;
1667: }
1668:
1.37 noro 1669: NODE gb_mod(NODE f,int m)
1.1 noro 1670: {
1671: int i,nh,prev;
1672: NODE r,g,gall;
1.37 noro 1673: DP_pairs d;
1.1 noro 1674: DP_pairs l;
1675: DP h,nf;
1.37 noro 1676: struct oEGT tnfm0,tnfm1,tpz0,tpz1;
1.1 noro 1677:
1678: prev = 1;
1679: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1680: i = (int)BDY(r);
1681: d = updpairs(d,g,i);
1682: g = updbase(g,i);
1683: gall = append_one(gall,i);
1684: }
1685: while ( d ) {
1.14 noro 1686: l = minp(d,&d);
1.1 noro 1687: if ( PCoeffs ) {
1688: dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1.14 noro 1689: get_eg(&tnfm0);
1.1 noro 1690: dp_nf_mod(gall,h,ps,m,!Top,&nf);
1691: } else {
1.6 noro 1692: _dp_sp_mod_dup(ps[l->dp1],ps[l->dp2],m,&h);
1.14 noro 1693: get_eg(&tnfm0);
1.6 noro 1694: _dp_nf_mod_destructive(gall,h,ps,m,!Top,&nf);
1.1 noro 1695: }
1696: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1697: if ( nf ) {
1698: NZR++;
1699: get_eg(&tpz0);
1700: prim_part(nf,m,&h);
1701: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1702: nh = newps_mod(h,m);
1703: if ( nh < 0 )
1704: return 0;
1705: d = updpairs(d,g,nh);
1706: g = updbase(g,nh);
1707: gall = append_one(gall,nh);
1.12 noro 1708: if ( DP_Print ) {
1.1 noro 1709: if ( !prev )
1710: fprintf(asir_out,"\n");
1711: print_split_eg(&tnfm0,&tnfm1); fflush(asir_out);
1712: 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]);
1713: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1.12 noro 1714: } else if ( DP_PrintShort ) {
1.6 noro 1715: fprintf(asir_out,"+"); fflush(asir_out);
1.1 noro 1716: }
1717: prev = 1;
1718: } else {
1719: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1720: ZR++;
1.12 noro 1721: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1722: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1723: }
1724: }
1725: }
1.12 noro 1726: if ( DP_Print || DP_PrintShort )
1.1 noro 1727: fprintf(asir_out,"gb_mod done\n");
1728: return g;
1729: }
1730:
1.37 noro 1731: DP_pairs updpairs( DP_pairs d, NODE /* of index */ g, int t)
1.1 noro 1732: {
1.14 noro 1733: register DP_pairs d1, dd, nd;
1734: int dl,dl1;
1.1 noro 1735:
1.14 noro 1736: if ( !g ) return d;
1737: if ( !NoCriB && d ) {
1738: dl = DPPlength(d);
1739: d = criterion_B( d, t );
1740: dl -= DPPlength(d); NBP += dl;
1.1 noro 1741: }
1.14 noro 1742: d1 = newpairs( g, t );
1743: if ( NEXT(d1) ) {
1744: dl = DPPlength(d1); TP += dl;
1745: d1 = criterion_M( d1 );
1746: dl1 = DPPlength(d1); NMP += (dl-dl1); dl = dl1;
1747: d1 = criterion_F( d1 );
1748: dl1 = DPPlength(d1); NFP += (dl-dl1); dl = dl1;
1749: } else
1750: dl = 1;
1751: if ( !do_weyl )
1752: for ( dd = 0; d1; d1 = nd ) {
1753: nd = NEXT(d1);
1754: if ( !criterion_2( d1->dp1, d1->dp2 ) ) {
1755: NEXT(d1) = dd;
1756: dd = d1;
1.1 noro 1757: }
1758: }
1.14 noro 1759: else
1760: dd = d1;
1761: dl1 = DPPlength(dd); NDP += (dl-dl1);
1762: if ( !(nd = d) ) return dd;
1763: while ( nd = NEXT(d1 = nd) ) ;
1764: NEXT(d1) = dd;
1765: return d;
1.1 noro 1766: }
1767:
1.37 noro 1768: DP_pairs newpairs( NODE /* of index */ g, int t )
1.14 noro 1769: {
1770: register NODE r;
1771: register DL tdl = psh[t];
1772: register int ts;
1773: register DP_pairs p, last;
1774: int dp;
1775: register DL dl;
1776: register int s;
1.1 noro 1777:
1.14 noro 1778: ts = pss[t] - tdl->td;
1779: for ( last = 0, r = g; r; r = NEXT(r) ) {
1780: NEXT(p = NEWDP_pairs) = last;
1781: last = p;
1782: dp = p->dp1 = (int)BDY(r); p->dp2 = t;
1783: p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
1784: #if 0
1785: if ( do_weyl )
1.41 noro 1786: p->sugar = dl_weyl_weight(p->lcm);
1.14 noro 1787: else
1788: #endif
1789: p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
1790: }
1791: return last;
1792: }
1.1 noro 1793:
1.37 noro 1794: DP_pairs criterion_B( DP_pairs d, int s )
1.1 noro 1795: {
1.14 noro 1796: register DP_pairs dd, p;
1797: register DL tij, t = psh[s], dltmp;
1.1 noro 1798:
1.14 noro 1799: if ( !d ) return 0;
1800: NEWDL( dltmp, CNVars );
1801: for ( dd = 0; d; d = p ) {
1802: p = NEXT(d),
1803: tij = d->lcm;
1804: if ( tij->td != lcm_of_DL(CNVars, tij, t, dltmp )->td
1805: || !dl_equal(CNVars, tij, dltmp )
1806: || (tij->td == lcm_of_DL(CNVars, psh[d->dp1], t, dltmp )->td
1807: && dl_equal(CNVars, dltmp, tij ))
1808: || (tij->td == lcm_of_DL(CNVars, psh[d->dp2], t, dltmp )->td
1809: && dl_equal(CNVars, dltmp, tij )) ) {
1810: NEXT(d) = dd;
1811: dd = d;
1812: }
1.1 noro 1813: }
1.14 noro 1814: return dd;
1815: }
1.1 noro 1816:
1.37 noro 1817: DP_pairs criterion_M( DP_pairs d1 )
1.14 noro 1818: {
1819: register DP_pairs dd, e, d3, d2, p;
1820: register DL itdl, jtdl;
1821: register int itdltd, jtdltd;
1.1 noro 1822:
1.14 noro 1823: for ( dd = 0, e = d1; e; e = d3 ) {
1824: if ( !(d2 = NEXT(e)) ) {
1825: NEXT(e) = dd;
1826: return e;
1827: }
1828: itdltd = (itdl = e->lcm)->td;
1829: for ( d3 = 0; d2; d2 = p ) {
1830: p = NEXT(d2),
1831: jtdltd = (jtdl = d2->lcm)->td;
1832: if ( jtdltd == itdltd )
1833: if ( dl_equal(CNVars, itdl, jtdl ) ) ;
1834: else if ( dl_redble( jtdl, itdl ) ) continue;
1835: else if ( dl_redble( itdl, jtdl ) ) goto delit;
1836: else ;
1837: else if ( jtdltd > itdltd )
1838: if ( dl_redble( jtdl, itdl ) ) continue;
1.1 noro 1839: else ;
1840: else if ( dl_redble( itdl, jtdl ) ) goto delit;
1841: NEXT(d2) = d3;
1842: d3 = d2;
1843: }
1844: NEXT(e) = dd;
1845: dd = e;
1846: continue;
1847: /**/
1848: delit: NEXT(d2) = d3;
1849: d3 = d2;
1850: for ( ; p; p = d2 ) {
1851: d2 = NEXT(p);
1852: NEXT(p) = d3;
1853: d3 = p;
1854: }
1855: }
1856: return dd;
1857: }
1858:
1.37 noro 1859: static DP_pairs collect_pairs_of_hdlcm( DP_pairs d1, DP_pairs *prest )
1.1 noro 1860: {
1861: register DP_pairs w, p, r, s;
1862: register DL ti;
1863: register int td;
1864:
1865: td = (ti = (w = d1)->lcm)->td;
1866: s = NEXT(w);
1867: NEXT(w) = 0;
1868: for ( r = 0; s; s = p ) {
1869: p = NEXT(s);
1870: if ( td == s->lcm->td && dl_equal(CNVars, ti, s->lcm ) )
1871: {
1872: NEXT(s) = w;
1873: w = s;
1874: } else {
1875: NEXT(s) = r;
1876: r = s;
1877: }
1878: }
1879: *prest = r;
1880: return w;
1881: }
1882:
1.37 noro 1883: int criterion_2( int dp1, int dp2 )
1.1 noro 1884: {
1885: register int i, *d1, *d2;
1886:
1887: d1 = psh[dp1]->d, d2 = psh[dp2]->d;
1888: for ( i = CNVars; --i >= 0; d1++, d2++ )
1889: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
1890: return 1;
1891: }
1892:
1.37 noro 1893: DP_pairs criterion_F( DP_pairs d1 )
1.1 noro 1894: {
1895: DP_pairs rest, head;
1896: register DP_pairs last, p, r, w;
1897: register int s;
1898:
1899: for ( head = last = 0, p = d1; NEXT(p); ) {
1900: s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
1901: while ( w = NEXT(w) )
1902: if ( criterion_2( w->dp1, w->dp2 ) ) {
1903: r = w;
1904: break;
1905: } else if ( w->sugar < s ) s = (r = w)->sugar;
1906: if ( last ) NEXT(last) = r;
1907: else head = r;
1908: NEXT(last = r) = 0;
1909: if ( !(p = rest) ) return head;
1910: }
1911: if ( !last ) return p;
1912: NEXT(last) = p;
1913: return head;
1914: }
1915:
1.37 noro 1916: NODE updbase(NODE g,int t)
1.1 noro 1917: {
1918: g = remove_reducibles(g,t);
1919: g = append_one(g,t);
1920: return g;
1921: }
1922:
1.37 noro 1923: NODE /* of index */ remove_reducibles(NODE /* of index */ nd, int newdp )
1.1 noro 1924: {
1925: register DL dl, dln;
1926: register NODE last, p, head;
1927: register int td;
1928:
1929: dl = psh[newdp];
1930: td = dl->td;
1931: for ( head = last = 0, p = nd; p; ) {
1932: dln = psh[(int)BDY(p)];
1933: if ( dln->td >= td && dl_redble( dln, dl ) ) {
1934: p = NEXT(p);
1935: if ( last ) NEXT(last) = p;
1936: } else {
1937: if ( !last ) head = p;
1938: p = NEXT(last = p);
1939: }
1940: }
1941: return head;
1942: }
1943:
1.37 noro 1944: int dl_redble(DL dl1,DL dl2)
1.1 noro 1945: {
1946: register int n, *d1, *d2;
1947:
1948: for ( d1 = dl1->d, d2 = dl2->d, n = CNVars; --n >= 0; d1++, d2++ )
1949: if ( *d1 < *d2 ) return 0;
1950: return 1;
1.5 noro 1951: }
1952:
1.41 noro 1953: #if 0
1954: int dl_weyl_weight(DL dl)
1.5 noro 1955: {
1956: int n,w,i;
1957:
1958: n = CNVars/2;
1959: for ( i = 0, w = 0; i < n; i++ )
1960: w += (-dl->d[i]+dl->d[n+i]);
1961: return w;
1.1 noro 1962: }
1.41 noro 1963: #endif
1.1 noro 1964:
1.37 noro 1965: int gbcheck(NODE f)
1.1 noro 1966: {
1967: int i;
1968: NODE r,g,gall;
1969: DP_pairs d,l;
1970: DP h,nf,dp1,dp2;
1971: struct oEGT tmp0,tmp1;
1972:
1973: if ( NoGC )
1974: return 1;
1975: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1976: i = (int)BDY(r);
1977: d = updpairs(d,g,i);
1978: g = updbase(g,i);
1979: gall = append_one(gall,i);
1980: }
1.12 noro 1981: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1982: fprintf(asir_out,"gbcheck total %d pairs\n",DPPlength(d)); fflush(asir_out);
1983: }
1984: while ( d ) {
1985: l = d; d = NEXT(d);
1986: get_eg(&tmp0);
1.17 noro 1987: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1988: dp_sp(dp1,dp2,&h);
1989: /* fprintf(stderr,"{%d,%d}",l->dp1,l->dp2); */
1.16 noro 1990: _dp_nf(gall,h,ps,1,&nf);
1.1 noro 1991: get_eg(&tmp1); add_eg(&eg_gc,&tmp0,&tmp1);
1.12 noro 1992: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1993: fprintf(asir_out,"."); fflush(asir_out);
1994: }
1995: if ( nf )
1996: return 0;
1997: }
1.12 noro 1998: if ( DP_Print || DP_PrintShort )
1.1 noro 1999: fprintf(asir_out,"\n");
2000: return 1;
1.38 noro 2001: }
2002:
1.40 noro 2003: void gbcheck_list(NODE f,int n,VECT *gp,LIST *pp)
1.38 noro 2004: {
2005: int i;
2006: NODE r,g,gall,u,u0,t;
1.39 noro 2007: VECT vect;
1.38 noro 2008: LIST pair;
2009: DP_pairs d,l;
2010: Q q1,q2;
2011:
1.40 noro 2012: /* we need the following settings */
2013: NVars = CNVars = n;
1.39 noro 2014: setup_arrays(f,0,&r);
2015: for ( gall = g = 0, d = 0; r; r = NEXT(r) ) {
1.38 noro 2016: i = (int)BDY(r);
2017: d = updpairs(d,g,i);
2018: g = updbase(g,i);
2019: gall = append_one(gall,i);
2020: }
1.39 noro 2021: NEWVECT(vect); vect->len = psn; vect->body = (pointer)ps;
2022: *gp = vect;
2023:
1.38 noro 2024: for ( u0 = 0, l = d; l; l = NEXT(l) ) {
2025: NEXTNODE(u0,u);
2026: STOQ(l->dp1,q1);
2027: STOQ(l->dp2,q2);
2028: t = mknode(2,q1,q2);
2029: MKLIST(pair,t);
2030: BDY(u) = (pointer)pair;
2031: }
2032: if ( u0 )
2033: NEXT(u) = 0;
1.39 noro 2034: MKLIST(*pp,u0);
1.1 noro 2035: }
2036:
1.37 noro 2037: int membercheck(NODE f,NODE x)
1.1 noro 2038: {
2039: DP g;
2040: struct oEGT tmp0,tmp1;
2041:
2042: if ( NoMC )
2043: return 1;
1.12 noro 2044: if ( DP_Print || DP_PrintShort ) {
1.1 noro 2045: fprintf(asir_out,"membercheck\n"); fflush(asir_out);
2046: }
2047: for ( ; f; f = NEXT(f) ) {
2048: get_eg(&tmp0);
1.16 noro 2049: _dp_nf(x,(DP)BDY(f),ps,1,&g);
1.1 noro 2050: get_eg(&tmp1); add_eg(&eg_mc,&tmp0,&tmp1);
1.12 noro 2051: if ( DP_Print ) {
1.1 noro 2052: print_split_eg(&tmp0,&tmp1); fflush(asir_out);
1.12 noro 2053: } else if ( DP_PrintShort ) {
1.1 noro 2054: fprintf(asir_out,"."); fflush(asir_out);
2055: }
2056: if ( g )
2057: return 0;
2058: }
1.12 noro 2059: if ( DP_Print || DP_PrintShort )
1.1 noro 2060: fprintf(asir_out,"\n");
2061: return 1;
2062: }
2063:
1.37 noro 2064: void dp_set_flag(Obj name,Obj value)
1.1 noro 2065: {
2066: char *n;
2067: int v;
1.43 noro 2068: Q ratio;
1.1 noro 2069:
2070: if ( OID(name) != O_STR )
2071: return;
2072: n = BDY((STRING)name);
2073: if ( !strcmp(n,"Demand") ) {
2074: Demand = value ? BDY((STRING)value) : 0; return;
2075: }
2076: if ( !strcmp(n,"Dist") ) {
2077: Dist = (LIST)value; return;
2078: }
1.43 noro 2079: if ( !strcmp(n,"Content") ) {
2080: ratio = (Q)value;
2081: if ( ratio ) {
2082: DP_Multiple = BD(NM(ratio))[0];
2083: Denominator = INT(ratio) ? 1 : BD(DN(ratio))[0];
2084: } else {
2085: DP_Multiple = 0;
2086: Denominator = 1;
2087: }
2088: }
1.1 noro 2089: if ( value && OID(value) != O_N )
2090: return;
2091: v = QTOS((Q)value);
2092: if ( !strcmp(n,"NoSugar") )
2093: NoSugar = v;
2094: else if ( !strcmp(n,"NoCriB") )
2095: NoCriB = v;
2096: else if ( !strcmp(n,"NoGC") )
2097: NoGC = v;
2098: else if ( !strcmp(n,"NoMC") )
2099: NoMC = v;
2100: else if ( !strcmp(n,"NoRA") )
2101: NoRA = v;
2102: else if ( !strcmp(n,"NoGCD") )
2103: NoGCD = v;
2104: else if ( !strcmp(n,"Top") )
2105: Top = v;
2106: else if ( !strcmp(n,"ShowMag") )
2107: ShowMag = v;
1.14 noro 2108: else if ( !strcmp(n,"PrintShort") )
1.12 noro 2109: DP_PrintShort = v;
1.14 noro 2110: else if ( !strcmp(n,"Print") )
1.12 noro 2111: DP_Print = v;
1.14 noro 2112: else if ( !strcmp(n,"NFStat") )
2113: DP_NFStat = v;
1.1 noro 2114: else if ( !strcmp(n,"Stat") )
2115: Stat = v;
2116: else if ( !strcmp(n,"Reverse") )
2117: Reverse = v;
2118: else if ( !strcmp(n,"Multiple") )
1.14 noro 2119: DP_Multiple = v;
1.1 noro 2120: else if ( !strcmp(n,"Denominator") )
2121: Denominator = v;
2122: else if ( !strcmp(n,"PtozpRA") )
2123: PtozpRA = v;
1.7 noro 2124: else if ( !strcmp(n,"GenTrace") )
2125: GenTrace = v;
2126: else if ( !strcmp(n,"OXCheck") )
2127: OXCheck = v;
1.1 noro 2128: }
2129:
1.37 noro 2130: void dp_make_flaglist(LIST *list)
1.1 noro 2131: {
1.43 noro 2132: Q v,nm,dn;
1.1 noro 2133: STRING name,path;
2134: NODE n,n1;
2135:
1.43 noro 2136: #if 0
1.14 noro 2137: STOQ(DP_Multiple,v); MKNODE(n,v,0); MKSTR(name,"DP_Multiple"); MKNODE(n1,name,n); n = n1;
1.1 noro 2138: STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
1.43 noro 2139: #else
2140: if ( DP_Multiple ) {
2141: STOQ(DP_Multiple,nm); STOQ(Denominator,dn); divq(nm,dn,&v);
2142: } else
2143: v = 0;
2144: MKNODE(n,v,0); MKSTR(name,"Content"); MKNODE(n1,name,n); n = n1;
2145: #endif
1.1 noro 2146: MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
2147: STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
2148: STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
1.14 noro 2149: STOQ(DP_Print,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Print"); MKNODE(n1,name,n); n = n1;
1.20 noro 2150: STOQ(DP_PrintShort,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PrintShort"); MKNODE(n1,name,n); n = n1;
1.14 noro 2151: STOQ(DP_NFStat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NFStat"); MKNODE(n1,name,n); n = n1;
1.7 noro 2152: STOQ(OXCheck,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"OXCheck"); MKNODE(n1,name,n); n = n1;
2153: STOQ(GenTrace,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenTrace"); MKNODE(n1,name,n); n = n1;
1.1 noro 2154: STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
2155: STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
2156: STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
2157: STOQ(NoGCD,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGCD"); MKNODE(n1,name,n); n = n1;
2158: STOQ(NoRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoRA"); MKNODE(n1,name,n); n = n1;
2159: STOQ(NoMC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoMC"); MKNODE(n1,name,n); n = n1;
2160: STOQ(NoGC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGC"); MKNODE(n1,name,n); n = n1;
2161: STOQ(NoCriB,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoCriB"); MKNODE(n1,name,n); n = n1;
2162: STOQ(NoSugar,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoSugar"); MKNODE(n1,name,n); n = n1;
2163: if ( Demand )
2164: MKSTR(path,Demand);
2165: else
2166: path = 0;
2167: MKNODE(n1,path,n); n = n1; MKSTR(name,"Demand"); MKNODE(n1,name,n); n = n1;
2168: MKLIST(*list,n);
2169: }
2170:
2171: #define DELIM '/'
2172:
1.37 noro 2173: void dp_save(int index,Obj p,char *prefix)
1.1 noro 2174: {
2175: FILE *fp;
2176: char path[BUFSIZ];
2177:
1.7 noro 2178: if ( prefix )
2179: sprintf(path,"%s%c%s%d",Demand,DELIM,prefix,index);
2180: else
2181: sprintf(path,"%s%c%d",Demand,DELIM,index);
1.1 noro 2182: if ( !(fp = fopen(path,"wb") ) )
2183: error("dp_save : cannot open a file");
1.7 noro 2184: savevl(fp,VC); saveobj(fp,p); fclose(fp);
1.1 noro 2185: }
2186:
1.37 noro 2187: void dp_load(int index,DP *p)
1.1 noro 2188: {
2189: FILE *fp;
2190: char path[BUFSIZ];
2191:
2192: if ( !Demand || ps[index] )
2193: *p = ps[index];
2194: else {
2195: sprintf(path,"%s%c%d",Demand,DELIM,index);
2196: if ( !(fp = fopen(path,"rb") ) )
2197: error("dp_load : cannot open a file");
1.45 noro 2198: if ( PCoeffs )
2199: loadvl(fp);
2200: else
2201: skipvl(fp);
2202: loadobj(fp,(Obj *)p); fclose(fp);
1.1 noro 2203: }
2204: }
2205:
1.37 noro 2206: int dp_load_t(int index,DP *p)
1.1 noro 2207: {
2208: FILE *fp;
2209: char path[BUFSIZ];
2210:
2211: sprintf(path,"%s%c%d",Demand,DELIM,index);
2212: if ( !(fp = fopen(path,"rb") ) )
2213: return 0;
2214: else {
1.45 noro 2215: if ( PCoeffs )
2216: loadvl(fp);
2217: else
2218: skipvl(fp);
2219: loadobj(fp,(Obj *)p); fclose(fp); return 1;
1.1 noro 2220: }
2221: }
2222:
2223: void init_stat() {
2224: init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
1.14 noro 2225: init_eg(&eg_pz); init_eg(&eg_np);
1.1 noro 2226: init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
1.35 noro 2227: ZR = NZR = TP = NMP = NBP = NFP = NDP = 0;
1.1 noro 2228: }
2229:
2230: void print_stat() {
1.12 noro 2231: if ( !DP_Print && !Stat )
1.1 noro 2232: return;
2233: print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
1.14 noro 2234: print_eg("PZ",&eg_pz); print_eg("NP",&eg_np);
1.1 noro 2235: print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
2236: 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);
2237: }
1.14 noro 2238:
2239: /*
2240: * dp_nf used in gb()
2241: *
2242: */
2243:
2244: double pz_t_e, pz_t_d, pz_t_d1, pz_t_c, im_t_s, im_t_r;
2245:
2246: extern int GenTrace;
2247: extern NODE TraceList;
1.15 noro 2248: extern int mpi_mag;
2249:
1.37 noro 2250: void dp_mulc_d(DP p,P c,DP *r)
1.15 noro 2251: {
2252: if ( Dist && BDY(Dist)
2253: && HMAG(p) > mpi_mag
2254: && p_mag((P)c) > mpi_mag ) {
2255: if ( DP_NFStat ) fprintf(asir_out,"~");
2256: dp_imul_d(p,(Q)c,r);
2257: } else {
2258: if ( DP_NFStat ) fprintf(asir_out,"_");
2259: muldc(CO,p,c,r);
2260: }
2261: }
1.14 noro 2262:
1.37 noro 2263: void _dp_nf(NODE b,DP g,DP *ps,int full,DP *rp)
1.16 noro 2264: {
2265: DP u,p,d,s,t,mult;
2266: P coef;
2267: NODE l;
2268: MP m,mr;
2269: int sugar,psugar;
2270:
2271: if ( !g ) {
2272: *rp = 0; return;
2273: }
2274: sugar = g->sugar;
2275: for ( d = 0; g; ) {
2276: for ( u = 0, l = b; l; l = NEXT(l) ) {
2277: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
2278: dp_load((int)BDY(l),&p);
2279: /* t+u = coef*(d+g) - mult*p (t = coef*d) */
2280: dp_red(d,g,p,&t,&u,&coef,&mult);
2281: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
2282: sugar = MAX(sugar,psugar);
2283: if ( GenTrace ) {
2284: LIST hist;
2285: Q cq;
2286: NODE node,node0;
2287:
2288: STOQ((int)BDY(l),cq);
2289: node0 = mknode(4,coef,cq,mult,ONE);
2290: MKLIST(hist,node0);
2291: MKNODE(node,hist,TraceList); TraceList = node;
2292: }
2293: if ( !u ) {
2294: if ( d )
2295: d->sugar = sugar;
2296: *rp = d; return;
2297: }
2298: d = t;
2299: break;
2300: }
2301: }
2302: if ( u )
2303: g = u;
2304: else if ( !full ) {
2305: if ( g ) {
2306: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
2307: }
2308: *rp = g; return;
2309: } else {
2310: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
2311: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
2312: addd(CO,d,t,&s); d = s;
2313: dp_rest(g,&t); g = t;
2314: }
2315: }
2316: if ( d )
2317: d->sugar = sugar;
2318: *rp = d;
2319: }
2320:
1.37 noro 2321: void _dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *r)
1.14 noro 2322: {
1.37 noro 2323: DP u,dp,rp,t,t1,red,shift;
1.14 noro 2324: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
2325: NODE l;
2326: int hmag,denom;
2327: int sugar,psugar;
2328: STRING imul;
1.37 noro 2329: double t_0,tt,t_p,t_m,t_g,t_a;
1.14 noro 2330: LIST hist;
2331: NODE node;
2332: Q rcred,mrcred;
2333:
2334: if ( !g ) {
2335: *r = 0; return;
2336: }
2337: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1.15 noro 2338: t_p = t_m = t_g = t_a = 0;
1.14 noro 2339:
2340: denom = Denominator?Denominator:1;
2341: hmag = multiple*HMAG(g)/denom;
2342: sugar = g->sugar;
2343:
2344: dc = 0; dp = 0; rc = ONE; rp = g;
2345: MKSTR(imul,"dp_imul_index");
2346:
2347: /* g = dc*dp+rc*rp */
2348: for ( ; rp; ) {
2349: for ( u = 0, l = b; l; l = NEXT(l) ) {
2350: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
2351: t_0 = get_rtime();
2352: dp_load((int)BDY(l),&red);
2353: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
2354: igcd_cofactor((Q)BDY(rp)->c,(Q)BDY(red)->c,&gcd,&cred,&cr);
1.15 noro 2355: tt = get_rtime(); t_p += tt-t_0;
1.14 noro 2356:
2357: dp_subd(rp,red,&shift);
1.37 noro 2358: dp_mulc_d(rp,(P)cr,&t);
1.14 noro 2359: chsgnp((P)cred,(P *)&mcred);
1.37 noro 2360: dp_mulc_d(red,(P)mcred,&t1);
1.15 noro 2361: muld(CO,shift,t1,&t1);
2362: addd(CO,t,t1,&u);
2363: t_m += get_rtime()-tt;
1.14 noro 2364:
2365: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
2366: sugar = MAX(sugar,psugar);
1.15 noro 2367:
1.14 noro 2368: if ( GenTrace ) {
2369: /* u = cr*rp + (-cred)*shift*red */
2370: STOQ((int)BDY(l),cq);
2371: node = mknode(4,cr,cq,0,0);
2372: mulq(cred,rc,&rcred);
2373: chsgnnum((Num)rcred,(Num *)&mrcred);
2374: muldc(CO,shift,(P)mrcred,(DP *)&ARG2(node));
2375: MKLIST(hist,node);
2376: }
1.15 noro 2377:
1.14 noro 2378: if ( !u ) {
2379: if ( dp )
2380: dp->sugar = sugar;
2381: *r = dp;
2382: if ( GenTrace ) {
2383: ARG3(BDY(hist)) = ONE;
2384: MKNODE(node,hist,TraceList); TraceList = node;
2385: }
2386: goto final;
2387: }
2388: break;
2389: }
2390: }
2391: if ( u ) {
2392: if ( multiple && HMAG(u) > hmag ) {
2393: t_0 = get_rtime();
1.16 noro 2394: dp_ptozp_d(u,&rp);
1.15 noro 2395: tt = get_rtime(); t_g += tt-t_0;
2396:
2397: divsq((Q)BDY(u)->c,(Q)BDY(rp)->c,&cont);
1.14 noro 2398: if ( !dp_fcoeffs && DP_NFStat ) {
1.15 noro 2399: fprintf(asir_out,
2400: "(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1.14 noro 2401: fflush(asir_out);
2402: }
1.15 noro 2403: mulq(cr,dc,&dcq); mulq(cont,rc,&rcq);
1.14 noro 2404: igcd_cofactor(dcq,rcq,&gcd,&dc,&rc);
1.15 noro 2405: t_a = get_rtime()-tt;
2406:
1.14 noro 2407: hmag = multiple*HMAG(rp)/denom;
2408: if ( GenTrace ) {
2409: ARG3(BDY(hist)) = (pointer)gcd;
2410: MKNODE(node,hist,TraceList); TraceList = node;
2411: }
2412: } else {
1.15 noro 2413: rp = u;
1.14 noro 2414: t_0 = get_rtime();
1.15 noro 2415: mulq(cr,dc,&dc);
2416: t_a += get_rtime()-t_0;
1.14 noro 2417: if ( GenTrace ) {
2418: ARG3(BDY(hist)) = (pointer)ONE;
2419: MKNODE(node,hist,TraceList); TraceList = node;
2420: }
2421: }
2422: } else if ( !full ) {
2423: if ( rp ) {
2424: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
2425: }
2426: *r = rp;
2427: goto final;
2428: } else {
2429: t_0 = get_rtime();
2430: mulq((Q)BDY(rp)->c,rc,&c);
1.15 noro 2431: igcd_cofactor(dc,c,&dc,&dcq,&cq);
2432: muldc(CO,dp,(P)dcq,&t);
2433: dp_hm(rp,&t1); BDY(t1)->c = (P)cq; addd(CO,t,t1,&dp);
2434: dp_rest(rp,&rp);
2435: t_a += get_rtime()-t_0;
1.14 noro 2436: }
2437: }
2438: if ( GenTrace ) {
2439: mulq(ARG3(BDY((LIST)BDY(TraceList))),dc,&cq);
2440: ARG3(BDY((LIST)BDY(TraceList))) = (pointer)cq;
2441: }
2442: if ( dp )
2443: dp->sugar = sugar;
2444: *r = dp;
2445: final:
2446: if ( DP_NFStat )
1.15 noro 2447: fprintf(asir_out,
2448: "(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
2449: t_p,t_m,t_g,t_a,
1.14 noro 2450: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
2451: }
2452:
2453: void imulv();
2454:
1.37 noro 2455: void dp_imul_d(DP p,Q q,DP *rp)
1.14 noro 2456: {
2457: int nsep,ndist,i,j,k,l,n;
2458: double t0,t1,t2;
2459: Q *s;
2460: pointer *b;
2461: VECT c,cs,ri;
2462: VECT *r;
2463: MP m;
1.37 noro 2464: NODE tn,dist,n0;
2465: Obj obj;
1.14 noro 2466: STRING imul;
2467: extern LIST Dist;
2468:
2469: if ( !p || !q ) {
2470: *rp = 0; return;
2471: }
2472: dist = BDY(Dist);
2473: for ( tn = dist, ndist = 0; tn; tn = NEXT(tn), ndist++ );
2474: nsep = ndist + 1;
2475: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
2476: if ( n <= nsep ) {
2477: muldc(CO,p,(P)q,rp); return;
2478: }
2479: MKSTR(imul,"imulv");
2480: t0 = get_rtime();
2481: dp_dtov(p,&c);
2482: sepvect(c,nsep,&cs);
2483: r = (VECT *)CALLOC(nsep,sizeof(VECT *));
2484: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
2485: n0 = mknode(4,BDY(tn),imul,b[i],q);
1.37 noro 2486: Pox_rpc(n0,&obj);
1.14 noro 2487: }
2488: t1 = get_rtime();
2489: im_t_s += t1 - t0;
2490: imulv(b[i],q,&r[i]);
2491: t1 = get_rtime();
2492: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
2493: MKNODE(n0,BDY(tn),0);
1.37 noro 2494: Pox_pop_local(n0,&obj); r[i] = (VECT)obj;
1.14 noro 2495: if ( OID(r[i]) == O_ERR ) {
2496: printexpr(CO,(Obj)r[i]);
2497: error("dp_imul_d : aborted");
2498: }
2499: }
2500: t2 = get_rtime();
2501: im_t_r += t2 - t1;
2502: s = (Q *)CALLOC(n,sizeof(Q));
2503: for ( i = j = 0; i < nsep; i++ ) {
2504: for ( k = 0, ri = r[i], l = ri->len; k < l; k++, j++ ) {
2505: s[j] = (Q)BDY(ri)[k];
2506: }
2507: }
2508: dp_vtod(s,p,rp);
2509: }
2510:
1.37 noro 2511: void imulv(VECT w,Q c,VECT *rp)
1.14 noro 2512: {
2513: int n,i;
2514: VECT r;
2515:
2516: n = w->len;
2517: MKVECT(r,n); *rp = r;
2518: for ( i = 0; i < n; i++ )
2519: mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
1.32 noro 2520: }
2521:
1.37 noro 2522: void dptoca(DP p,unsigned int **rp)
1.32 noro 2523: {
2524: int i;
2525: MP m;
2526: unsigned int *r;
2527:
2528: if ( !p )
2529: *rp = 0;
2530: else {
2531: for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ );
2532: *rp = r = (unsigned int *)MALLOC_ATOMIC(i*sizeof(unsigned int));
2533: for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ )
2534: r[i] = ITOS(C(m));
2535: }
1.14 noro 2536: }
2537:
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