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