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