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