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