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