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