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