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