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