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