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