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