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