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