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