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