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