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