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