Annotation of OpenXM_contrib2/asir2000/builtin/gr.c, Revision 1.50
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.50 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.49 2003/06/21 02:09:15 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.50 ! noro 622: struct oEGT tmp0,tmp1,eg_split_symb,eg_split_conv,eg_split_elim1,eg_split_elim2;
! 623: extern struct oEGT eg_symb,eg_conv,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.50 ! noro 629: init_eg(&eg_symb); init_eg(&eg_conv); 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.50 ! noro 676: get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
! 677: init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
! 678:
1.35 noro 679: if ( DP_Print )
680: fprintf(asir_out,"number of reducers : %d\n",nred);
1.1 noro 681:
1.50 ! noro 682: get_eg(&tmp0);
1.1 noro 683: /* the first nred polys in blist are reducers */
684: /* row = the number of all the polys */
685: for ( r = blist, row = 0; r; r = NEXT(r), row++ );
1.5 noro 686:
687: /* col = number of all terms */
1.29 noro 688: for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.5 noro 689:
690: /* head terms of all terms */
1.1 noro 691: at = (DL *)MALLOC(col*sizeof(DL));
1.29 noro 692: for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
693: at[i] = (DL)BDY(s);
1.5 noro 694:
695: /* store coefficients separately in spmat and redmat */
696: nsp = row-nred;
697:
698: /* reducer matrix */
1.25 noro 699: /* indred : register the position of the head term */
1.24 noro 700: redmat = (CDP *)MALLOC(nred*sizeof(CDP));
1.5 noro 701: for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ )
1.30 noro 702: _tf_to_vect_compress(BDY(r),at,&redmat[i]);
1.32 noro 703:
1.5 noro 704: /* register the position of the head term */
1.31 noro 705: indred = (int *)MALLOC_ATOMIC(nred*sizeof(int));
1.5 noro 706: bzero(indred,nred*sizeof(int));
1.31 noro 707: isred = (int *)MALLOC_ATOMIC(col*sizeof(int));
1.5 noro 708: bzero(isred,col*sizeof(int));
709: for ( i = 0; i < nred; i++ ) {
710: ri = redmat[i];
1.32 noro 711: indred[i] = ri->body[0];
1.24 noro 712: isred[indred[i]] = 1;
1.5 noro 713: }
714:
715: spcol = col-nred;
716: /* head terms not in ht */
717: st = (DL *)MALLOC(spcol*sizeof(DL));
718: for ( j = 0, k = 0; j < col; j++ )
719: if ( !isred[j] )
720: st[k++] = at[j];
1.50 ! noro 721: get_eg(&tmp1); add_eg(&eg_conv,&tmp0,&tmp1);
! 722: init_eg(&eg_split_conv); add_eg(&eg_split_conv,&tmp0,&tmp1);
1.5 noro 723:
1.24 noro 724: get_eg(&tmp1);
1.5 noro 725: /* spoly matrix; stored in reduced form; terms in ht[] are omitted */
1.25 noro 726: spmat = (int **)MALLOC(nsp*sizeof(int *));
1.31 noro 727: w = (int *)MALLOC_ATOMIC(col*sizeof(int));
1.25 noro 728:
729: /* skip reducers in blist */
730: for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ );
731: for ( i = 0; r; r = NEXT(r) ) {
1.5 noro 732: bzero(w,col*sizeof(int));
733: _dpmod_to_vect(BDY(r),at,w);
1.24 noro 734: reduce_sp_by_red_mod_compress(w,redmat,indred,nred,col,m);
1.25 noro 735: for ( j = 0; j < col; j++ )
736: if ( w[j] )
737: break;
738: if ( j < col ) {
739: w1 = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
740: for ( j = 0, k = 0; j < col; j++ )
741: if ( !isred[j] )
742: w1[k++] = w[j];
743: spmat[i] = w1;
744: i++;
745: }
1.5 noro 746: }
1.25 noro 747: /* update nsp */
1.35 noro 748: nsp0 = nsp;
1.25 noro 749: nsp = i;
1.5 noro 750:
1.30 noro 751: /* XXX free redmat explicitly */
752: for ( k = 0; k < nred; k++ ) {
753: GC_free(BDY(redmat[k]));
754: GC_free(redmat[k]);
755: }
756:
1.5 noro 757: get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);
758: init_eg(&eg_split_elim1); add_eg(&eg_split_elim1,&tmp1,&tmp0);
759:
760: colstat = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
1.25 noro 761: bzero(colstat,spcol*sizeof(int));
1.5 noro 762: for ( i = 0, nonzero=0; i < nsp; i++ )
763: for ( j = 0; j < spcol; j++ )
764: if ( spmat[i][j] )
1.1 noro 765: nonzero++;
1.17 noro 766: if ( DP_Print && nsp )
1.5 noro 767: fprintf(asir_out,"spmat : %d x %d (nonzero=%f%%)...",
768: nsp,spcol,((double)nonzero*100)/(nsp*spcol));
1.17 noro 769: if ( nsp )
770: rank = generic_gauss_elim_mod(spmat,nsp,spcol,m,colstat);
771: else
772: rank = 0;
1.5 noro 773: get_eg(&tmp1); add_eg(&eg_elim2,&tmp0,&tmp1);
774: init_eg(&eg_split_elim2); add_eg(&eg_split_elim2,&tmp0,&tmp1);
775:
1.12 noro 776: if ( DP_Print ) {
1.1 noro 777: fprintf(asir_out,"done rank = %d\n",rank,row,col);
778: print_eg("Symb",&eg_split_symb);
1.50 ! noro 779: print_eg("Conv",&eg_split_conv);
1.5 noro 780: print_eg("Elim1",&eg_split_elim1);
781: print_eg("Elim2",&eg_split_elim2);
1.1 noro 782: fprintf(asir_out,"\n");
783: }
1.25 noro 784:
1.35 noro 785: NZR += rank;
786: ZR += nsp0-rank;
787:
1.25 noro 788: if ( !rank )
789: continue;
790:
1.5 noro 791: for ( j = 0, i = 0; j < spcol; j++ )
1.1 noro 792: if ( colstat[j] ) {
1.5 noro 793: mp0 = 0;
1.13 noro 794: NEXTMP(mp0,mp); mp->dl = st[j]; mp->c = STOI(1);
1.5 noro 795: for ( k = j+1; k < spcol; k++ )
796: if ( !colstat[k] && spmat[i][k] ) {
1.13 noro 797: NEXTMP(mp0,mp); mp->dl = st[k];
1.5 noro 798: mp->c = STOI(spmat[i][k]);
1.1 noro 799: }
1.5 noro 800: NEXT(mp) = 0;
801: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
802: nh = newps_mod(nf,m);
1.32 noro 803: if ( nh == pscalen ) {
804: psca = (unsigned int **)
805: REALLOC(psca,2*pscalen*sizeof(unsigned int *));
806: pscalen *= 2;
807: }
808: dptoca(ps[nh],&psca[nh]);
1.5 noro 809: d = updpairs(d,g,nh);
810: g = updbase(g,nh);
811: gall = append_one(gall,nh);
1.1 noro 812: i++;
813: }
1.30 noro 814:
815: /* XXX free spmat[] explicitly */
816: for ( j = 0; j < nsp; j++ ) {
817: GC_free(spmat[j]);
818: }
1.36 noro 819: }
820: if ( DP_Print ) {
821: print_eg("Symb",&eg_symb);
1.50 ! noro 822: print_eg("Conv",&eg_conv);
1.36 noro 823: print_eg("Elim1",&eg_elim1);
824: print_eg("Elim2",&eg_elim2);
825: fflush(asir_out);
826: }
827: return g;
828: }
829:
1.37 noro 830: NODE gb_f4_mod_old(NODE f,int m)
1.36 noro 831: {
832: int i,j,k,nh,row,col,nv;
833: NODE r,g,gall;
834: NODE s,s0;
835: DP_pairs d,dm,dr,t;
1.37 noro 836: DP nf,f2,sp,sd,sdm,tdp;
1.36 noro 837: MP mp,mp0;
1.37 noro 838: NODE blist,bt;
1.36 noro 839: DL *ht,*at,*st;
840: int **spmat,**redmat;
841: int *colstat,*w;
842: int rank,nred,nsp,nonzero,spcol;
843: int *indred,*isred,*ri;
1.37 noro 844: struct oEGT tmp0,tmp1,eg_split_symb,eg_split_elim1,eg_split_elim2;
1.36 noro 845: extern struct oEGT eg_symb,eg_elim1,eg_elim2;
846:
847: init_eg(&eg_symb); init_eg(&eg_elim1); init_eg(&eg_elim2);
848: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
849: i = (int)BDY(r);
850: d = updpairs(d,g,i);
851: g = updbase(g,i);
852: gall = append_one(gall,i);
853: }
854: if ( gall )
855: nv = ((DP)ps[(int)BDY(gall)])->nv;
856: while ( d ) {
857: get_eg(&tmp0);
858: minsugar(d,&dm,&dr); d = dr;
859: if ( DP_Print )
860: fprintf(asir_out,"sugar=%d\n",dm->sugar);
861: blist = 0; s0 = 0;
862: /* asph : sum of all head terms of spoly */
863: for ( t = dm; t; t = NEXT(t) ) {
864: _dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp);
865: if ( sp ) {
866: MKNODE(bt,sp,blist); blist = bt;
867: s0 = symb_merge(s0,dp_dllist(sp),nv);
868: }
869: }
870: /* s0 : all the terms appeared in symbolic redunction */
871: for ( s = s0, nred = 0; s; s = NEXT(s) ) {
872: for ( r = gall; r; r = NEXT(r) )
873: if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
874: break;
875: if ( r ) {
876: dltod(BDY(s),nv,&tdp);
877: dp_subd(tdp,ps[(int)BDY(r)],&sd);
878: _dp_mod(sd,m,0,&sdm);
879: mulmd_dup(m,sdm,ps[(int)BDY(r)],&f2);
880: MKNODE(bt,f2,blist); blist = bt;
881: s = symb_merge(s,dp_dllist(f2),nv);
882: nred++;
883: }
884: }
885:
886: get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
887: init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
888:
889: /* the first nred polys in blist are reducers */
890: /* row = the number of all the polys */
891: for ( r = blist, row = 0; r; r = NEXT(r), row++ );
892:
893: /* head terms of reducers */
894: ht = (DL *)MALLOC(nred*sizeof(DL));
895: for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
896: ht[i] = BDY((DP)BDY(r))->dl;
897:
898: /* col = number of all terms */
899: for ( s = s0, col = 0; s; s = NEXT(s), col++ );
900:
901: /* head terms of all terms */
902: at = (DL *)MALLOC(col*sizeof(DL));
903: for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
904: at[i] = (DL)BDY(s);
905:
906: /* store coefficients separately in spmat and redmat */
907: nsp = row-nred;
908:
909: /* reducer matrix */
910: redmat = (int **)almat(nred,col);
911: for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ )
912: _dpmod_to_vect(BDY(r),at,redmat[i]);
913: /* XXX */
914: /* reduce_reducers_mod(redmat,nred,col,m); */
915: /* register the position of the head term */
916: indred = (int *)MALLOC(nred*sizeof(int));
917: bzero(indred,nred*sizeof(int));
918: isred = (int *)MALLOC(col*sizeof(int));
919: bzero(isred,col*sizeof(int));
920: for ( i = 0; i < nred; i++ ) {
921: ri = redmat[i];
922: for ( j = 0; j < col && !ri[j]; j++ );
923: indred[i] = j;
924: isred[j] = 1;
925: }
926:
927: spcol = col-nred;
928: /* head terms not in ht */
929: st = (DL *)MALLOC(spcol*sizeof(DL));
930: for ( j = 0, k = 0; j < col; j++ )
931: if ( !isred[j] )
932: st[k++] = at[j];
933:
934: /* spoly matrix; stored in reduced form; terms in ht[] are omitted */
935: spmat = almat(nsp,spcol);
936: w = (int *)MALLOC(col*sizeof(int));
937: for ( ; i < row; r = NEXT(r), i++ ) {
938: bzero(w,col*sizeof(int));
939: _dpmod_to_vect(BDY(r),at,w);
940: reduce_sp_by_red_mod(w,redmat,indred,nred,col,m);
941: for ( j = 0, k = 0; j < col; j++ )
942: if ( !isred[j] )
943: spmat[i-nred][k++] = w[j];
944: }
945:
946: get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);
947: init_eg(&eg_split_elim1); add_eg(&eg_split_elim1,&tmp1,&tmp0);
948:
949: colstat = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
950: for ( i = 0, nonzero=0; i < nsp; i++ )
951: for ( j = 0; j < spcol; j++ )
952: if ( spmat[i][j] )
953: nonzero++;
954: if ( DP_Print && nsp )
955: fprintf(asir_out,"spmat : %d x %d (nonzero=%f%%)...",
956: nsp,spcol,((double)nonzero*100)/(nsp*spcol));
957: if ( nsp )
958: rank = generic_gauss_elim_mod(spmat,nsp,spcol,m,colstat);
959: else
960: rank = 0;
961: get_eg(&tmp1); add_eg(&eg_elim2,&tmp0,&tmp1);
962: init_eg(&eg_split_elim2); add_eg(&eg_split_elim2,&tmp0,&tmp1);
963:
964: if ( DP_Print ) {
965: fprintf(asir_out,"done rank = %d\n",rank,row,col);
966: print_eg("Symb",&eg_split_symb);
967: print_eg("Elim1",&eg_split_elim1);
968: print_eg("Elim2",&eg_split_elim2);
969: fprintf(asir_out,"\n");
970: }
971: for ( j = 0, i = 0; j < spcol; j++ )
972: if ( colstat[j] ) {
973: mp0 = 0;
974: NEXTMP(mp0,mp); mp->dl = st[j]; mp->c = STOI(1);
975: for ( k = j+1; k < spcol; k++ )
976: if ( !colstat[k] && spmat[i][k] ) {
977: NEXTMP(mp0,mp); mp->dl = st[k];
978: mp->c = STOI(spmat[i][k]);
979: }
980: NEXT(mp) = 0;
981: MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
982: nh = newps_mod(nf,m);
983: d = updpairs(d,g,nh);
984: g = updbase(g,nh);
985: gall = append_one(gall,nh);
986: i++;
987: }
1.1 noro 988: }
1.12 noro 989: if ( DP_Print ) {
1.1 noro 990: print_eg("Symb",&eg_symb);
1.5 noro 991: print_eg("Elim1",&eg_elim1);
992: print_eg("Elim2",&eg_elim2);
1.1 noro 993: fflush(asir_out);
994: }
995: return g;
996: }
997:
1.37 noro 998: int DPPlength(DP_pairs n)
1.1 noro 999: {
1000: int i;
1001:
1002: for ( i = 0; n; n = NEXT(n), i++ );
1003: return i;
1004: }
1005:
1.37 noro 1006: void printdl(DL dl)
1.1 noro 1007: {
1008: int i;
1009:
1010: fprintf(asir_out,"<<");
1011: for ( i = 0; i < CNVars-1; i++ )
1012: fprintf(asir_out,"%d,",dl->d[i]);
1013: fprintf(asir_out,"%d>>",dl->d[i]);
1014: }
1015:
1.37 noro 1016: void pltovl(LIST l,VL *vl)
1.1 noro 1017: {
1018: NODE n;
1019: VL r,r0;
1020:
1021: n = BDY(l);
1022: for ( r0 = 0; n; n = NEXT(n) ) {
1023: NEXTVL(r0,r); r->v = VR((P)BDY(n));
1024: }
1025: if ( r0 ) NEXT(r) = 0;
1026: *vl = r0;
1027: }
1028:
1.37 noro 1029: void makesubst(VL v,NODE *s)
1.1 noro 1030: {
1031: NODE r,r0;
1032: Q q;
1033: unsigned int n;
1034:
1035: for ( r0 = 0; v; v = NEXT(v) ) {
1036: NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
1037: #if defined(_PA_RISC1_1)
1038: n = mrand48()&BMASK; UTOQ(n,q);
1039: #else
1040: n = random(); UTOQ(n,q);
1041: #endif
1042: NEXTNODE(r0,r); BDY(r) = (pointer)q;
1043: }
1044: if ( r0 ) NEXT(r) = 0;
1045: *s = r0;
1046: }
1047:
1.37 noro 1048: void printsubst(NODE s)
1.1 noro 1049: {
1050: fputc('[',asir_out);
1051: while ( s ) {
1052: printv(CO,(V)BDY(s)); s = NEXT(s);
1053: fprintf(asir_out,"->%d",QTOS((Q)BDY(s)));
1054: if ( NEXT(s) ) {
1055: fputc(',',asir_out); s = NEXT(s);
1056: } else
1057: break;
1058: }
1059: fprintf(asir_out,"]\n"); return;
1060: }
1061:
1.37 noro 1062: void vlminus(VL v,VL w,VL *d)
1.1 noro 1063: {
1064: int i,j,n,m;
1065: V *va,*wa;
1066: V a;
1067: VL r,r0;
1068: VL t;
1069:
1070: for ( n = 0, t = v; t; t = NEXT(t), n++ );
1071: va = (V *)ALLOCA(n*sizeof(V));
1072: for ( i = 0, t = v; t; t = NEXT(t), i++ )
1073: va[i] = t->v;
1074: for ( m = 0, t = w; t; t = NEXT(t), m++ );
1075: wa = (V *)ALLOCA(m*sizeof(V));
1076: for ( i = 0, t = w; t; t = NEXT(t), i++ )
1077: wa[i] = t->v;
1078: for ( i = 0; i < n; i++ ) {
1079: a = va[i];
1080: for ( j = 0; j < m; j++ )
1081: if ( a == wa[j] )
1082: break;
1083: if ( j < m )
1084: va[i] = 0;
1085: }
1086: for ( r0 = 0, i = 0; i < n; i++ )
1087: if ( va[i] ) { NEXTVL(r0,r); r->v = va[i]; }
1088: if ( r0 ) NEXT(r) = 0;
1089: *d = r0;
1090: }
1091:
1.37 noro 1092: int validhc(P a,int m,NODE s)
1.1 noro 1093: {
1094: P c,c1;
1095: V v;
1096:
1097: if ( !a )
1098: return 0;
1099: for ( c = a; s; s = NEXT(s) ) {
1100: v = (V)BDY(s); s = NEXT(s);
1101: substp(CO,c,v,(P)BDY(s),&c1); c = c1;
1102: }
1103: ptomp(m,c,&c1);
1104: return c1 ? 1 : 0;
1105: }
1106:
1.37 noro 1107: void setup_arrays(NODE f,int m,NODE *r)
1.1 noro 1108: {
1109: int i;
1.7 noro 1110: NODE s,s0,f0;
1.1 noro 1111:
1.17 noro 1112: #if 1
1.7 noro 1113: f0 = f = NODE_sortb(f,1);
1.17 noro 1114: #else
1115: f0 = f;
1116: #endif
1.1 noro 1117: psn = length(f); pslen = 2*psn;
1118: ps = (DP *)MALLOC(pslen*sizeof(DP));
1119: psh = (DL *)MALLOC(pslen*sizeof(DL));
1120: pss = (int *)MALLOC(pslen*sizeof(int));
1121: psc = (P *)MALLOC(pslen*sizeof(P));
1122: for ( i = 0; i < psn; i++, f = NEXT(f) ) {
1123: prim_part((DP)BDY(f),m,&ps[i]);
1124: if ( Demand )
1.7 noro 1125: dp_save(i,(Obj)ps[i],0);
1.1 noro 1126: psh[i] = BDY(ps[i])->dl;
1127: pss[i] = ps[i]->sugar;
1128: psc[i] = BDY(ps[i])->c;
1129: }
1.23 noro 1130: if ( GenTrace ) {
1.7 noro 1131: Q q;
1132: STRING fname;
1133: LIST input;
1.23 noro 1134: NODE arg,t,t1;
1.37 noro 1135: Obj obj;
1.23 noro 1136:
1137: t = 0;
1138: for ( i = psn-1; i >= 0; i-- ) {
1139: MKNODE(t1,ps[i],t);
1140: t = t1;
1141: }
1142: MKLIST(input,t);
1.7 noro 1143:
1.23 noro 1144: if ( OXCheck >= 0 ) {
1145: STOQ(OXCheck,q);
1146: MKSTR(fname,"register_input");
1147: arg = mknode(3,q,fname,input);
1.37 noro 1148: Pox_cmo_rpc(arg,&obj);
1.23 noro 1149: } else if ( OXCheck < 0 ) {
1150: MKNODE(AllTraceList,input,0);
1151: }
1.7 noro 1152: }
1.1 noro 1153: for ( s0 = 0, i = 0; i < psn; i++ ) {
1154: NEXTNODE(s0,s); BDY(s) = (pointer)i;
1155: }
1156: if ( s0 ) NEXT(s) = 0;
1157: *r = s0;
1158: }
1159:
1.37 noro 1160: void prim_part(DP f,int m,DP *r)
1.1 noro 1161: {
1.7 noro 1162: P d,t;
1163:
1.1 noro 1164: if ( m > 0 ) {
1165: if ( PCoeffs )
1166: dp_prim_mod(f,m,r);
1167: else
1.49 noro 1168: _dp_monic(f,m,r);
1.1 noro 1169: } else {
1.44 noro 1170: if ( dp_fcoeffs || PCoeffs )
1.1 noro 1171: dp_prim(f,r);
1172: else
1173: dp_ptozp(f,r);
1.7 noro 1174: if ( GenTrace && TraceList ) {
1.23 noro 1175: /* adust the denominator according to the final
1176: content reduction */
1.7 noro 1177: divsp(CO,BDY(f)->c,BDY(*r)->c,&d);
1178: mulp(CO,(P)ARG3(BDY((LIST)BDY(TraceList))),d,&t);
1179: ARG3(BDY((LIST)BDY(TraceList))) = t;
1180: }
1.1 noro 1181: }
1182: }
1183:
1.37 noro 1184: NODE /* of DP */ NODE_sortb_insert( DP newdp, NODE /* of DP */ nd, int dec )
1.1 noro 1185: {
1186: register NODE last, p;
1187: register DL newdl = BDY(newdp)->dl;
1188: register int (*cmpfun)() = cmpdl, nv = CNVars;
1189: NODE newnd;
1190: int sgn = dec ? 1 : -1;
1191: MKNODE( newnd, newdp, 0 );
1192: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(last))->dl ) > 0 ) {
1193: NEXT(newnd) = last;
1194: return newnd;
1195: }
1196: for ( ; p = NEXT(last); last = p )
1197: if ( sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(p))->dl ) > 0 ) break;
1198: if ( p ) NEXT(NEXT(last) = newnd) = p;
1199: else NEXT(last) = newnd;
1200: return nd;
1201: }
1202:
1.37 noro 1203: NODE NODE_sortb( NODE node, int dec )
1.1 noro 1204: {
1205: register NODE nd, ans;
1206:
1207: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
1208: ans = NODE_sortb_insert( (DP) BDY(nd), ans, dec );
1209: return ans;
1210: }
1211:
1.37 noro 1212: NODE /* of index */ NODE_sortbi_insert( int newdpi, NODE /* of index */ nd, int dec )
1.1 noro 1213: {
1214: register NODE last, p;
1215: register DL newdl = psh[newdpi];
1216: register int (*cmpfun)() = cmpdl, nv = CNVars;
1217: NODE newnd;
1218: int sgn = dec ? 1 : -1;
1219: MKNODE( newnd, newdpi, 0 );
1220: if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(last)] ) > 0 ) {
1221: NEXT(newnd) = last;
1222: return newnd;
1223: }
1224: for ( ; p = NEXT(last); last = p )
1225: if ( sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(p)] ) > 0 ) break;
1226: if ( p ) NEXT(NEXT(last) = newnd) = p;
1227: else NEXT(last) = newnd;
1228: return nd;
1229: }
1230:
1.37 noro 1231: NODE NODE_sortbi( NODE node, int dec )
1.1 noro 1232: {
1233: register NODE nd, ans;
1234:
1235: for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
1236: ans = NODE_sortbi_insert( (int) BDY(nd), ans, dec );
1237: return ans;
1238: }
1239:
1.37 noro 1240: void reduceall(NODE in,NODE *h)
1.1 noro 1241: {
1242: NODE r,t,top;
1243: int n,i,j;
1244: int *w;
1245: DP g,g1;
1246: struct oEGT tmp0,tmp1;
1247:
1248: if ( NoRA ) {
1249: *h = in; return;
1250: }
1.12 noro 1251: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1252: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1253: }
1254: r = NODE_sortbi(in,0);
1255: n = length(r);
1256: w = (int *)ALLOCA(n*sizeof(int));
1257: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1258: w[i] = (int)BDY(t);
1259: for ( i = 0; i < n; i++ ) {
1260: for ( top = 0, j = n-1; j >= 0; j-- )
1261: if ( j != i ) {
1262: MKNODE(t,(pointer)w[j],top); top = t;
1263: }
1264: get_eg(&tmp0);
1265: dp_load(w[i],&ps[w[i]]);
1266:
1.7 noro 1267: if ( GenTrace ) {
1268: Q q;
1269: NODE node;
1270: LIST hist;
1271:
1272: STOQ(w[i],q);
1273: node = mknode(4,ONE,q,ONE,ONE);
1274: MKLIST(hist,node);
1275: MKNODE(TraceList,hist,0);
1276: }
1.16 noro 1277: _dp_nf(top,ps[w[i]],ps,1,&g);
1.1 noro 1278: prim_part(g,0,&g1);
1279: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12 noro 1280: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1281: fprintf(asir_out,"."); fflush(asir_out);
1282: }
1283: w[i] = newps(g1,0,(NODE)0);
1284: }
1285: for ( top = 0, j = n-1; j >= 0; j-- ) {
1286: MKNODE(t,(pointer)w[j],top); top = t;
1287: }
1288: *h = top;
1.12 noro 1289: if ( DP_Print || DP_PrintShort )
1.1 noro 1290: fprintf(asir_out,"\n");
1291: }
1292:
1.37 noro 1293: void reduceall_mod(NODE in,int m,NODE *h)
1.1 noro 1294: {
1295: NODE r,t,top;
1296: int n,i,j;
1297: int *w;
1.12 noro 1298: DP g,p;
1.1 noro 1299: struct oEGT tmp0,tmp1;
1300:
1301: if ( NoRA ) {
1302: *h = in; return;
1303: }
1.12 noro 1304: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1305: fprintf(asir_out,"reduceall\n"); fflush(asir_out);
1306: }
1307: r = NODE_sortbi(in,0);
1308: n = length(r);
1309: w = (int *)ALLOCA(n*sizeof(int));
1310: for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
1311: w[i] = (int)BDY(t);
1312: for ( i = 0; i < n; i++ ) {
1313: for ( top = 0, j = n-1; j >= 0; j-- )
1314: if ( j != i ) {
1315: MKNODE(t,(pointer)w[j],top); top = t;
1316: }
1317: get_eg(&tmp0);
1318: if ( PCoeffs )
1319: dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1.12 noro 1320: else {
1321: dpto_dp(ps[w[i]],&p);
1322: _dp_nf_mod_destructive(top,p,ps,m,1,&g);
1323: }
1.1 noro 1324: get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12 noro 1325: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1326: fprintf(asir_out,"."); fflush(asir_out);
1327: }
1328: w[i] = newps_mod(g,m);
1329: }
1330: for ( top = 0, j = n-1; j >= 0; j-- ) {
1331: MKNODE(t,(pointer)w[j],top); top = t;
1332: }
1333: *h = top;
1.12 noro 1334: if ( DP_Print || DP_PrintShort )
1.1 noro 1335: fprintf(asir_out,"\n");
1336: }
1337:
1.37 noro 1338: int newps(DP a,int m,NODE subst)
1.1 noro 1339: {
1340: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1341: return -1;
1342: if ( psn == pslen ) {
1343: pslen *= 2;
1344: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1345: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1346: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1347: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1348: if ( m )
1349: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1350: }
1351: if ( Demand ) {
1352: if ( doing_f4 )
1353: ps[psn] = a;
1354: else
1355: ps[psn] = 0;
1.7 noro 1356: dp_save(psn,(Obj)a,0);
1.1 noro 1357: } else
1358: ps[psn] = a;
1359: psh[psn] = BDY(a)->dl;
1360: pss[psn] = a->sugar;
1361: psc[psn] = BDY(a)->c;
1362: if ( m )
1363: _dp_mod(a,m,subst,&psm[psn]);
1.7 noro 1364: if ( GenTrace ) {
1365: NODE tn,tr,tr1;
1.23 noro 1366: LIST trace,trace1;
1367: NODE arg;
1368: Q q1,q2;
1369: STRING fname;
1.37 noro 1370: Obj obj;
1.7 noro 1371:
1372: /* reverse the TraceList */
1373: tn = TraceList;
1374: for ( tr = 0; tn; tn = NEXT(tn) ) {
1375: MKNODE(tr1,BDY(tn),tr); tr = tr1;
1376: }
1377: MKLIST(trace,tr);
1378: if ( OXCheck >= 0 ) {
1379: STOQ(OXCheck,q1);
1380: MKSTR(fname,"check_trace");
1381: STOQ(psn,q2);
1382: arg = mknode(5,q1,fname,a,q2,trace);
1.37 noro 1383: Pox_cmo_rpc(arg,&obj);
1.23 noro 1384: } else if ( OXCheck < 0 ) {
1385: STOQ(psn,q1);
1386: tn = mknode(2,q1,trace);
1387: MKLIST(trace1,tn);
1388: MKNODE(tr,trace1,AllTraceList);
1389: AllTraceList = tr;
1.7 noro 1390: } else
1391: dp_save(psn,(Obj)trace,"t");
1392: TraceList = 0;
1393: }
1.1 noro 1394: return psn++;
1395: }
1396:
1.37 noro 1397: int newps_nosave(DP a,int m,NODE subst)
1.1 noro 1398: {
1399: if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
1400: return -1;
1401: if ( psn == pslen ) {
1402: pslen *= 2;
1403: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1404: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1405: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1406: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
1407: if ( m )
1408: psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
1409: }
1410: ps[psn] = 0;
1411: psh[psn] = BDY(a)->dl;
1412: pss[psn] = a->sugar;
1413: psc[psn] = BDY(a)->c;
1414: if ( m )
1415: _dp_mod(a,m,subst,&psm[psn]);
1416: return psn++;
1417: }
1418:
1.37 noro 1419: int newps_mod(DP a,int m)
1.1 noro 1420: {
1421: if ( psn == pslen ) {
1422: pslen *= 2;
1423: ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
1424: psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
1425: pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
1426: psc = (P *)REALLOC((char *)psc,pslen*sizeof(P)); /* XXX */
1427: }
1428: ps[psn] = a;
1429: psh[psn] = BDY(ps[psn])->dl;
1430: pss[psn] = ps[psn]->sugar;
1431: return psn++;
1432: }
1433:
1.37 noro 1434: void reducebase_dehomo(NODE f,NODE *g)
1.1 noro 1435: {
1436: int n,i,j,k;
1437: int *r;
1438: DL *w,d;
1439: DP u;
1440: NODE t,top;
1441:
1442: n = length(f);
1443: w = (DL *)ALLOCA(n*sizeof(DL));
1444: r = (int *)ALLOCA(n*sizeof(int));
1445: for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {
1446: r[i] = (int)BDY(t); w[i] = psh[r[i]];
1447: }
1448: for ( i = 0; i < n; i++ ) {
1449: for ( j = 0, d = w[i]; j < n; j++ ) {
1450: if ( j != i ) {
1451: for ( k = 0; k < NVars; k++ )
1452: if ( d->d[k] < w[j]->d[k] )
1453: break;
1454: if ( k == NVars )
1455: break;
1456: }
1457: }
1458: if ( j != n )
1459: r[i] = -1;
1460: }
1461: for ( top = 0, i = n-1; i >= 0; i-- )
1462: if ( r[i] >= 0 ) {
1.7 noro 1463: dp_load(r[i],&ps[r[i]]); dp_dehomo(ps[r[i]],&u);
1464: if ( GenTrace ) {
1465: Q q;
1466: LIST hist;
1467: NODE node;
1468:
1469: STOQ(r[i],q);
1470: node = mknode(4,0,q,0,0);
1471: MKLIST(hist,node);
1472: MKNODE(TraceList,hist,0);
1473: }
1474: j = newps(u,0,0);
1.1 noro 1475: MKNODE(t,j,top); top = t;
1476: }
1477: *g = top;
1478: }
1479:
1.37 noro 1480: NODE append_one(NODE f,int n)
1.1 noro 1481: {
1482: NODE t;
1483:
1484: if ( Reverse || !f ) {
1485: MKNODE(t,(pointer)n,f); return t;
1486: } else {
1487: for ( t = f; NEXT(t); t = NEXT(t) );
1488: MKNODE(NEXT(t),(pointer)n,0);
1489: return f;
1490: }
1491: }
1492:
1.37 noro 1493: DP_pairs minp( DP_pairs d, DP_pairs *prest )
1.1 noro 1494: {
1495: register DP_pairs m, ml, p, l;
1496: register DL lcm;
1497: register int s, nv = CNVars;
1498: register int (*cmpfun)() = cmpdl;
1499:
1500: if ( !(p = NEXT(m = d)) ) {
1501: *prest = p;
1502: NEXT(m) = 0;
1503: return m;
1504: }
1505: for ( lcm = m->lcm, s = m->sugar, ml = 0, l = m; p; p = NEXT(l = p) )
1506: if ( NoSugar ? (*cmpfun)( nv, lcm, p->lcm ) >= 0 :
1507: (s > p->sugar || s == p->sugar && (*cmpfun)( nv, lcm, p->lcm ) >= 0) )
1508: ml = l, lcm = (m = p)->lcm, s = p->sugar;
1509: if ( !ml ) *prest = NEXT(m);
1510: else {
1511: NEXT(ml) = NEXT(m);
1512: *prest = d;
1513: }
1514: NEXT(m) = 0;
1515: return m;
1516: }
1517:
1.37 noro 1518: void minsugar(DP_pairs d,DP_pairs *dm,DP_pairs *dr)
1.1 noro 1519: {
1520: int msugar;
1521: DP_pairs t,dm0,dr0,dmt,drt;
1522:
1523: for ( msugar = d->sugar, t = NEXT(d); t; t = NEXT(t) )
1524: if ( t->sugar < msugar )
1525: msugar = t->sugar;
1526: dm0 = 0; dr0 = 0;
1527: for ( t = d; t; t = NEXT(t) ) {
1528: if ( t->sugar == msugar ) {
1529: NEXTDPP(dm0,dmt);
1530: dmt->dp1 = t->dp1; dmt->dp2 = t->dp2;
1531: dmt->lcm = t->lcm; dmt->sugar = t->sugar;
1532: } else {
1533: NEXTDPP(dr0,drt);
1534: drt->dp1 = t->dp1; drt->dp2 = t->dp2;
1535: drt->lcm = t->lcm; drt->sugar = t->sugar;
1536: }
1537: }
1538: if ( dm0 ) NEXT(dmt) = 0;
1539: if ( dr0 ) NEXT(drt) = 0;
1540: *dm = dm0; *dr = dr0;
1541: }
1542:
1.37 noro 1543: NODE gb(NODE f,int m,NODE subst)
1.1 noro 1544: {
1.42 noro 1545: int i,nh,prev,mag,mag0,magt;
1.1 noro 1546: NODE r,g,gall;
1.37 noro 1547: DP_pairs d;
1.1 noro 1548: DP_pairs l;
1549: DP h,nf,nfm,dp1,dp2;
1550: MP mp;
1.37 noro 1551: struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tnp0,tnp1;
1.1 noro 1552: int skip_nf_flag;
1553: double t_0;
1.7 noro 1554: Q q;
1.6 noro 1555: int new_sugar;
1.1 noro 1556: static prev_sugar = -1;
1557:
1558: Max_mag = 0;
1.42 noro 1559: Max_coef = 0;
1.1 noro 1560: prev = 1;
1561: doing_f4 = 0;
1562: if ( m ) {
1563: psm = (DP *)MALLOC(pslen*sizeof(DP));
1564: for ( i = 0; i < psn; i++ )
1565: if ( psh[i] && !validhc(psc[i],m,subst) )
1566: return 0;
1567: else
1568: _dp_mod(ps[i],m,subst,&psm[i]);
1569: }
1570: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1571: i = (int)BDY(r);
1572: d = updpairs(d,g,i);
1573: g = updbase(g,i);
1574: gall = append_one(gall,i);
1575: }
1576: while ( d ) {
1.14 noro 1577: l = minp(d,&d);
1.1 noro 1578: if ( m ) {
1.6 noro 1579: _dp_sp_mod_dup(psm[l->dp1],psm[l->dp2],m,&h);
1.10 noro 1580: if ( h )
1581: new_sugar = h->sugar;
1.1 noro 1582: get_eg(&tnfm0);
1.6 noro 1583: _dp_nf_mod_destructive(gall,h,psm,m,0,&nfm);
1.1 noro 1584: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1585: } else
1586: nfm = (DP)1;
1587: if ( nfm ) {
1588: if ( Demand ) {
1589: if ( dp_load_t(psn,&nf) ) {
1590: skip_nf_flag = 1;
1591: goto skip_nf;
1592: } else {
1593: skip_nf_flag = 0;
1594: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1595: dp_sp(dp1,dp2,&h);
1596: }
1597: } else
1598: dp_sp(ps[l->dp1],ps[l->dp2],&h);
1.7 noro 1599: if ( GenTrace ) {
1600: STOQ(l->dp1,q); ARG1(BDY((LIST)BDY(NEXT(TraceList)))) = q;
1601: STOQ(l->dp2,q); ARG1(BDY((LIST)BDY(TraceList))) = q;
1602: }
1.10 noro 1603: if ( h )
1604: new_sugar = h->sugar;
1.1 noro 1605: get_eg(&tnf0);
1606: t_0 = get_rtime();
1.18 noro 1607: if ( PCoeffs || dp_fcoeffs )
1.16 noro 1608: _dp_nf(gall,h,ps,!Top,&nf);
1609: else
1.21 noro 1610: _dp_nf_z(gall,h,ps,!Top,DP_Multiple,&nf);
1.50 ! noro 1611: if ( DP_Print && nf )
1.1 noro 1612: fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
1613: get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
1614: } else
1615: nf = 0;
1616: skip_nf:
1617: if ( nf ) {
1618: NZR++;
1619: get_eg(&tpz0);
1620: prim_part(nf,0,&h);
1621: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1622: get_eg(&tnp0);
1623: if ( Demand && skip_nf_flag )
1624: nh = newps_nosave(h,m,subst);
1625: else
1626: nh = newps(h,m,subst);
1627: get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
1628: if ( nh < 0 )
1629: return 0;
1630: d = updpairs(d,g,nh);
1631: g = updbase(g,nh);
1632: gall = append_one(gall,nh);
1633: if ( !dp_fcoeffs && ShowMag ) {
1.42 noro 1634: for ( mag = 0, mag0 = 0, mp = BDY(h); mp; mp = NEXT(mp) ) {
1635: magt = p_mag((P)mp->c);
1636: mag0 = MAX(mag0,magt);
1637: mag += magt;
1638: }
1639: Max_coef = MAX(Max_coef,mag0);
1.1 noro 1640: Max_mag = MAX(Max_mag,mag);
1641: }
1.12 noro 1642: if ( DP_Print ) {
1.1 noro 1643: if ( !prev )
1644: fprintf(asir_out,"\n");
1645: print_split_e(&tnf0,&tnf1); print_split_e(&tpz0,&tpz1);
1646: printdl(psh[nh]);
1647: fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",
1648: l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
1649: pss[nh]);
1650: if ( ShowMag )
1.42 noro 1651: fprintf(asir_out,",mag=(%d,%d)",mag,mag0);
1.1 noro 1652: fprintf(asir_out,"\n"); fflush(asir_out);
1.12 noro 1653: } else if ( DP_PrintShort ) {
1.1 noro 1654: fprintf(asir_out,"+"); fflush(asir_out);
1655: }
1656: prev = 1;
1657: } else {
1658: if ( m )
1659: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1660: ZR++;
1.12 noro 1661: if ( DP_Print || DP_PrintShort ) {
1.6 noro 1662: if ( new_sugar != prev_sugar ) {
1663: fprintf(asir_out,"[%d]",new_sugar);
1664: prev_sugar = new_sugar;
1.1 noro 1665: }
1666: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1667: }
1668: }
1669: }
1.12 noro 1670: if ( DP_Print || DP_PrintShort )
1.1 noro 1671: fprintf(asir_out,"gb done\n");
1672: return g;
1673: }
1674:
1.37 noro 1675: NODE gb_mod(NODE f,int m)
1.1 noro 1676: {
1677: int i,nh,prev;
1678: NODE r,g,gall;
1.37 noro 1679: DP_pairs d;
1.1 noro 1680: DP_pairs l;
1681: DP h,nf;
1.37 noro 1682: struct oEGT tnfm0,tnfm1,tpz0,tpz1;
1.1 noro 1683:
1684: prev = 1;
1685: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1686: i = (int)BDY(r);
1687: d = updpairs(d,g,i);
1688: g = updbase(g,i);
1689: gall = append_one(gall,i);
1690: }
1691: while ( d ) {
1.14 noro 1692: l = minp(d,&d);
1.1 noro 1693: if ( PCoeffs ) {
1694: dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1.14 noro 1695: get_eg(&tnfm0);
1.1 noro 1696: dp_nf_mod(gall,h,ps,m,!Top,&nf);
1697: } else {
1.6 noro 1698: _dp_sp_mod_dup(ps[l->dp1],ps[l->dp2],m,&h);
1.14 noro 1699: get_eg(&tnfm0);
1.6 noro 1700: _dp_nf_mod_destructive(gall,h,ps,m,!Top,&nf);
1.1 noro 1701: }
1702: get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
1703: if ( nf ) {
1704: NZR++;
1705: get_eg(&tpz0);
1706: prim_part(nf,m,&h);
1707: get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
1708: nh = newps_mod(h,m);
1709: if ( nh < 0 )
1710: return 0;
1711: d = updpairs(d,g,nh);
1712: g = updbase(g,nh);
1713: gall = append_one(gall,nh);
1.12 noro 1714: if ( DP_Print ) {
1.1 noro 1715: if ( !prev )
1716: fprintf(asir_out,"\n");
1717: print_split_eg(&tnfm0,&tnfm1); fflush(asir_out);
1718: 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]);
1719: printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1.12 noro 1720: } else if ( DP_PrintShort ) {
1.6 noro 1721: fprintf(asir_out,"+"); fflush(asir_out);
1.1 noro 1722: }
1723: prev = 1;
1724: } else {
1725: add_eg(&eg_znfm,&tnfm0,&tnfm1);
1726: ZR++;
1.12 noro 1727: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1728: fprintf(asir_out,"."); fflush(asir_out); prev = 0;
1729: }
1730: }
1731: }
1.12 noro 1732: if ( DP_Print || DP_PrintShort )
1.1 noro 1733: fprintf(asir_out,"gb_mod done\n");
1734: return g;
1735: }
1736:
1.37 noro 1737: DP_pairs updpairs( DP_pairs d, NODE /* of index */ g, int t)
1.1 noro 1738: {
1.14 noro 1739: register DP_pairs d1, dd, nd;
1740: int dl,dl1;
1.1 noro 1741:
1.14 noro 1742: if ( !g ) return d;
1743: if ( !NoCriB && d ) {
1744: dl = DPPlength(d);
1745: d = criterion_B( d, t );
1746: dl -= DPPlength(d); NBP += dl;
1.1 noro 1747: }
1.14 noro 1748: d1 = newpairs( g, t );
1749: if ( NEXT(d1) ) {
1750: dl = DPPlength(d1); TP += dl;
1751: d1 = criterion_M( d1 );
1752: dl1 = DPPlength(d1); NMP += (dl-dl1); dl = dl1;
1753: d1 = criterion_F( d1 );
1754: dl1 = DPPlength(d1); NFP += (dl-dl1); dl = dl1;
1755: } else
1756: dl = 1;
1757: if ( !do_weyl )
1758: for ( dd = 0; d1; d1 = nd ) {
1759: nd = NEXT(d1);
1760: if ( !criterion_2( d1->dp1, d1->dp2 ) ) {
1761: NEXT(d1) = dd;
1762: dd = d1;
1.1 noro 1763: }
1764: }
1.14 noro 1765: else
1766: dd = d1;
1767: dl1 = DPPlength(dd); NDP += (dl-dl1);
1768: if ( !(nd = d) ) return dd;
1769: while ( nd = NEXT(d1 = nd) ) ;
1770: NEXT(d1) = dd;
1771: return d;
1.1 noro 1772: }
1773:
1.37 noro 1774: DP_pairs newpairs( NODE /* of index */ g, int t )
1.14 noro 1775: {
1776: register NODE r;
1777: register DL tdl = psh[t];
1778: register int ts;
1779: register DP_pairs p, last;
1780: int dp;
1781: register DL dl;
1782: register int s;
1.1 noro 1783:
1.14 noro 1784: ts = pss[t] - tdl->td;
1785: for ( last = 0, r = g; r; r = NEXT(r) ) {
1786: NEXT(p = NEWDP_pairs) = last;
1787: last = p;
1788: dp = p->dp1 = (int)BDY(r); p->dp2 = t;
1789: p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
1790: #if 0
1791: if ( do_weyl )
1.41 noro 1792: p->sugar = dl_weyl_weight(p->lcm);
1.14 noro 1793: else
1794: #endif
1795: p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
1796: }
1797: return last;
1798: }
1.1 noro 1799:
1.37 noro 1800: DP_pairs criterion_B( DP_pairs d, int s )
1.1 noro 1801: {
1.14 noro 1802: register DP_pairs dd, p;
1803: register DL tij, t = psh[s], dltmp;
1.1 noro 1804:
1.14 noro 1805: if ( !d ) return 0;
1806: NEWDL( dltmp, CNVars );
1807: for ( dd = 0; d; d = p ) {
1808: p = NEXT(d),
1809: tij = d->lcm;
1810: if ( tij->td != lcm_of_DL(CNVars, tij, t, dltmp )->td
1811: || !dl_equal(CNVars, tij, dltmp )
1812: || (tij->td == lcm_of_DL(CNVars, psh[d->dp1], t, dltmp )->td
1813: && dl_equal(CNVars, dltmp, tij ))
1814: || (tij->td == lcm_of_DL(CNVars, psh[d->dp2], t, dltmp )->td
1815: && dl_equal(CNVars, dltmp, tij )) ) {
1816: NEXT(d) = dd;
1817: dd = d;
1818: }
1.1 noro 1819: }
1.14 noro 1820: return dd;
1821: }
1.1 noro 1822:
1.37 noro 1823: DP_pairs criterion_M( DP_pairs d1 )
1.14 noro 1824: {
1825: register DP_pairs dd, e, d3, d2, p;
1826: register DL itdl, jtdl;
1827: register int itdltd, jtdltd;
1.1 noro 1828:
1.14 noro 1829: for ( dd = 0, e = d1; e; e = d3 ) {
1830: if ( !(d2 = NEXT(e)) ) {
1831: NEXT(e) = dd;
1832: return e;
1833: }
1834: itdltd = (itdl = e->lcm)->td;
1835: for ( d3 = 0; d2; d2 = p ) {
1836: p = NEXT(d2),
1837: jtdltd = (jtdl = d2->lcm)->td;
1838: if ( jtdltd == itdltd )
1839: if ( dl_equal(CNVars, itdl, jtdl ) ) ;
1840: else if ( dl_redble( jtdl, itdl ) ) continue;
1841: else if ( dl_redble( itdl, jtdl ) ) goto delit;
1842: else ;
1843: else if ( jtdltd > itdltd )
1844: if ( dl_redble( jtdl, itdl ) ) continue;
1.1 noro 1845: else ;
1846: else if ( dl_redble( itdl, jtdl ) ) goto delit;
1847: NEXT(d2) = d3;
1848: d3 = d2;
1849: }
1850: NEXT(e) = dd;
1851: dd = e;
1852: continue;
1853: /**/
1854: delit: NEXT(d2) = d3;
1855: d3 = d2;
1856: for ( ; p; p = d2 ) {
1857: d2 = NEXT(p);
1858: NEXT(p) = d3;
1859: d3 = p;
1860: }
1861: }
1862: return dd;
1863: }
1864:
1.37 noro 1865: static DP_pairs collect_pairs_of_hdlcm( DP_pairs d1, DP_pairs *prest )
1.1 noro 1866: {
1867: register DP_pairs w, p, r, s;
1868: register DL ti;
1869: register int td;
1870:
1871: td = (ti = (w = d1)->lcm)->td;
1872: s = NEXT(w);
1873: NEXT(w) = 0;
1874: for ( r = 0; s; s = p ) {
1875: p = NEXT(s);
1876: if ( td == s->lcm->td && dl_equal(CNVars, ti, s->lcm ) )
1877: {
1878: NEXT(s) = w;
1879: w = s;
1880: } else {
1881: NEXT(s) = r;
1882: r = s;
1883: }
1884: }
1885: *prest = r;
1886: return w;
1887: }
1888:
1.37 noro 1889: int criterion_2( int dp1, int dp2 )
1.1 noro 1890: {
1891: register int i, *d1, *d2;
1892:
1893: d1 = psh[dp1]->d, d2 = psh[dp2]->d;
1894: for ( i = CNVars; --i >= 0; d1++, d2++ )
1895: if ( (*d1 <= *d2 ? *d1 : *d2) > 0 ) return 0;
1896: return 1;
1897: }
1898:
1.37 noro 1899: DP_pairs criterion_F( DP_pairs d1 )
1.1 noro 1900: {
1901: DP_pairs rest, head;
1902: register DP_pairs last, p, r, w;
1903: register int s;
1904:
1905: for ( head = last = 0, p = d1; NEXT(p); ) {
1906: s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
1907: while ( w = NEXT(w) )
1908: if ( criterion_2( w->dp1, w->dp2 ) ) {
1909: r = w;
1910: break;
1911: } else if ( w->sugar < s ) s = (r = w)->sugar;
1912: if ( last ) NEXT(last) = r;
1913: else head = r;
1914: NEXT(last = r) = 0;
1915: if ( !(p = rest) ) return head;
1916: }
1917: if ( !last ) return p;
1918: NEXT(last) = p;
1919: return head;
1920: }
1921:
1.37 noro 1922: NODE updbase(NODE g,int t)
1.1 noro 1923: {
1924: g = remove_reducibles(g,t);
1925: g = append_one(g,t);
1926: return g;
1927: }
1928:
1.37 noro 1929: NODE /* of index */ remove_reducibles(NODE /* of index */ nd, int newdp )
1.1 noro 1930: {
1931: register DL dl, dln;
1932: register NODE last, p, head;
1933: register int td;
1934:
1935: dl = psh[newdp];
1936: td = dl->td;
1937: for ( head = last = 0, p = nd; p; ) {
1938: dln = psh[(int)BDY(p)];
1939: if ( dln->td >= td && dl_redble( dln, dl ) ) {
1940: p = NEXT(p);
1941: if ( last ) NEXT(last) = p;
1942: } else {
1943: if ( !last ) head = p;
1944: p = NEXT(last = p);
1945: }
1946: }
1947: return head;
1948: }
1949:
1.37 noro 1950: int dl_redble(DL dl1,DL dl2)
1.1 noro 1951: {
1952: register int n, *d1, *d2;
1953:
1954: for ( d1 = dl1->d, d2 = dl2->d, n = CNVars; --n >= 0; d1++, d2++ )
1955: if ( *d1 < *d2 ) return 0;
1956: return 1;
1.5 noro 1957: }
1958:
1.41 noro 1959: #if 0
1960: int dl_weyl_weight(DL dl)
1.5 noro 1961: {
1962: int n,w,i;
1963:
1964: n = CNVars/2;
1965: for ( i = 0, w = 0; i < n; i++ )
1966: w += (-dl->d[i]+dl->d[n+i]);
1967: return w;
1.1 noro 1968: }
1.41 noro 1969: #endif
1.1 noro 1970:
1.37 noro 1971: int gbcheck(NODE f)
1.1 noro 1972: {
1973: int i;
1974: NODE r,g,gall;
1975: DP_pairs d,l;
1976: DP h,nf,dp1,dp2;
1977: struct oEGT tmp0,tmp1;
1978:
1979: if ( NoGC )
1980: return 1;
1981: for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
1982: i = (int)BDY(r);
1983: d = updpairs(d,g,i);
1984: g = updbase(g,i);
1985: gall = append_one(gall,i);
1986: }
1.12 noro 1987: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1988: fprintf(asir_out,"gbcheck total %d pairs\n",DPPlength(d)); fflush(asir_out);
1989: }
1990: while ( d ) {
1991: l = d; d = NEXT(d);
1992: get_eg(&tmp0);
1.17 noro 1993: dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
1994: dp_sp(dp1,dp2,&h);
1995: /* fprintf(stderr,"{%d,%d}",l->dp1,l->dp2); */
1.16 noro 1996: _dp_nf(gall,h,ps,1,&nf);
1.1 noro 1997: get_eg(&tmp1); add_eg(&eg_gc,&tmp0,&tmp1);
1.12 noro 1998: if ( DP_Print || DP_PrintShort ) {
1.1 noro 1999: fprintf(asir_out,"."); fflush(asir_out);
2000: }
2001: if ( nf )
2002: return 0;
2003: }
1.12 noro 2004: if ( DP_Print || DP_PrintShort )
1.1 noro 2005: fprintf(asir_out,"\n");
2006: return 1;
1.38 noro 2007: }
2008:
1.40 noro 2009: void gbcheck_list(NODE f,int n,VECT *gp,LIST *pp)
1.38 noro 2010: {
2011: int i;
2012: NODE r,g,gall,u,u0,t;
1.39 noro 2013: VECT vect;
1.38 noro 2014: LIST pair;
2015: DP_pairs d,l;
2016: Q q1,q2;
2017:
1.40 noro 2018: /* we need the following settings */
2019: NVars = CNVars = n;
1.39 noro 2020: setup_arrays(f,0,&r);
2021: for ( gall = g = 0, d = 0; r; r = NEXT(r) ) {
1.38 noro 2022: i = (int)BDY(r);
2023: d = updpairs(d,g,i);
2024: g = updbase(g,i);
2025: gall = append_one(gall,i);
2026: }
1.39 noro 2027: NEWVECT(vect); vect->len = psn; vect->body = (pointer)ps;
2028: *gp = vect;
2029:
1.38 noro 2030: for ( u0 = 0, l = d; l; l = NEXT(l) ) {
2031: NEXTNODE(u0,u);
2032: STOQ(l->dp1,q1);
2033: STOQ(l->dp2,q2);
2034: t = mknode(2,q1,q2);
2035: MKLIST(pair,t);
2036: BDY(u) = (pointer)pair;
2037: }
2038: if ( u0 )
2039: NEXT(u) = 0;
1.39 noro 2040: MKLIST(*pp,u0);
1.1 noro 2041: }
2042:
1.37 noro 2043: int membercheck(NODE f,NODE x)
1.1 noro 2044: {
2045: DP g;
2046: struct oEGT tmp0,tmp1;
2047:
2048: if ( NoMC )
2049: return 1;
1.12 noro 2050: if ( DP_Print || DP_PrintShort ) {
1.1 noro 2051: fprintf(asir_out,"membercheck\n"); fflush(asir_out);
2052: }
2053: for ( ; f; f = NEXT(f) ) {
2054: get_eg(&tmp0);
1.16 noro 2055: _dp_nf(x,(DP)BDY(f),ps,1,&g);
1.1 noro 2056: get_eg(&tmp1); add_eg(&eg_mc,&tmp0,&tmp1);
1.12 noro 2057: if ( DP_Print ) {
1.1 noro 2058: print_split_eg(&tmp0,&tmp1); fflush(asir_out);
1.12 noro 2059: } else if ( DP_PrintShort ) {
1.1 noro 2060: fprintf(asir_out,"."); fflush(asir_out);
2061: }
2062: if ( g )
2063: return 0;
2064: }
1.12 noro 2065: if ( DP_Print || DP_PrintShort )
1.1 noro 2066: fprintf(asir_out,"\n");
2067: return 1;
2068: }
2069:
1.37 noro 2070: void dp_set_flag(Obj name,Obj value)
1.1 noro 2071: {
2072: char *n;
2073: int v;
1.43 noro 2074: Q ratio;
1.1 noro 2075:
2076: if ( OID(name) != O_STR )
2077: return;
2078: n = BDY((STRING)name);
2079: if ( !strcmp(n,"Demand") ) {
2080: Demand = value ? BDY((STRING)value) : 0; return;
2081: }
2082: if ( !strcmp(n,"Dist") ) {
2083: Dist = (LIST)value; return;
2084: }
1.43 noro 2085: if ( !strcmp(n,"Content") ) {
2086: ratio = (Q)value;
2087: if ( ratio ) {
2088: DP_Multiple = BD(NM(ratio))[0];
2089: Denominator = INT(ratio) ? 1 : BD(DN(ratio))[0];
2090: } else {
2091: DP_Multiple = 0;
2092: Denominator = 1;
2093: }
2094: }
1.1 noro 2095: if ( value && OID(value) != O_N )
2096: return;
2097: v = QTOS((Q)value);
2098: if ( !strcmp(n,"NoSugar") )
2099: NoSugar = v;
2100: else if ( !strcmp(n,"NoCriB") )
2101: NoCriB = v;
2102: else if ( !strcmp(n,"NoGC") )
2103: NoGC = v;
2104: else if ( !strcmp(n,"NoMC") )
2105: NoMC = v;
2106: else if ( !strcmp(n,"NoRA") )
2107: NoRA = v;
2108: else if ( !strcmp(n,"NoGCD") )
2109: NoGCD = v;
2110: else if ( !strcmp(n,"Top") )
2111: Top = v;
2112: else if ( !strcmp(n,"ShowMag") )
2113: ShowMag = v;
1.14 noro 2114: else if ( !strcmp(n,"PrintShort") )
1.12 noro 2115: DP_PrintShort = v;
1.14 noro 2116: else if ( !strcmp(n,"Print") )
1.12 noro 2117: DP_Print = v;
1.14 noro 2118: else if ( !strcmp(n,"NFStat") )
2119: DP_NFStat = v;
1.1 noro 2120: else if ( !strcmp(n,"Stat") )
2121: Stat = v;
2122: else if ( !strcmp(n,"Reverse") )
2123: Reverse = v;
2124: else if ( !strcmp(n,"Multiple") )
1.14 noro 2125: DP_Multiple = v;
1.1 noro 2126: else if ( !strcmp(n,"Denominator") )
2127: Denominator = v;
2128: else if ( !strcmp(n,"PtozpRA") )
2129: PtozpRA = v;
1.7 noro 2130: else if ( !strcmp(n,"GenTrace") )
2131: GenTrace = v;
2132: else if ( !strcmp(n,"OXCheck") )
2133: OXCheck = v;
1.1 noro 2134: }
2135:
1.37 noro 2136: void dp_make_flaglist(LIST *list)
1.1 noro 2137: {
1.43 noro 2138: Q v,nm,dn;
1.1 noro 2139: STRING name,path;
2140: NODE n,n1;
2141:
1.43 noro 2142: #if 0
1.14 noro 2143: STOQ(DP_Multiple,v); MKNODE(n,v,0); MKSTR(name,"DP_Multiple"); MKNODE(n1,name,n); n = n1;
1.1 noro 2144: STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
1.43 noro 2145: #else
2146: if ( DP_Multiple ) {
2147: STOQ(DP_Multiple,nm); STOQ(Denominator,dn); divq(nm,dn,&v);
2148: } else
2149: v = 0;
2150: MKNODE(n,v,0); MKSTR(name,"Content"); MKNODE(n1,name,n); n = n1;
2151: #endif
1.1 noro 2152: MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
2153: STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
2154: STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
1.14 noro 2155: STOQ(DP_Print,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Print"); MKNODE(n1,name,n); n = n1;
1.20 noro 2156: STOQ(DP_PrintShort,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PrintShort"); MKNODE(n1,name,n); n = n1;
1.14 noro 2157: STOQ(DP_NFStat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NFStat"); MKNODE(n1,name,n); n = n1;
1.7 noro 2158: STOQ(OXCheck,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"OXCheck"); MKNODE(n1,name,n); n = n1;
2159: STOQ(GenTrace,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenTrace"); MKNODE(n1,name,n); n = n1;
1.1 noro 2160: STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
2161: STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
2162: STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
2163: STOQ(NoGCD,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGCD"); MKNODE(n1,name,n); n = n1;
2164: STOQ(NoRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoRA"); MKNODE(n1,name,n); n = n1;
2165: STOQ(NoMC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoMC"); MKNODE(n1,name,n); n = n1;
2166: STOQ(NoGC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGC"); MKNODE(n1,name,n); n = n1;
2167: STOQ(NoCriB,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoCriB"); MKNODE(n1,name,n); n = n1;
2168: STOQ(NoSugar,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoSugar"); MKNODE(n1,name,n); n = n1;
2169: if ( Demand )
2170: MKSTR(path,Demand);
2171: else
2172: path = 0;
2173: MKNODE(n1,path,n); n = n1; MKSTR(name,"Demand"); MKNODE(n1,name,n); n = n1;
2174: MKLIST(*list,n);
2175: }
2176:
2177: #define DELIM '/'
2178:
1.37 noro 2179: void dp_save(int index,Obj p,char *prefix)
1.1 noro 2180: {
2181: FILE *fp;
2182: char path[BUFSIZ];
2183:
1.7 noro 2184: if ( prefix )
2185: sprintf(path,"%s%c%s%d",Demand,DELIM,prefix,index);
2186: else
2187: sprintf(path,"%s%c%d",Demand,DELIM,index);
1.1 noro 2188: if ( !(fp = fopen(path,"wb") ) )
2189: error("dp_save : cannot open a file");
1.7 noro 2190: savevl(fp,VC); saveobj(fp,p); fclose(fp);
1.1 noro 2191: }
2192:
1.37 noro 2193: void dp_load(int index,DP *p)
1.1 noro 2194: {
2195: FILE *fp;
2196: char path[BUFSIZ];
2197:
2198: if ( !Demand || ps[index] )
2199: *p = ps[index];
2200: else {
2201: sprintf(path,"%s%c%d",Demand,DELIM,index);
2202: if ( !(fp = fopen(path,"rb") ) )
2203: error("dp_load : cannot open a file");
1.45 noro 2204: if ( PCoeffs )
2205: loadvl(fp);
2206: else
2207: skipvl(fp);
2208: loadobj(fp,(Obj *)p); fclose(fp);
1.1 noro 2209: }
2210: }
2211:
1.37 noro 2212: int dp_load_t(int index,DP *p)
1.1 noro 2213: {
2214: FILE *fp;
2215: char path[BUFSIZ];
2216:
2217: sprintf(path,"%s%c%d",Demand,DELIM,index);
2218: if ( !(fp = fopen(path,"rb") ) )
2219: return 0;
2220: else {
1.45 noro 2221: if ( PCoeffs )
2222: loadvl(fp);
2223: else
2224: skipvl(fp);
2225: loadobj(fp,(Obj *)p); fclose(fp); return 1;
1.1 noro 2226: }
2227: }
2228:
2229: void init_stat() {
2230: init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
1.14 noro 2231: init_eg(&eg_pz); init_eg(&eg_np);
1.1 noro 2232: init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
1.35 noro 2233: ZR = NZR = TP = NMP = NBP = NFP = NDP = 0;
1.1 noro 2234: }
2235:
2236: void print_stat() {
1.12 noro 2237: if ( !DP_Print && !Stat )
1.1 noro 2238: return;
2239: print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
1.14 noro 2240: print_eg("PZ",&eg_pz); print_eg("NP",&eg_np);
1.1 noro 2241: print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
2242: 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);
2243: }
1.14 noro 2244:
2245: /*
2246: * dp_nf used in gb()
2247: *
2248: */
2249:
2250: double pz_t_e, pz_t_d, pz_t_d1, pz_t_c, im_t_s, im_t_r;
2251:
2252: extern int GenTrace;
2253: extern NODE TraceList;
1.15 noro 2254: extern int mpi_mag;
2255:
1.37 noro 2256: void dp_mulc_d(DP p,P c,DP *r)
1.15 noro 2257: {
2258: if ( Dist && BDY(Dist)
2259: && HMAG(p) > mpi_mag
2260: && p_mag((P)c) > mpi_mag ) {
2261: if ( DP_NFStat ) fprintf(asir_out,"~");
2262: dp_imul_d(p,(Q)c,r);
2263: } else {
2264: if ( DP_NFStat ) fprintf(asir_out,"_");
2265: muldc(CO,p,c,r);
2266: }
2267: }
1.14 noro 2268:
1.37 noro 2269: void _dp_nf(NODE b,DP g,DP *ps,int full,DP *rp)
1.16 noro 2270: {
2271: DP u,p,d,s,t,mult;
2272: P coef;
2273: NODE l;
2274: MP m,mr;
2275: int sugar,psugar;
2276:
2277: if ( !g ) {
2278: *rp = 0; return;
2279: }
2280: sugar = g->sugar;
2281: for ( d = 0; g; ) {
2282: for ( u = 0, l = b; l; l = NEXT(l) ) {
2283: if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
2284: dp_load((int)BDY(l),&p);
2285: /* t+u = coef*(d+g) - mult*p (t = coef*d) */
2286: dp_red(d,g,p,&t,&u,&coef,&mult);
2287: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
2288: sugar = MAX(sugar,psugar);
2289: if ( GenTrace ) {
2290: LIST hist;
2291: Q cq;
2292: NODE node,node0;
2293:
2294: STOQ((int)BDY(l),cq);
2295: node0 = mknode(4,coef,cq,mult,ONE);
2296: MKLIST(hist,node0);
2297: MKNODE(node,hist,TraceList); TraceList = node;
2298: }
2299: if ( !u ) {
2300: if ( d )
2301: d->sugar = sugar;
2302: *rp = d; return;
2303: }
2304: d = t;
2305: break;
2306: }
2307: }
2308: if ( u )
2309: g = u;
2310: else if ( !full ) {
2311: if ( g ) {
2312: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
2313: }
2314: *rp = g; return;
2315: } else {
2316: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
2317: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
2318: addd(CO,d,t,&s); d = s;
2319: dp_rest(g,&t); g = t;
2320: }
2321: }
2322: if ( d )
2323: d->sugar = sugar;
2324: *rp = d;
2325: }
2326:
1.37 noro 2327: void _dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *r)
1.14 noro 2328: {
1.37 noro 2329: DP u,dp,rp,t,t1,red,shift;
1.14 noro 2330: Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
2331: NODE l;
2332: int hmag,denom;
2333: int sugar,psugar;
2334: STRING imul;
1.37 noro 2335: double t_0,tt,t_p,t_m,t_g,t_a;
1.14 noro 2336: LIST hist;
2337: NODE node;
2338: Q rcred,mrcred;
2339:
2340: if ( !g ) {
2341: *r = 0; return;
2342: }
2343: pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1.15 noro 2344: t_p = t_m = t_g = t_a = 0;
1.14 noro 2345:
2346: denom = Denominator?Denominator:1;
2347: hmag = multiple*HMAG(g)/denom;
2348: sugar = g->sugar;
2349:
2350: dc = 0; dp = 0; rc = ONE; rp = g;
2351: MKSTR(imul,"dp_imul_index");
2352:
2353: /* g = dc*dp+rc*rp */
2354: for ( ; rp; ) {
2355: for ( u = 0, l = b; l; l = NEXT(l) ) {
2356: if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
2357: t_0 = get_rtime();
2358: dp_load((int)BDY(l),&red);
2359: hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
2360: igcd_cofactor((Q)BDY(rp)->c,(Q)BDY(red)->c,&gcd,&cred,&cr);
1.15 noro 2361: tt = get_rtime(); t_p += tt-t_0;
1.14 noro 2362:
2363: dp_subd(rp,red,&shift);
1.37 noro 2364: dp_mulc_d(rp,(P)cr,&t);
1.14 noro 2365: chsgnp((P)cred,(P *)&mcred);
1.37 noro 2366: dp_mulc_d(red,(P)mcred,&t1);
1.15 noro 2367: muld(CO,shift,t1,&t1);
2368: addd(CO,t,t1,&u);
2369: t_m += get_rtime()-tt;
1.14 noro 2370:
2371: psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
2372: sugar = MAX(sugar,psugar);
1.15 noro 2373:
1.14 noro 2374: if ( GenTrace ) {
2375: /* u = cr*rp + (-cred)*shift*red */
2376: STOQ((int)BDY(l),cq);
2377: node = mknode(4,cr,cq,0,0);
2378: mulq(cred,rc,&rcred);
2379: chsgnnum((Num)rcred,(Num *)&mrcred);
2380: muldc(CO,shift,(P)mrcred,(DP *)&ARG2(node));
2381: MKLIST(hist,node);
2382: }
1.15 noro 2383:
1.14 noro 2384: if ( !u ) {
2385: if ( dp )
2386: dp->sugar = sugar;
2387: *r = dp;
2388: if ( GenTrace ) {
2389: ARG3(BDY(hist)) = ONE;
2390: MKNODE(node,hist,TraceList); TraceList = node;
2391: }
2392: goto final;
2393: }
2394: break;
2395: }
2396: }
2397: if ( u ) {
2398: if ( multiple && HMAG(u) > hmag ) {
2399: t_0 = get_rtime();
1.16 noro 2400: dp_ptozp_d(u,&rp);
1.15 noro 2401: tt = get_rtime(); t_g += tt-t_0;
2402:
2403: divsq((Q)BDY(u)->c,(Q)BDY(rp)->c,&cont);
1.14 noro 2404: if ( !dp_fcoeffs && DP_NFStat ) {
1.15 noro 2405: fprintf(asir_out,
2406: "(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1.14 noro 2407: fflush(asir_out);
2408: }
1.15 noro 2409: mulq(cr,dc,&dcq); mulq(cont,rc,&rcq);
1.14 noro 2410: igcd_cofactor(dcq,rcq,&gcd,&dc,&rc);
1.15 noro 2411: t_a = get_rtime()-tt;
2412:
1.14 noro 2413: hmag = multiple*HMAG(rp)/denom;
2414: if ( GenTrace ) {
2415: ARG3(BDY(hist)) = (pointer)gcd;
2416: MKNODE(node,hist,TraceList); TraceList = node;
2417: }
2418: } else {
1.15 noro 2419: rp = u;
1.14 noro 2420: t_0 = get_rtime();
1.15 noro 2421: mulq(cr,dc,&dc);
2422: t_a += get_rtime()-t_0;
1.14 noro 2423: if ( GenTrace ) {
2424: ARG3(BDY(hist)) = (pointer)ONE;
2425: MKNODE(node,hist,TraceList); TraceList = node;
2426: }
2427: }
2428: } else if ( !full ) {
2429: if ( rp ) {
2430: MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
2431: }
2432: *r = rp;
2433: goto final;
2434: } else {
2435: t_0 = get_rtime();
2436: mulq((Q)BDY(rp)->c,rc,&c);
1.15 noro 2437: igcd_cofactor(dc,c,&dc,&dcq,&cq);
2438: muldc(CO,dp,(P)dcq,&t);
2439: dp_hm(rp,&t1); BDY(t1)->c = (P)cq; addd(CO,t,t1,&dp);
2440: dp_rest(rp,&rp);
2441: t_a += get_rtime()-t_0;
1.14 noro 2442: }
2443: }
2444: if ( GenTrace ) {
2445: mulq(ARG3(BDY((LIST)BDY(TraceList))),dc,&cq);
2446: ARG3(BDY((LIST)BDY(TraceList))) = (pointer)cq;
2447: }
2448: if ( dp )
2449: dp->sugar = sugar;
2450: *r = dp;
2451: final:
2452: if ( DP_NFStat )
1.15 noro 2453: fprintf(asir_out,
2454: "(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
2455: t_p,t_m,t_g,t_a,
1.14 noro 2456: pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
2457: }
2458:
2459: void imulv();
2460:
1.37 noro 2461: void dp_imul_d(DP p,Q q,DP *rp)
1.14 noro 2462: {
2463: int nsep,ndist,i,j,k,l,n;
2464: double t0,t1,t2;
2465: Q *s;
2466: pointer *b;
2467: VECT c,cs,ri;
2468: VECT *r;
2469: MP m;
1.37 noro 2470: NODE tn,dist,n0;
2471: Obj obj;
1.14 noro 2472: STRING imul;
2473: extern LIST Dist;
2474:
2475: if ( !p || !q ) {
2476: *rp = 0; return;
2477: }
2478: dist = BDY(Dist);
2479: for ( tn = dist, ndist = 0; tn; tn = NEXT(tn), ndist++ );
2480: nsep = ndist + 1;
2481: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
2482: if ( n <= nsep ) {
2483: muldc(CO,p,(P)q,rp); return;
2484: }
2485: MKSTR(imul,"imulv");
2486: t0 = get_rtime();
2487: dp_dtov(p,&c);
2488: sepvect(c,nsep,&cs);
2489: r = (VECT *)CALLOC(nsep,sizeof(VECT *));
2490: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
2491: n0 = mknode(4,BDY(tn),imul,b[i],q);
1.37 noro 2492: Pox_rpc(n0,&obj);
1.14 noro 2493: }
2494: t1 = get_rtime();
2495: im_t_s += t1 - t0;
2496: imulv(b[i],q,&r[i]);
2497: t1 = get_rtime();
2498: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
2499: MKNODE(n0,BDY(tn),0);
1.37 noro 2500: Pox_pop_local(n0,&obj); r[i] = (VECT)obj;
1.14 noro 2501: if ( OID(r[i]) == O_ERR ) {
2502: printexpr(CO,(Obj)r[i]);
2503: error("dp_imul_d : aborted");
2504: }
2505: }
2506: t2 = get_rtime();
2507: im_t_r += t2 - t1;
2508: s = (Q *)CALLOC(n,sizeof(Q));
2509: for ( i = j = 0; i < nsep; i++ ) {
2510: for ( k = 0, ri = r[i], l = ri->len; k < l; k++, j++ ) {
2511: s[j] = (Q)BDY(ri)[k];
2512: }
2513: }
2514: dp_vtod(s,p,rp);
2515: }
2516:
1.37 noro 2517: void imulv(VECT w,Q c,VECT *rp)
1.14 noro 2518: {
2519: int n,i;
2520: VECT r;
2521:
2522: n = w->len;
2523: MKVECT(r,n); *rp = r;
2524: for ( i = 0; i < n; i++ )
2525: mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
1.32 noro 2526: }
2527:
1.37 noro 2528: void dptoca(DP p,unsigned int **rp)
1.32 noro 2529: {
2530: int i;
2531: MP m;
2532: unsigned int *r;
2533:
2534: if ( !p )
2535: *rp = 0;
2536: else {
2537: for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ );
2538: *rp = r = (unsigned int *)MALLOC_ATOMIC(i*sizeof(unsigned int));
2539: for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ )
2540: r[i] = ITOS(C(m));
2541: }
1.14 noro 2542: }
2543:
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