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