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