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