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