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