![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dp-supp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin|
Return to dp-supp.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.2 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.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 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/dp-supp.c,v 1.65 2017/03/27 09:05:46 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
1.16 noro 52: #include "inline.h"
1.1 noro 53: #include "parse.h"
54: #include "ox.h"
55:
1.66 ! noro 56: #define HMAG(p) (p_mag((P)BDY(p)->c))
1.5 noro 57:
1.1 noro 58: extern int (*cmpdl)();
1.5 noro 59: extern double pz_t_e,pz_t_d,pz_t_d1,pz_t_c;
60: extern int dp_nelim,dp_fcoeffs;
1.7 noro 61: extern int NoGCD;
62: extern int GenTrace;
63: extern NODE TraceList;
64:
1.37 noro 65: int show_orderspec;
66:
67: void print_composite_order_spec(struct order_spec *spec);
1.66 ! noro 68: void dpm_rest(DPM,DPM *);
1.37 noro 69:
1.7 noro 70: /*
71: * content reduction
72: *
73: */
74:
1.46 noro 75: static NODE RatDenomList;
76:
77: void init_denomlist()
78: {
79: RatDenomList = 0;
80: }
81:
82: void add_denomlist(P f)
83: {
84: NODE n;
85:
86: if ( OID(f)==O_P ) {
87: MKNODE(n,f,RatDenomList); RatDenomList = n;
88: }
89: }
90:
91: LIST get_denomlist()
92: {
93: LIST l;
94:
95: MKLIST(l,RatDenomList); RatDenomList = 0;
96: return l;
97: }
98:
1.20 noro 99: void dp_ptozp(DP p,DP *rp)
1.7 noro 100: {
101: MP m,mr,mr0;
102: int i,n;
103: Q *w;
104: Q dvr;
105: P t;
106:
107: if ( !p )
108: *rp = 0;
109: else {
110: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
111: w = (Q *)ALLOCA(n*sizeof(Q));
112: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
113: if ( NUM(m->c) )
114: w[i] = (Q)m->c;
115: else
1.66 ! noro 116: ptozp((P)m->c,1,&w[i],&t);
1.7 noro 117: sortbynm(w,n);
118: qltozl(w,n,&dvr);
119: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1.66 ! noro 120: NEXTMP(mr0,mr); divsp(CO,(P)m->c,(P)dvr,(P *)&mr->c); mr->dl = m->dl;
1.7 noro 121: }
122: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
123: }
124: }
125:
1.20 noro 126: void dp_ptozp2(DP p0,DP p1,DP *hp,DP *rp)
1.7 noro 127: {
128: DP t,s,h,r;
129: MP m,mr,mr0,m0;
130:
131: addd(CO,p0,p1,&t); dp_ptozp(t,&s);
132: if ( !p0 ) {
133: h = 0; r = s;
134: } else if ( !p1 ) {
135: h = s; r = 0;
136: } else {
137: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
138: m = NEXT(m), m0 = NEXT(m0) ) {
139: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
140: }
141: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
142: }
143: if ( h )
144: h->sugar = p0->sugar;
145: if ( r )
146: r->sugar = p1->sugar;
147: *hp = h; *rp = r;
1.39 ohara 148: }
149:
1.66 ! noro 150: void dpm_ptozp(DPM p,DPM *rp)
! 151: {
! 152: DMM m,mr,mr0;
! 153: int i,n;
! 154: Q *w;
! 155: Q dvr;
! 156: P t;
! 157:
! 158: if ( !p )
! 159: *rp = 0;
! 160: else {
! 161: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
! 162: w = (Q *)ALLOCA(n*sizeof(Q));
! 163: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
! 164: if ( NUM(m->c) )
! 165: w[i] = (Q)m->c;
! 166: else
! 167: ptozp((P)m->c,1,&w[i],&t);
! 168: sortbynm(w,n);
! 169: qltozl(w,n,&dvr);
! 170: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
! 171: NEXTDMM(mr0,mr); divsp(CO,(P)m->c,(P)dvr,(P *)&mr->c); mr->dl = m->dl; mr->pos = m->pos;
! 172: }
! 173: NEXT(mr) = 0; MKDPM(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
! 174: }
! 175: }
! 176:
! 177: void dpm_ptozp2(DPM p0,DPM p1,DPM *hp,DPM *rp)
! 178: {
! 179: DPM t,s,h,r;
! 180: DMM m,mr,mr0,m0;
! 181:
! 182: adddpm(CO,p0,p1,&t); dpm_ptozp(t,&s);
! 183: if ( !p0 ) {
! 184: h = 0; r = s;
! 185: } else if ( !p1 ) {
! 186: h = s; r = 0;
! 187: } else {
! 188: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
! 189: m = NEXT(m), m0 = NEXT(m0) ) {
! 190: NEXTDMM(mr0,mr); mr->c = m->c; mr->dl = m->dl; mr->pos = m->pos;
! 191: }
! 192: NEXT(mr) = 0; MKDPM(p0->nv,mr0,h); MKDPM(p0->nv,m,r);
! 193: }
! 194: if ( h )
! 195: h->sugar = p0->sugar;
! 196: if ( r )
! 197: r->sugar = p1->sugar;
! 198: *hp = h; *rp = r;
! 199: }
! 200:
! 201:
1.39 ohara 202: void dp_ptozp3(DP p,Q *dvr,DP *rp)
203: {
204: MP m,mr,mr0;
205: int i,n;
206: Q *w;
207: P t;
208:
209: if ( !p ) {
210: *rp = 0; *dvr = 0;
211: }else {
212: for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ );
213: w = (Q *)ALLOCA(n*sizeof(Q));
214: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
215: if ( NUM(m->c) )
216: w[i] = (Q)m->c;
217: else
1.66 ! noro 218: ptozp((P)m->c,1,&w[i],&t);
1.39 ohara 219: sortbynm(w,n);
220: qltozl(w,n,dvr);
221: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1.66 ! noro 222: NEXTMP(mr0,mr); divsp(CO,(P)m->c,(P)(*dvr),(P *)&mr->c); mr->dl = m->dl;
1.39 ohara 223: }
224: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
225: }
1.7 noro 226: }
1.1 noro 227:
1.20 noro 228: void dp_idiv(DP p,Q c,DP *rp)
1.1 noro 229: {
230: Q t;
231: N nm,q;
232: int sgn,s;
233: MP mr0,m,mr;
234:
235: if ( !p )
236: *rp = 0;
237: else if ( MUNIQ((Q)c) )
238: *rp = p;
239: else if ( MUNIQ((Q)c) )
240: chsgnd(p,rp);
241: else {
242: nm = NM(c); sgn = SGN(c);
243: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
244: NEXTMP(mr0,mr);
245:
246: divsn(NM((Q)(m->c)),nm,&q);
247: s = sgn*SGN((Q)(m->c));
248: NTOQ(q,s,t);
1.66 ! noro 249: mr->c = (Obj)t;
1.1 noro 250: mr->dl = m->dl;
251: }
252: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
253: if ( *rp )
254: (*rp)->sugar = p->sugar;
255: }
256: }
257:
1.20 noro 258: void dp_mbase(NODE hlist,NODE *mbase)
1.1 noro 259: {
260: DL *dl;
261: DL d;
1.65 noro 262: int *t;
263: int i,j,k,n,nvar,td;
1.1 noro 264:
265: n = length(hlist); nvar = ((DP)BDY(hlist))->nv;
266: dl = (DL *)MALLOC(n*sizeof(DL));
1.65 noro 267: NEWDL(d,nvar); *mbase = 0;
268: for ( i = 0; i < n; i++, hlist = NEXT(hlist) ) {
1.1 noro 269: dl[i] = BDY((DP)BDY(hlist))->dl;
1.65 noro 270: /* trivial ideal check */
271: if ( (*cmpdl)(nvar,d,dl[i]) == 0 ) {
272: return;
273: }
274: }
275: /* zero-dim. ideal check */
276: for ( i = 0; i < nvar; i++ ) {
277: for ( j = 0; j < n; j++ ) {
278: for ( k = 0, t = dl[j]->d; k < nvar; k++ )
279: if ( k != i && t[k] != 0 ) break;
280: if ( k == nvar ) break;
281: }
282: if ( j == n )
283: error("dp_mbase : input ideal is not zero-dimensional");
284: }
1.1 noro 285: while ( 1 ) {
286: insert_to_node(d,mbase,nvar);
287: for ( i = nvar-1; i >= 0; ) {
1.21 noro 288: d->d[i]++;
289: d->td += MUL_WEIGHT(1,i);
1.1 noro 290: for ( j = 0; j < n; j++ ) {
291: if ( _dl_redble(dl[j],d,nvar) )
292: break;
293: }
294: if ( j < n ) {
295: for ( j = nvar-1; j >= i; j-- )
296: d->d[j] = 0;
297: for ( j = 0, td = 0; j < i; j++ )
1.21 noro 298: td += MUL_WEIGHT(d->d[j],j);
1.1 noro 299: d->td = td;
300: i--;
301: } else
302: break;
303: }
304: if ( i < 0 )
305: break;
306: }
307: }
308:
1.20 noro 309: int _dl_redble(DL d1,DL d2,int nvar)
1.1 noro 310: {
311: int i;
312:
313: if ( d1->td > d2->td )
314: return 0;
315: for ( i = 0; i < nvar; i++ )
316: if ( d1->d[i] > d2->d[i] )
317: break;
318: if ( i < nvar )
319: return 0;
320: else
321: return 1;
322: }
323:
1.20 noro 324: void insert_to_node(DL d,NODE *n,int nvar)
1.1 noro 325: {
326: DL d1;
327: MP m;
328: DP dp;
329: NODE n0,n1,n2;
330:
331: NEWDL(d1,nvar); d1->td = d->td;
332: bcopy((char *)d->d,(char *)d1->d,nvar*sizeof(int));
1.66 ! noro 333: NEWMP(m); m->dl = d1; m->c = (Obj)ONE; NEXT(m) = 0;
1.1 noro 334: MKDP(nvar,m,dp); dp->sugar = d->td;
335: if ( !(*n) ) {
336: MKNODE(n1,dp,0); *n = n1;
337: } else {
338: for ( n1 = *n, n0 = 0; n1; n0 = n1, n1 = NEXT(n1) )
339: if ( (*cmpdl)(nvar,d,BDY((DP)BDY(n1))->dl) > 0 ) {
340: MKNODE(n2,dp,n1);
341: if ( !n0 )
342: *n = n2;
343: else
344: NEXT(n0) = n2;
345: break;
346: }
347: if ( !n1 ) {
348: MKNODE(n2,dp,0); NEXT(n0) = n2;
349: }
350: }
351: }
352:
1.20 noro 353: void dp_vtod(Q *c,DP p,DP *rp)
1.1 noro 354: {
355: MP mr0,m,mr;
356: int i;
357:
358: if ( !p )
359: *rp = 0;
360: else {
361: for ( mr0 = 0, m = BDY(p), i = 0; m; m = NEXT(m), i++ ) {
1.66 ! noro 362: NEXTMP(mr0,mr); mr->c = (Obj)c[i]; mr->dl = m->dl;
1.1 noro 363: }
364: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp);
365: (*rp)->sugar = p->sugar;
366: }
367: }
368:
1.8 noro 369: extern int mpi_mag;
370: extern int PCoeffs;
371:
1.20 noro 372: void dp_ptozp_d(DP p,DP *rp)
1.1 noro 373: {
374: int i,j,k,l,n,nsep;
375: MP m;
376: NODE tn,n0,n1,n2,n3;
377: struct oVECT v;
378: VECT c,cs;
379: VECT qi,ri;
380: LIST *qr;
381: Obj dmy;
382: Q d0,d1,gcd,a,u,u1;
383: Q *q,*r;
384: STRING iqr_v;
385: pointer *b;
386: N qn,gn;
387: double get_rtime();
388: int blen;
1.8 noro 389: NODE dist;
390: int ndist;
1.1 noro 391: double t0;
392: double t_e,t_d,t_d1,t_c;
1.8 noro 393: extern int DP_NFStat;
394: extern LIST Dist;
1.20 noro 395: void Pox_rpc();
396: void Pox_pop_local();
1.1 noro 397:
398: if ( !p )
399: *rp = 0;
400: else {
1.8 noro 401: if ( PCoeffs ) {
402: dp_ptozp(p,rp); return;
403: }
1.66 ! noro 404: if ( !Dist || p_mag((P)BDY(p)->c) <= mpi_mag ) {
1.8 noro 405: dist = 0; ndist = 0;
406: if ( DP_NFStat ) fprintf(asir_out,"L");
407: } else {
408: dist = BDY(Dist); ndist = length(dist);
409: if ( DP_NFStat ) fprintf(asir_out,"D");
410: }
1.1 noro 411: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
412: nsep = ndist + 1;
413: if ( n <= nsep ) {
414: dp_ptozp(p,rp); return;
415: }
416: t0 = get_rtime();
417: dp_dtov(p,&c);
418: igcdv_estimate(c,&d0);
419: t_e = get_rtime()-t0;
420: t0 = get_rtime();
421: dp_dtov(p,&c);
422: sepvect(c,nsep,&cs);
423: MKSTR(iqr_v,"iqr");
424: qr = (LIST *)CALLOC(nsep,sizeof(LIST));
425: q = (Q *)CALLOC(n,sizeof(Q));
426: r = (Q *)CALLOC(n,sizeof(Q));
427: for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
428: MKNODE(n3,d0,0); MKNODE(n2,b[i],n3);
429: MKNODE(n1,iqr_v,n2); MKNODE(n0,BDY(tn),n1);
430: Pox_rpc(n0,&dmy);
431: }
432: iqrv(b[i],d0,&qr[i]);
433: dp_dtov(p,&c);
434: for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
435: Pox_pop_local(tn,&qr[i]);
436: if ( OID(qr[i]) == O_ERR ) {
437: printexpr(CO,(Obj)qr[i]);
438: error("dp_ptozp_d : aborted");
439: }
440: }
441: t_d = get_rtime()-t0;
442: t_d1 = t_d/n;
443: t0 = get_rtime();
444: for ( i = j = 0; i < nsep; i++ ) {
445: tn = BDY(qr[i]); qi = (VECT)BDY(tn); ri = (VECT)BDY(NEXT(tn));
446: for ( k = 0, l = qi->len; k < l; k++, j++ ) {
447: q[j] = (Q)BDY(qi)[k]; r[j] = (Q)BDY(ri)[k];
448: }
449: }
450: v.id = O_VECT; v.len = n; v.body = (pointer *)r; igcdv(&v,&d1);
451: if ( d1 ) {
452: gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
453: divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
454: for ( i = 0; i < n; i++ ) {
455: mulq(a,q[i],&u);
456: if ( r[i] ) {
457: divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
458: addq(u,u1,&q[i]);
459: } else
460: q[i] = u;
461: }
462: } else
463: gcd = d0;
464: dp_vtod(q,p,rp);
465: t_c = get_rtime()-t0;
466: blen=p_mag((P)gcd);
467: pz_t_e += t_e; pz_t_d += t_d; pz_t_d1 += t_d1; pz_t_c += t_c;
468: if ( 0 )
469: fprintf(stderr,"(%d,%d)",p_mag((P)d0)-blen,blen);
470: }
471: }
472:
1.20 noro 473: void dp_ptozp2_d(DP p0,DP p1,DP *hp,DP *rp)
1.1 noro 474: {
475: DP t,s,h,r;
476: MP m,mr,mr0,m0;
477:
1.8 noro 478: addd(CO,p0,p1,&t); dp_ptozp_d(t,&s);
1.1 noro 479: if ( !p0 ) {
480: h = 0; r = s;
481: } else if ( !p1 ) {
482: h = s; r = 0;
483: } else {
484: for ( mr0 = 0, m = BDY(s), m0 = BDY(p0); m0;
485: m = NEXT(m), m0 = NEXT(m0) ) {
486: NEXTMP(mr0,mr); mr->c = m->c; mr->dl = m->dl;
487: }
488: NEXT(mr) = 0; MKDP(p0->nv,mr0,h); MKDP(p0->nv,m,r);
489: }
490: if ( h )
491: h->sugar = p0->sugar;
492: if ( r )
493: r->sugar = p1->sugar;
494: *hp = h; *rp = r;
1.5 noro 495: }
496:
1.22 noro 497: int have_sf_coef(P p)
498: {
499: DCP dc;
500:
501: if ( !p )
502: return 0;
503: else if ( NUM(p) )
504: return NID((Num)p) == N_GFS ? 1 : 0;
505: else {
506: for ( dc = DC(p); dc; dc = NEXT(dc) )
507: if ( have_sf_coef(COEF(dc)) )
508: return 1;
509: return 0;
510: }
511: }
512:
1.25 noro 513: void head_coef(P p,Num *c)
514: {
515: if ( !p )
516: *c = 0;
517: else if ( NUM(p) )
518: *c = (Num)p;
519: else
520: head_coef(COEF(DC(p)),c);
521: }
522:
523: void dp_monic_sf(DP p,DP *rp)
524: {
525: Num c;
526:
527: if ( !p )
528: *rp = 0;
529: else {
1.66 ! noro 530: head_coef((P)BDY(p)->c,&c);
1.25 noro 531: divsdc(CO,p,(P)c,rp);
532: }
533: }
534:
1.20 noro 535: void dp_prim(DP p,DP *rp)
1.5 noro 536: {
1.7 noro 537: P t,g;
538: DP p1;
539: MP m,mr,mr0;
540: int i,n;
541: P *w;
542: Q *c;
543: Q dvr;
1.46 noro 544: NODE tn;
1.5 noro 545:
1.7 noro 546: if ( !p )
547: *rp = 0;
1.23 noro 548: else if ( dp_fcoeffs == N_GFS ) {
549: for ( m = BDY(p); m; m = NEXT(m) )
1.22 noro 550: if ( OID(m->c) == O_N ) {
551: /* GCD of coeffs = 1 */
1.25 noro 552: dp_monic_sf(p,rp);
1.22 noro 553: return;
1.23 noro 554: } else break;
555: /* compute GCD over the finite fieid */
556: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
557: w = (P *)ALLOCA(n*sizeof(P));
558: for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ )
1.66 ! noro 559: w[i] = (P)m->c;
1.23 noro 560: gcdsf(CO,w,n,&g);
561: if ( NUM(g) )
1.25 noro 562: dp_monic_sf(p,rp);
1.23 noro 563: else {
564: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1.66 ! noro 565: NEXTMP(mr0,mr); divsp(CO,(P)m->c,g,(P *)&mr->c); mr->dl = m->dl;
1.22 noro 566: }
1.25 noro 567: NEXT(mr) = 0; MKDP(p->nv,mr0,p1); p1->sugar = p->sugar;
568: dp_monic_sf(p1,rp);
1.22 noro 569: }
1.23 noro 570: return;
571: } else if ( dp_fcoeffs )
1.7 noro 572: *rp = p;
1.23 noro 573: else if ( NoGCD )
1.7 noro 574: dp_ptozp(p,rp);
575: else {
576: dp_ptozp(p,&p1); p = p1;
577: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
578: if ( n == 1 ) {
579: m = BDY(p);
1.66 ! noro 580: NEWMP(mr); mr->dl = m->dl; mr->c = (Obj)ONE; NEXT(mr) = 0;
1.7 noro 581: MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
582: return;
583: }
584: w = (P *)ALLOCA(n*sizeof(P));
585: c = (Q *)ALLOCA(n*sizeof(Q));
586: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
587: if ( NUM(m->c) ) {
588: c[i] = (Q)m->c; w[i] = (P)ONE;
589: } else
1.66 ! noro 590: ptozp((P)m->c,1,&c[i],&w[i]);
1.7 noro 591: qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
592: if ( NUM(g) )
593: *rp = p;
594: else {
595: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1.66 ! noro 596: NEXTMP(mr0,mr); divsp(CO,(P)m->c,g,(P *)&mr->c); mr->dl = m->dl;
1.7 noro 597: }
598: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.46 noro 599: add_denomlist(g);
1.5 noro 600: }
1.7 noro 601: }
1.5 noro 602: }
603:
1.20 noro 604: void heu_nezgcdnpz(VL vl,P *pl,int m,P *pr)
1.5 noro 605: {
606: int i,r;
607: P gcd,t,s1,s2,u;
608: Q rq;
1.40 noro 609: DCP dc;
610: extern int DP_Print;
611:
1.5 noro 612: while ( 1 ) {
613: for ( i = 0, s1 = 0; i < m; i++ ) {
614: r = random(); UTOQ(r,rq);
615: mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
616: }
617: for ( i = 0, s2 = 0; i < m; i++ ) {
618: r = random(); UTOQ(r,rq);
619: mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
620: }
621: ezgcdp(vl,s1,s2,&gcd);
1.40 noro 622: if ( DP_Print > 2 )
623: { fprintf(asir_out,"(%d)",nmonop(gcd)); fflush(asir_out); }
1.5 noro 624: for ( i = 0; i < m; i++ ) {
625: if ( !divtpz(vl,pl[i],gcd,&t) )
626: break;
627: }
628: if ( i == m )
629: break;
630: }
631: *pr = gcd;
632: }
633:
1.20 noro 634: void dp_prim_mod(DP p,int mod,DP *rp)
1.5 noro 635: {
636: P t,g;
637: MP m,mr,mr0;
638:
639: if ( !p )
640: *rp = 0;
641: else if ( NoGCD )
642: *rp = p;
643: else {
1.66 ! noro 644: for ( m = BDY(p), g = (P)m->c, m = NEXT(m); m; m = NEXT(m) ) {
! 645: gcdprsmp(CO,mod,g,(P)m->c,&t); g = t;
1.5 noro 646: }
647: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1.66 ! noro 648: NEXTMP(mr0,mr); divsmp(CO,mod,(P)m->c,g,(P *)&mr->c); mr->dl = m->dl;
1.5 noro 649: }
650: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
651: }
652: }
653:
1.20 noro 654: void dp_cont(DP p,Q *rp)
1.5 noro 655: {
1.7 noro 656: VECT v;
1.5 noro 657:
1.7 noro 658: dp_dtov(p,&v); igcdv(v,rp);
1.5 noro 659: }
660:
1.20 noro 661: void dp_dtov(DP dp,VECT *rp)
1.5 noro 662: {
1.7 noro 663: MP m,t;
664: int i,n;
665: VECT v;
666: pointer *p;
1.5 noro 667:
1.7 noro 668: m = BDY(dp);
669: for ( t = m, n = 0; t; t = NEXT(t), n++ );
670: MKVECT(v,n);
671: for ( i = 0, p = BDY(v), t = m; i < n; t = NEXT(t), i++ )
672: p[i] = (pointer)(t->c);
673: *rp = v;
1.5 noro 674: }
675:
1.7 noro 676: /*
677: * s-poly computation
678: *
679: */
1.5 noro 680:
1.20 noro 681: void dp_sp(DP p1,DP p2,DP *rp)
1.5 noro 682: {
1.7 noro 683: int i,n,td;
684: int *w;
685: DL d1,d2,d;
686: MP m;
687: DP t,s1,s2,u;
688: Q c,c1,c2;
689: N gn,tn;
1.5 noro 690:
1.7 noro 691: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
692: w = (int *)ALLOCA(n*sizeof(int));
693: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 694: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.5 noro 695: }
1.7 noro 696:
697: NEWDL(d,n); d->td = td - d1->td;
698: for ( i = 0; i < n; i++ )
699: d->d[i] = w[i] - d1->d[i];
700: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
701: if ( INT(c1) && INT(c2) ) {
702: gcdn(NM(c1),NM(c2),&gn);
703: if ( !UNIN(gn) ) {
704: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
705: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1.5 noro 706: }
707: }
1.7 noro 708:
1.66 ! noro 709: NEWMP(m); m->dl = d; m->c = (Obj)c2; NEXT(m) = 0;
1.7 noro 710: MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
711:
712: NEWDL(d,n); d->td = td - d2->td;
713: for ( i = 0; i < n; i++ )
714: d->d[i] = w[i] - d2->d[i];
1.66 ! noro 715: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0;
1.7 noro 716: MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
717:
718: subd(CO,t,u,rp);
1.14 noro 719: if ( GenTrace ) {
720: LIST hist;
721: NODE node;
722:
1.58 noro 723: node = mknode(4,ONE,NULLP,s1,ONE);
1.14 noro 724: MKLIST(hist,node);
725: MKNODE(TraceList,hist,0);
726:
1.58 noro 727: node = mknode(4,ONE,NULLP,NULLP,ONE);
1.14 noro 728: chsgnd(s2,(DP *)&ARG2(node));
729: MKLIST(hist,node);
730: MKNODE(node,hist,TraceList); TraceList = node;
731: }
732: }
733:
1.66 ! noro 734: void dpm_sp(DPM p1,DPM p2,DPM *rp)
! 735: {
! 736: int i,n,td;
! 737: int *w;
! 738: DL d1,d2,d;
! 739: MP m;
! 740: DP s1,s2;
! 741: DPM t,u;
! 742: Q c,c1,c2;
! 743: N gn,tn;
! 744:
! 745: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 746: if ( BDY(p1)->pos != BDY(p2)->pos ) {
! 747: *rp = 0;
! 748: return;
! 749: }
! 750: w = (int *)ALLOCA(n*sizeof(int));
! 751: for ( i = 0, td = 0; i < n; i++ ) {
! 752: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
! 753: }
! 754:
! 755: NEWDL(d,n); d->td = td - d1->td;
! 756: for ( i = 0; i < n; i++ )
! 757: d->d[i] = w[i] - d1->d[i];
! 758: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 759: if ( INT(c1) && INT(c2) ) {
! 760: gcdn(NM(c1),NM(c2),&gn);
! 761: if ( !UNIN(gn) ) {
! 762: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 763: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 764: }
! 765: }
! 766:
! 767: NEWMP(m); m->dl = d; m->c = (Obj)c2; NEXT(m) = 0;
! 768: MKDP(n,m,s1); s1->sugar = d->td; mulobjdpm(CO,(Obj)s1,p1,&t);
! 769:
! 770: NEWDL(d,n); d->td = td - d2->td;
! 771: for ( i = 0; i < n; i++ )
! 772: d->d[i] = w[i] - d2->d[i];
! 773: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0;
! 774: MKDP(n,m,s2); s2->sugar = d->td; mulobjdpm(CO,(Obj)s2,p2,&u);
! 775:
! 776: subdpm(CO,t,u,rp);
! 777: if ( GenTrace ) {
! 778: LIST hist;
! 779: NODE node;
! 780:
! 781: node = mknode(4,ONE,NULLP,s1,ONE);
! 782: MKLIST(hist,node);
! 783: MKNODE(TraceList,hist,0);
! 784:
! 785: node = mknode(4,ONE,NULLP,NULLP,ONE);
! 786: chsgnd(s2,(DP *)&ARG2(node));
! 787: MKLIST(hist,node);
! 788: MKNODE(node,hist,TraceList); TraceList = node;
! 789: }
! 790: }
! 791:
1.20 noro 792: void _dp_sp_dup(DP p1,DP p2,DP *rp)
1.14 noro 793: {
794: int i,n,td;
795: int *w;
796: DL d1,d2,d;
797: MP m;
798: DP t,s1,s2,u;
799: Q c,c1,c2;
800: N gn,tn;
801:
802: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
803: w = (int *)ALLOCA(n*sizeof(int));
804: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 805: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.14 noro 806: }
807:
808: _NEWDL(d,n); d->td = td - d1->td;
809: for ( i = 0; i < n; i++ )
810: d->d[i] = w[i] - d1->d[i];
811: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
812: if ( INT(c1) && INT(c2) ) {
813: gcdn(NM(c1),NM(c2),&gn);
814: if ( !UNIN(gn) ) {
815: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
816: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
817: }
818: }
819:
1.66 ! noro 820: _NEWMP(m); m->dl = d; m->c = (Obj)c2; NEXT(m) = 0;
1.14 noro 821: _MKDP(n,m,s1); s1->sugar = d->td; _muld_dup(CO,s1,p1,&t); _free_dp(s1);
822:
823: _NEWDL(d,n); d->td = td - d2->td;
824: for ( i = 0; i < n; i++ )
825: d->d[i] = w[i] - d2->d[i];
1.66 ! noro 826: _NEWMP(m); m->dl = d; chsgnp((P)c1,(P *)&m->c); NEXT(m) = 0;
1.14 noro 827: _MKDP(n,m,s2); s2->sugar = d->td; _muld_dup(CO,s2,p2,&u); _free_dp(s2);
828:
829: _addd_destructive(CO,t,u,rp);
1.7 noro 830: if ( GenTrace ) {
831: LIST hist;
832: NODE node;
833:
1.58 noro 834: node = mknode(4,ONE,NULLP,s1,ONE);
1.7 noro 835: MKLIST(hist,node);
836: MKNODE(TraceList,hist,0);
837:
1.58 noro 838: node = mknode(4,ONE,NULLP,NULLP,ONE);
1.7 noro 839: chsgnd(s2,(DP *)&ARG2(node));
840: MKLIST(hist,node);
841: MKNODE(node,hist,TraceList); TraceList = node;
842: }
843: }
844:
1.20 noro 845: void dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7 noro 846: {
847: int i,n,td;
848: int *w;
849: DL d1,d2,d;
850: MP m;
851: DP t,s,u;
852:
853: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
854: w = (int *)ALLOCA(n*sizeof(int));
855: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 856: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7 noro 857: }
1.18 noro 858: NEWDL_NOINIT(d,n); d->td = td - d1->td;
1.7 noro 859: for ( i = 0; i < n; i++ )
860: d->d[i] = w[i] - d1->d[i];
1.66 ! noro 861: NEWMP(m); m->dl = d; m->c = (Obj)BDY(p2)->c; NEXT(m) = 0;
1.7 noro 862: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
1.18 noro 863: NEWDL_NOINIT(d,n); d->td = td - d2->td;
1.7 noro 864: for ( i = 0; i < n; i++ )
865: d->d[i] = w[i] - d2->d[i];
1.66 ! noro 866: NEWMP(m); m->dl = d; m->c = (Obj)BDY(p1)->c; NEXT(m) = 0;
1.7 noro 867: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
868: submd(CO,mod,t,u,rp);
869: }
870:
1.20 noro 871: void _dp_sp_mod_dup(DP p1,DP p2,int mod,DP *rp)
1.7 noro 872: {
873: int i,n,td;
874: int *w;
875: DL d1,d2,d;
876: MP m;
877: DP t,s,u;
878:
879: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
880: w = (int *)ALLOCA(n*sizeof(int));
881: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 882: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7 noro 883: }
884: _NEWDL(d,n); d->td = td - d1->td;
885: for ( i = 0; i < n; i++ )
886: d->d[i] = w[i] - d1->d[i];
887: _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
888: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
889: _NEWDL(d,n); d->td = td - d2->td;
890: for ( i = 0; i < n; i++ )
891: d->d[i] = w[i] - d2->d[i];
1.66 ! noro 892: _NEWMP(m); m->dl = d; m->c = (Obj)STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
1.7 noro 893: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
894: _addmd_destructive(mod,t,u,rp);
895: }
896:
1.20 noro 897: void _dp_sp_mod(DP p1,DP p2,int mod,DP *rp)
1.7 noro 898: {
899: int i,n,td;
900: int *w;
901: DL d1,d2,d;
902: MP m;
903: DP t,s,u;
904:
905: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
906: w = (int *)ALLOCA(n*sizeof(int));
907: for ( i = 0, td = 0; i < n; i++ ) {
1.21 noro 908: w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i);
1.7 noro 909: }
910: NEWDL(d,n); d->td = td - d1->td;
911: for ( i = 0; i < n; i++ )
912: d->d[i] = w[i] - d1->d[i];
913: NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
914: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p1,&t);
915: NEWDL(d,n); d->td = td - d2->td;
916: for ( i = 0; i < n; i++ )
917: d->d[i] = w[i] - d2->d[i];
1.66 ! noro 918: NEWMP(m); m->dl = d; m->c = (Obj)STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
1.7 noro 919: MKDP(n,m,s); s->sugar = d->td; mulmd_dup(mod,s,p2,&u);
920: addmd_destructive(mod,t,u,rp);
921: }
922:
923: /*
924: * m-reduction
1.13 noro 925: * do content reduction over Z or Q(x,...)
926: * do nothing over finite fields
1.7 noro 927: *
928: */
929:
1.20 noro 930: void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp,DP *multp)
1.7 noro 931: {
932: int i,n;
933: DL d1,d2,d;
934: MP m;
935: DP t,s,r,h;
936: Q c,c1,c2;
937: N gn,tn;
938: P g,a;
1.23 noro 939: P p[2];
1.7 noro 940:
941: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
942: NEWDL(d,n); d->td = d1->td - d2->td;
943: for ( i = 0; i < n; i++ )
944: d->d[i] = d1->d[i]-d2->d[i];
945: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
1.23 noro 946: if ( dp_fcoeffs == N_GFS ) {
947: p[0] = (P)c1; p[1] = (P)c2;
948: gcdsf(CO,p,2,&g);
949: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
950: } else if ( dp_fcoeffs ) {
1.7 noro 951: /* do nothing */
952: } else if ( INT(c1) && INT(c2) ) {
953: gcdn(NM(c1),NM(c2),&gn);
954: if ( !UNIN(gn) ) {
955: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
956: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
957: }
958: } else {
959: ezgcdpz(CO,(P)c1,(P)c2,&g);
960: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
1.46 noro 961: add_denomlist(g);
1.7 noro 962: }
1.66 ! noro 963: NEWMP(m); m->dl = d; chsgnp((P)c1,(P *)&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
1.7 noro 964: *multp = s;
1.66 ! noro 965: muld(CO,s,p2,&t); muldc(CO,p1,(Obj)c2,&s); addd(CO,s,t,&r);
! 966: muldc(CO,p0,(Obj)c2,&h);
1.7 noro 967: *head = h; *rest = r; *dnp = (P)c2;
968: }
969:
1.66 ! noro 970: void dpm_red(DPM p0,DPM p1,DPM p2,DPM *head,DPM *rest,P *dnp,DP *multp)
! 971: {
! 972: int i,n,pos;
! 973: DL d1,d2,d;
! 974: MP m;
! 975: DP s;
! 976: DPM t,r,h,u,w;
! 977: Q c,c1,c2;
! 978: N gn,tn;
! 979: P g,a;
! 980: P p[2];
! 981:
! 982: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; pos = BDY(p1)->pos;
! 983: if ( pos != BDY(p2)->pos )
! 984: error("dpm_red : cannot happen");
! 985: NEWDL(d,n); d->td = d1->td - d2->td;
! 986: for ( i = 0; i < n; i++ )
! 987: d->d[i] = d1->d[i]-d2->d[i];
! 988: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
! 989: if ( dp_fcoeffs == N_GFS ) {
! 990: p[0] = (P)c1; p[1] = (P)c2;
! 991: gcdsf(CO,p,2,&g);
! 992: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 993: } else if ( dp_fcoeffs ) {
! 994: /* do nothing */
! 995: } else if ( INT(c1) && INT(c2) ) {
! 996: gcdn(NM(c1),NM(c2),&gn);
! 997: if ( !UNIN(gn) ) {
! 998: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
! 999: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
! 1000: }
! 1001: } else {
! 1002: ezgcdpz(CO,(P)c1,(P)c2,&g);
! 1003: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
! 1004: add_denomlist(g);
! 1005: }
! 1006: NEWMP(m); m->dl = d; chsgnp((P)c1,(P *)&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
! 1007: *multp = s;
! 1008: mulobjdpm(CO,(Obj)s,p2,&u); mulobjdpm(CO,(Obj)c2,p1,&w); adddpm(CO,u,w,&r);
! 1009: mulobjdpm(CO,(Obj)c2,p0,&h);
! 1010: *head = h; *rest = r; *dnp = (P)c2;
! 1011: }
! 1012:
! 1013:
1.41 noro 1014: /*
1015: * m-reduction by a marked poly
1016: * do content reduction over Z or Q(x,...)
1017: * do nothing over finite fields
1018: *
1019: */
1020:
1021:
1022: void dp_red_marked(DP p0,DP p1,DP p2,DP hp2,DP *head,DP *rest,P *dnp,DP *multp)
1023: {
1024: int i,n;
1025: DL d1,d2,d;
1026: MP m;
1027: DP t,s,r,h;
1028: Q c,c1,c2;
1029: N gn,tn;
1030: P g,a;
1031: P p[2];
1032:
1033: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl;
1034: NEWDL(d,n); d->td = d1->td - d2->td;
1035: for ( i = 0; i < n; i++ )
1036: d->d[i] = d1->d[i]-d2->d[i];
1037: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(hp2)->c;
1038: if ( dp_fcoeffs == N_GFS ) {
1039: p[0] = (P)c1; p[1] = (P)c2;
1040: gcdsf(CO,p,2,&g);
1041: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
1042: } else if ( dp_fcoeffs ) {
1043: /* do nothing */
1044: } else if ( INT(c1) && INT(c2) ) {
1045: gcdn(NM(c1),NM(c2),&gn);
1046: if ( !UNIN(gn) ) {
1047: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
1048: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1049: }
1050: } else {
1051: ezgcdpz(CO,(P)c1,(P)c2,&g);
1052: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
1053: }
1.66 ! noro 1054: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
1.41 noro 1055: *multp = s;
1.66 ! noro 1056: muld(CO,s,p2,&t); muldc(CO,p1,(Obj)c2,&s); subd(CO,s,t,&r);
! 1057: muldc(CO,p0,(Obj)c2,&h);
1.41 noro 1058: *head = h; *rest = r; *dnp = (P)c2;
1059: }
1060:
1.55 noro 1061: void dp_red_marked_mod(DP p0,DP p1,DP p2,DP hp2,int mod,DP *head,DP *rest,P *dnp,DP *multp)
1.44 noro 1062: {
1063: int i,n;
1064: DL d1,d2,d;
1065: MP m;
1066: DP t,s,r,h;
1067: P c1,c2,g,u;
1068:
1069: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl;
1070: NEWDL(d,n); d->td = d1->td - d2->td;
1071: for ( i = 0; i < n; i++ )
1072: d->d[i] = d1->d[i]-d2->d[i];
1073: c1 = (P)BDY(p1)->c; c2 = (P)BDY(hp2)->c;
1074: gcdprsmp(CO,mod,c1,c2,&g);
1075: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
1076: if ( NUM(c2) ) {
1077: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
1078: }
1.66 ! noro 1079: NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0;
1.55 noro 1080: MKDP(n,m,s); s->sugar = d->td;
1081: *multp = s;
1082: mulmd(CO,mod,s,p2,&t);
1.44 noro 1083: if ( NUM(c2) ) {
1.55 noro 1084: submd(CO,mod,p1,t,&r); h = p0;
1.44 noro 1085: } else {
1.55 noro 1086: mulmdc(CO,mod,p1,c2,&s); submd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
1.44 noro 1087: }
1088: *head = h; *rest = r; *dnp = c2;
1089: }
1090:
1.13 noro 1091: /* m-reduction over a field */
1092:
1.20 noro 1093: void dp_red_f(DP p1,DP p2,DP *rest)
1.13 noro 1094: {
1095: int i,n;
1096: DL d1,d2,d;
1097: MP m;
1.20 noro 1098: DP t,s;
1.13 noro 1099: Obj a,b;
1100:
1101: n = p1->nv;
1102: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1103:
1104: NEWDL(d,n); d->td = d1->td - d2->td;
1105: for ( i = 0; i < n; i++ )
1106: d->d[i] = d1->d[i]-d2->d[i];
1107:
1108: NEWMP(m); m->dl = d;
1109: divr(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); chsgnr(a,&b);
1.66 ! noro 1110: C(m) = (Obj)b;
1.13 noro 1111: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
1112:
1113: muld(CO,s,p2,&t); addd(CO,p1,t,rest);
1114: }
1115:
1.66 ! noro 1116: void dpm_red_f(DPM p1,DPM p2,DPM *rest)
! 1117: {
! 1118: int i,n;
! 1119: DL d1,d2,d;
! 1120: MP m;
! 1121: DPM t;
! 1122: DP s;
! 1123: Obj a,b;
! 1124:
! 1125: n = p1->nv;
! 1126: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 1127:
! 1128: NEWDL(d,n); d->td = d1->td - d2->td;
! 1129: for ( i = 0; i < n; i++ )
! 1130: d->d[i] = d1->d[i]-d2->d[i];
! 1131:
! 1132: NEWMP(m); m->dl = d;
! 1133: arf_div(CO,(Obj)BDY(p1)->c,(Obj)BDY(p2)->c,&a); arf_chsgn(a,&b);
! 1134: C(m) = b;
! 1135: NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
! 1136:
! 1137: mulobjdpm(CO,(Obj)s,p2,&t); adddpm(CO,p1,t,rest);
! 1138: }
! 1139:
! 1140:
1.20 noro 1141: void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp)
1.7 noro 1142: {
1143: int i,n;
1144: DL d1,d2,d;
1145: MP m;
1146: DP t,s,r,h;
1147: P c1,c2,g,u;
1148:
1149: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1150: NEWDL(d,n); d->td = d1->td - d2->td;
1151: for ( i = 0; i < n; i++ )
1152: d->d[i] = d1->d[i]-d2->d[i];
1153: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
1154: gcdprsmp(CO,mod,c1,c2,&g);
1155: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
1156: if ( NUM(c2) ) {
1157: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
1158: }
1.66 ! noro 1159: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,(P *)&m->c); NEXT(m) = 0;
1.11 noro 1160: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t);
1.7 noro 1161: if ( NUM(c2) ) {
1162: addmd(CO,mod,p1,t,&r); h = p0;
1163: } else {
1164: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
1165: }
1166: *head = h; *rest = r; *dnp = c2;
1167: }
1168:
1.10 noro 1169: struct oEGT eg_red_mod;
1170:
1.20 noro 1171: void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp)
1.7 noro 1172: {
1173: int i,n;
1174: DL d1,d2,d;
1175: MP m;
1176: DP t,s;
1.16 noro 1177: int c,c1,c2;
1178: extern int do_weyl;
1.7 noro 1179:
1180: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1181: _NEWDL(d,n); d->td = d1->td - d2->td;
1182: for ( i = 0; i < n; i++ )
1183: d->d[i] = d1->d[i]-d2->d[i];
1.16 noro 1184: c = invm(ITOS(BDY(p2)->c),mod);
1185: c2 = ITOS(BDY(p1)->c);
1186: DMAR(c,c2,0,mod,c1);
1.66 ! noro 1187: _NEWMP(m); m->dl = d; m->c = (Obj)STOI(mod-c1); NEXT(m) = 0;
1.16 noro 1188: #if 0
1.7 noro 1189: _MKDP(n,m,s); s->sugar = d->td;
1190: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 1191: #else
1192: if ( do_weyl ) {
1.19 noro 1193: _MKDP(n,m,s); s->sugar = d->td;
1194: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
1.16 noro 1195: } else {
1196: _mulmdm_dup(mod,p2,m,&t); _FREEMP(m);
1197: }
1198: #endif
1.10 noro 1199: /* get_eg(&t0); */
1.7 noro 1200: _addmd_destructive(mod,p1,t,rp);
1.10 noro 1201: /* get_eg(&t1); add_eg(&eg_red_mod,&t0,&t1); */
1.7 noro 1202: }
1203:
1204: /*
1205: * normal form computation
1206: *
1207: */
1.5 noro 1208:
1.20 noro 1209: void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp)
1.5 noro 1210: {
1211: DP u,p,d,s,t,dmy;
1212: NODE l;
1213: MP m,mr;
1214: int i,n;
1215: int *wb;
1216: int sugar,psugar;
1217: P dn,tdn,tdn1;
1218:
1219: dn = (P)ONE;
1220: if ( !g ) {
1221: *rp = 0; *dnp = dn; return;
1222: }
1223: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1224: wb = (int *)ALLOCA(n*sizeof(int));
1225: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1226: wb[i] = QTOS((Q)BDY(l));
1227: sugar = g->sugar;
1228: for ( d = 0; g; ) {
1229: for ( u = 0, i = 0; i < n; i++ ) {
1230: if ( dp_redble(g,p = ps[wb[i]]) ) {
1231: dp_red(d,g,p,&t,&u,&tdn,&dmy);
1232: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1233: sugar = MAX(sugar,psugar);
1234: if ( !u ) {
1235: if ( d )
1236: d->sugar = sugar;
1237: *rp = d; *dnp = dn; return;
1238: } else {
1239: d = t;
1240: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
1241: }
1242: break;
1243: }
1244: }
1245: if ( u )
1246: g = u;
1247: else if ( !full ) {
1248: if ( g ) {
1249: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1250: }
1251: *rp = g; *dnp = dn; return;
1252: } else {
1253: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1254: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1255: addd(CO,d,t,&s); d = s;
1256: dp_rest(g,&t); g = t;
1257: }
1258: }
1259: if ( d )
1260: d->sugar = sugar;
1261: *rp = d; *dnp = dn;
1262: }
1263:
1.43 noro 1264: void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Q *contp)
1265: {
1266: struct oVECT v;
1267: int i,n1,n2,n;
1268: MP m,m0,t;
1269: Q *w;
1270: Q h;
1271:
1272: if ( p1 ) {
1273: for ( i = 0, m = BDY(p1); m; m = NEXT(m), i++ );
1274: n1 = i;
1275: } else
1276: n1 = 0;
1277: if ( p2 ) {
1278: for ( i = 0, m = BDY(p2); m; m = NEXT(m), i++ );
1279: n2 = i;
1280: } else
1281: n2 = 0;
1282: n = n1+n2;
1283: if ( !n ) {
1284: *r1p = 0; *r2p = 0; *contp = ONE; return;
1285: }
1286: w = (Q *)ALLOCA(n*sizeof(Q));
1287: v.len = n;
1288: v.body = (pointer *)w;
1289: i = 0;
1290: if ( p1 )
1291: for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) w[i] = (Q)m->c;
1292: if ( p2 )
1293: for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = (Q)m->c;
1294: h = w[0]; removecont_array((P *)w,n,1); divq(h,w[0],contp);
1295: i = 0;
1296: if ( p1 ) {
1297: for ( m0 = 0, t = BDY(p1); i < n1; i++, t = NEXT(t) ) {
1.66 ! noro 1298: NEXTMP(m0,m); m->c = (Obj)w[i]; m->dl = t->dl;
1.43 noro 1299: }
1300: NEXT(m) = 0;
1301: MKDP(p1->nv,m0,*r1p); (*r1p)->sugar = p1->sugar;
1302: } else
1303: *r1p = 0;
1304: if ( p2 ) {
1305: for ( m0 = 0, t = BDY(p2); i < n; i++, t = NEXT(t) ) {
1.66 ! noro 1306: NEXTMP(m0,m); m->c = (Obj)w[i]; m->dl = t->dl;
1.43 noro 1307: }
1308: NEXT(m) = 0;
1309: MKDP(p2->nv,m0,*r2p); (*r2p)->sugar = p2->sugar;
1310: } else
1311: *r2p = 0;
1312: }
1313:
1.41 noro 1314: /* true nf by a marked GB */
1315:
1.43 noro 1316: void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP *rp,P *nmp,P *dnp)
1.41 noro 1317: {
1318: DP u,p,d,s,t,dmy,hp;
1319: NODE l;
1320: MP m,mr;
1.43 noro 1321: int i,n,hmag;
1.41 noro 1322: int *wb;
1.43 noro 1323: int sugar,psugar,multiple;
1324: P nm,tnm1,dn,tdn,tdn1;
1325: Q cont;
1.41 noro 1326:
1.43 noro 1327: multiple = 0;
1328: hmag = multiple*HMAG(g);
1329: nm = (P)ONE;
1.41 noro 1330: dn = (P)ONE;
1331: if ( !g ) {
1332: *rp = 0; *dnp = dn; return;
1333: }
1334: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1335: wb = (int *)ALLOCA(n*sizeof(int));
1336: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1337: wb[i] = QTOS((Q)BDY(l));
1338: sugar = g->sugar;
1339: for ( d = 0; g; ) {
1340: for ( u = 0, i = 0; i < n; i++ ) {
1341: if ( dp_redble(g,hp = hps[wb[i]]) ) {
1342: p = ps[wb[i]];
1343: dp_red_marked(d,g,p,hp,&t,&u,&tdn,&dmy);
1344: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1345: sugar = MAX(sugar,psugar);
1346: if ( !u ) {
1.43 noro 1347: goto last;
1.41 noro 1348: } else {
1349: d = t;
1350: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
1351: }
1352: break;
1353: }
1354: }
1.43 noro 1355: if ( u ) {
1.41 noro 1356: g = u;
1.43 noro 1357: if ( multiple && ((d && HMAG(d)>hmag) || (HMAG(g)>hmag)) ) {
1358: dp_removecont2(d,g,&t,&u,&cont); d = t; g = u;
1359: mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
1360: if ( d )
1361: hmag = multiple*HMAG(d);
1362: else
1363: hmag = multiple*HMAG(g);
1364: }
1365: } else {
1.41 noro 1366: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1367: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1368: addd(CO,d,t,&s); d = s;
1369: dp_rest(g,&t); g = t;
1370: }
1371: }
1.43 noro 1372: last:
1373: if ( d ) {
1374: dp_removecont2(d,0,&t,&u,&cont); d = t;
1375: mulp(CO,nm,(P)cont,&tnm1); nm = tnm1;
1.41 noro 1376: d->sugar = sugar;
1.43 noro 1377: }
1378: *rp = d; *nmp = nm; *dnp = dn;
1.41 noro 1379: }
1380:
1.44 noro 1381: void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
1382: {
1.55 noro 1383: DP hp,u,p,d,s,t,dmy;
1.44 noro 1384: NODE l;
1385: MP m,mr;
1386: int i,n;
1387: int *wb;
1388: int sugar,psugar;
1389: P dn,tdn,tdn1;
1390:
1391: dn = (P)ONEM;
1392: if ( !g ) {
1393: *rp = 0; *dnp = dn; return;
1394: }
1395: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1396: wb = (int *)ALLOCA(n*sizeof(int));
1397: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1398: wb[i] = QTOS((Q)BDY(l));
1399: sugar = g->sugar;
1400: for ( d = 0; g; ) {
1401: for ( u = 0, i = 0; i < n; i++ ) {
1402: if ( dp_redble(g,hp = hps[wb[i]]) ) {
1403: p = ps[wb[i]];
1.55 noro 1404: dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&dmy);
1.44 noro 1405: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1406: sugar = MAX(sugar,psugar);
1407: if ( !u ) {
1408: if ( d )
1409: d->sugar = sugar;
1410: *rp = d; *dnp = dn; return;
1411: } else {
1412: d = t;
1413: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1414: }
1415: break;
1416: }
1417: }
1418: if ( u )
1419: g = u;
1420: else {
1421: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1422: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1423: addmd(CO,mod,d,t,&s); d = s;
1424: dp_rest(g,&t); g = t;
1425: }
1426: }
1427: if ( d )
1428: d->sugar = sugar;
1429: *rp = d; *dnp = dn;
1430: }
1431:
1.47 noro 1432: /* true nf by a marked GB and collect quotients */
1433:
1434: DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps,DP *hps,DP *rp,P *dnp)
1435: {
1436: DP u,p,d,s,t,dmy,hp,mult;
1437: DP *q;
1438: NODE l;
1439: MP m,mr;
1440: int i,n,j;
1441: int *wb;
1442: int sugar,psugar,multiple;
1443: P nm,tnm1,dn,tdn,tdn1;
1444: Q cont;
1445:
1446: dn = (P)ONE;
1447: if ( !g ) {
1.59 noro 1448: *rp = 0; *dnp = dn; return 0;
1.47 noro 1449: }
1450: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1451: wb = (int *)ALLOCA(n*sizeof(int));
1452: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1453: wb[i] = QTOS((Q)BDY(l));
1454: q = (DP *)MALLOC(n*sizeof(DP));
1455: for ( i = 0; i < n; i++ ) q[i] = 0;
1456: sugar = g->sugar;
1457: for ( d = 0; g; ) {
1458: for ( u = 0, i = 0; i < n; i++ ) {
1459: if ( dp_redble(g,hp = hps[wb[i]]) ) {
1460: p = ps[wb[i]];
1461: dp_red_marked(d,g,p,hp,&t,&u,&tdn,&mult);
1462: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1463: sugar = MAX(sugar,psugar);
1464: for ( j = 0; j < n; j++ ) {
1.66 ! noro 1465: muldc(CO,q[j],(Obj)tdn,&dmy); q[j] = dmy;
1.47 noro 1466: }
1467: addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
1468: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
1469: d = t;
1470: if ( !u ) goto last;
1471: break;
1472: }
1473: }
1474: if ( u ) {
1475: g = u;
1476: } else {
1477: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1478: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1479: addd(CO,d,t,&s); d = s;
1480: dp_rest(g,&t); g = t;
1481: }
1482: }
1483: last:
1484: if ( d ) d->sugar = sugar;
1485: *rp = d; *dnp = dn;
1486: return q;
1487: }
1488:
1.55 noro 1489: DP *dp_true_nf_and_quotient_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp)
1490: {
1491: DP u,p,d,s,t,dmy,hp,mult;
1492: DP *q;
1493: NODE l;
1494: MP m,mr;
1495: int i,n,j;
1496: int *wb;
1497: int sugar,psugar;
1498: P dn,tdn,tdn1;
1499:
1500: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1501: q = (DP *)MALLOC(n*sizeof(DP));
1502: for ( i = 0; i < n; i++ ) q[i] = 0;
1503: dn = (P)ONEM;
1504: if ( !g ) {
1.59 noro 1505: *rp = 0; *dnp = dn; return 0;
1.55 noro 1506: }
1507: wb = (int *)ALLOCA(n*sizeof(int));
1508: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1509: wb[i] = QTOS((Q)BDY(l));
1510: sugar = g->sugar;
1511: for ( d = 0; g; ) {
1512: for ( u = 0, i = 0; i < n; i++ ) {
1513: if ( dp_redble(g,hp = hps[wb[i]]) ) {
1514: p = ps[wb[i]];
1515: dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&mult);
1516: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1517: sugar = MAX(sugar,psugar);
1518: for ( j = 0; j < n; j++ ) {
1519: mulmdc(CO,mod,q[j],(P)tdn,&dmy); q[j] = dmy;
1520: }
1521: addmd(CO,mod,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
1522: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1523: d = t;
1524: if ( !u ) goto last;
1525: break;
1526: }
1527: }
1528: if ( u )
1529: g = u;
1530: else {
1531: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1532: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1533: addmd(CO,mod,d,t,&s); d = s;
1534: dp_rest(g,&t); g = t;
1535: }
1536: }
1537: last:
1538: if ( d )
1539: d->sugar = sugar;
1540: *rp = d; *dnp = dn;
1541: return q;
1542: }
1543:
1.13 noro 1544: /* nf computation over Z */
1545:
1.20 noro 1546: void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
1.5 noro 1547: {
1548: DP u,p,d,s,t,dmy1;
1549: P dmy;
1550: NODE l;
1551: MP m,mr;
1552: int i,n;
1553: int *wb;
1554: int hmag;
1555: int sugar,psugar;
1556:
1557: if ( !g ) {
1558: *rp = 0; return;
1559: }
1560: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1561: wb = (int *)ALLOCA(n*sizeof(int));
1562: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1563: wb[i] = QTOS((Q)BDY(l));
1.12 noro 1564:
1.13 noro 1565: hmag = multiple*HMAG(g);
1.5 noro 1566: sugar = g->sugar;
1.12 noro 1567:
1.5 noro 1568: for ( d = 0; g; ) {
1569: for ( u = 0, i = 0; i < n; i++ ) {
1570: if ( dp_redble(g,p = ps[wb[i]]) ) {
1571: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
1572: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1573: sugar = MAX(sugar,psugar);
1574: if ( !u ) {
1575: if ( d )
1576: d->sugar = sugar;
1577: *rp = d; return;
1578: }
1579: d = t;
1580: break;
1581: }
1582: }
1583: if ( u ) {
1584: g = u;
1585: if ( d ) {
1.13 noro 1586: if ( multiple && HMAG(d) > hmag ) {
1.5 noro 1587: dp_ptozp2(d,g,&t,&u); d = t; g = u;
1588: hmag = multiple*HMAG(d);
1589: }
1590: } else {
1.13 noro 1591: if ( multiple && HMAG(g) > hmag ) {
1.5 noro 1592: dp_ptozp(g,&t); g = t;
1593: hmag = multiple*HMAG(g);
1594: }
1595: }
1596: }
1597: else if ( !full ) {
1598: if ( g ) {
1599: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1600: }
1601: *rp = g; return;
1602: } else {
1603: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1604: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1605: addd(CO,d,t,&s); d = s;
1606: dp_rest(g,&t); g = t;
1607:
1608: }
1609: }
1610: if ( d )
1611: d->sugar = sugar;
1612: *rp = d;
1613: }
1614:
1.66 ! noro 1615: void dpm_nf_z(NODE b,DPM g,DPM *ps,int full,int multiple,DPM *rp)
! 1616: {
! 1617: DPM u,p,d,s,t;
! 1618: DP dmy1;
! 1619: P dmy;
! 1620: NODE l;
! 1621: DMM m,mr;
! 1622: int i,n;
! 1623: int *wb;
! 1624: int hmag;
! 1625: int sugar,psugar;
! 1626:
! 1627: if ( !g ) {
! 1628: *rp = 0; return;
! 1629: }
! 1630: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1631: wb = (int *)ALLOCA(n*sizeof(int));
! 1632: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1633: wb[i] = QTOS((Q)BDY(l));
! 1634:
! 1635: hmag = multiple*HMAG(g);
! 1636: sugar = g->sugar;
! 1637:
! 1638: for ( d = 0; g; ) {
! 1639: for ( u = 0, i = 0; i < n; i++ ) {
! 1640: if ( dpm_redble(g,p = ps[wb[i]]) ) {
! 1641: dpm_red(d,g,p,&t,&u,&dmy,&dmy1);
! 1642: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1643: sugar = MAX(sugar,psugar);
! 1644: if ( !u ) {
! 1645: if ( d )
! 1646: d->sugar = sugar;
! 1647: *rp = d; return;
! 1648: }
! 1649: d = t;
! 1650: break;
! 1651: }
! 1652: }
! 1653: if ( u ) {
! 1654: g = u;
! 1655: if ( d ) {
! 1656: if ( multiple && HMAG(d) > hmag ) {
! 1657: dpm_ptozp2(d,g,&t,&u); d = t; g = u;
! 1658: hmag = multiple*HMAG(d);
! 1659: }
! 1660: } else {
! 1661: if ( multiple && HMAG(g) > hmag ) {
! 1662: dpm_ptozp(g,&t); g = t;
! 1663: hmag = multiple*HMAG(g);
! 1664: }
! 1665: }
! 1666: }
! 1667: else if ( !full ) {
! 1668: if ( g ) {
! 1669: MKDPM(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1670: }
! 1671: *rp = g; return;
! 1672: } else {
! 1673: m = BDY(g); NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos;
! 1674: NEXT(mr) = 0; MKDPM(g->nv,mr,t); t->sugar = mr->dl->td;
! 1675: adddpm(CO,d,t,&s); d = s;
! 1676: dpm_rest(g,&t); g = t;
! 1677: }
! 1678: }
! 1679: if ( d )
! 1680: d->sugar = sugar;
! 1681: *rp = d;
! 1682: }
! 1683:
1.13 noro 1684: /* nf computation over a field */
1685:
1.20 noro 1686: void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp)
1.13 noro 1687: {
1688: DP u,p,d,s,t;
1689: NODE l;
1690: MP m,mr;
1691: int i,n;
1692: int *wb;
1693: int sugar,psugar;
1694:
1695: if ( !g ) {
1696: *rp = 0; return;
1697: }
1698: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1699: wb = (int *)ALLOCA(n*sizeof(int));
1700: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1701: wb[i] = QTOS((Q)BDY(l));
1702:
1703: sugar = g->sugar;
1704: for ( d = 0; g; ) {
1705: for ( u = 0, i = 0; i < n; i++ ) {
1706: if ( dp_redble(g,p = ps[wb[i]]) ) {
1707: dp_red_f(g,p,&u);
1708: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1709: sugar = MAX(sugar,psugar);
1710: if ( !u ) {
1711: if ( d )
1712: d->sugar = sugar;
1713: *rp = d; return;
1714: }
1715: break;
1716: }
1717: }
1718: if ( u )
1719: g = u;
1720: else if ( !full ) {
1721: if ( g ) {
1722: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1723: }
1724: *rp = g; return;
1725: } else {
1726: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1727: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1728: addd(CO,d,t,&s); d = s;
1729: dp_rest(g,&t); g = t;
1730: }
1731: }
1732: if ( d )
1733: d->sugar = sugar;
1734: *rp = d;
1735: }
1736:
1.66 ! noro 1737: void dpm_nf_f(NODE b,DPM g,DPM *ps,int full,DPM *rp)
! 1738: {
! 1739: DPM u,p,d,s,t;
! 1740: NODE l;
! 1741: DMM m,mr;
! 1742: int i,n;
! 1743: int *wb;
! 1744: int sugar,psugar;
! 1745:
! 1746: if ( !g ) {
! 1747: *rp = 0; return;
! 1748: }
! 1749: for ( n = 0, l = b; l; l = NEXT(l), n++ );
! 1750: wb = (int *)ALLOCA(n*sizeof(int));
! 1751: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
! 1752: wb[i] = QTOS((Q)BDY(l));
! 1753:
! 1754: sugar = g->sugar;
! 1755: for ( d = 0; g; ) {
! 1756: for ( u = 0, i = 0; i < n; i++ ) {
! 1757: if ( dpm_redble(g,p = ps[wb[i]]) ) {
! 1758: dpm_red_f(g,p,&u);
! 1759: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
! 1760: sugar = MAX(sugar,psugar);
! 1761: if ( !u ) {
! 1762: if ( d )
! 1763: d->sugar = sugar;
! 1764: *rp = d; return;
! 1765: }
! 1766: break;
! 1767: }
! 1768: }
! 1769: if ( u )
! 1770: g = u;
! 1771: else if ( !full ) {
! 1772: if ( g ) {
! 1773: MKDPM(g->nv,BDY(g),t); t->sugar = sugar; g = t;
! 1774: }
! 1775: *rp = g; return;
! 1776: } else {
! 1777: m = BDY(g); NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos;
! 1778: NEXT(mr) = 0; MKDPM(g->nv,mr,t); t->sugar = mr->dl->td;
! 1779: adddpm(CO,d,t,&s); d = s;
! 1780: dpm_rest(g,&t); g = t;
! 1781: }
! 1782: }
! 1783: if ( d )
! 1784: d->sugar = sugar;
! 1785: *rp = d;
! 1786: }
! 1787:
1.13 noro 1788: /* nf computation over GF(mod) (only for internal use) */
1789:
1.20 noro 1790: void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1791: {
1792: DP u,p,d,s,t;
1793: P dmy;
1794: NODE l;
1795: MP m,mr;
1796: int sugar,psugar;
1797:
1798: if ( !g ) {
1799: *rp = 0; return;
1800: }
1801: sugar = g->sugar;
1802: for ( d = 0; g; ) {
1803: for ( u = 0, l = b; l; l = NEXT(l) ) {
1804: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1805: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
1806: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1807: sugar = MAX(sugar,psugar);
1808: if ( !u ) {
1809: if ( d )
1810: d->sugar = sugar;
1811: *rp = d; return;
1812: }
1813: d = t;
1814: break;
1815: }
1816: }
1817: if ( u )
1818: g = u;
1819: else if ( !full ) {
1820: if ( g ) {
1821: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1822: }
1823: *rp = g; return;
1824: } else {
1825: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1826: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1827: addmd(CO,mod,d,t,&s); d = s;
1828: dp_rest(g,&t); g = t;
1829: }
1830: }
1831: if ( d )
1832: d->sugar = sugar;
1833: *rp = d;
1834: }
1835:
1.20 noro 1836: void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp)
1.5 noro 1837: {
1838: DP u,p,d,s,t;
1839: NODE l;
1840: MP m,mr;
1841: int i,n;
1842: int *wb;
1843: int sugar,psugar;
1844: P dn,tdn,tdn1;
1845:
1846: dn = (P)ONEM;
1847: if ( !g ) {
1848: *rp = 0; *dnp = dn; return;
1849: }
1850: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1851: wb = (int *)ALLOCA(n*sizeof(int));
1852: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1853: wb[i] = QTOS((Q)BDY(l));
1854: sugar = g->sugar;
1855: for ( d = 0; g; ) {
1856: for ( u = 0, i = 0; i < n; i++ ) {
1857: if ( dp_redble(g,p = ps[wb[i]]) ) {
1858: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
1859: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1860: sugar = MAX(sugar,psugar);
1861: if ( !u ) {
1862: if ( d )
1863: d->sugar = sugar;
1864: *rp = d; *dnp = dn; return;
1865: } else {
1866: d = t;
1867: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1868: }
1869: break;
1870: }
1871: }
1872: if ( u )
1873: g = u;
1874: else if ( !full ) {
1875: if ( g ) {
1876: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1877: }
1878: *rp = g; *dnp = dn; return;
1879: } else {
1880: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1881: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1882: addmd(CO,mod,d,t,&s); d = s;
1883: dp_rest(g,&t); g = t;
1884: }
1885: }
1886: if ( d )
1887: d->sugar = sugar;
1888: *rp = d; *dnp = dn;
1889: }
1890:
1.20 noro 1891: void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp)
1.5 noro 1892: {
1.20 noro 1893: DP u,p,d;
1.7 noro 1894: NODE l;
1.20 noro 1895: MP m,mrd;
1896: int sugar,psugar,n,h_reducible;
1.5 noro 1897:
1.7 noro 1898: if ( !g ) {
1899: *rp = 0; return;
1.5 noro 1900: }
1.7 noro 1901: sugar = g->sugar;
1902: n = g->nv;
1903: for ( d = 0; g; ) {
1904: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
1905: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
1906: h_reducible = 1;
1907: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1908: _dp_red_mod_destructive(g,p,mod,&u); g = u;
1909: sugar = MAX(sugar,psugar);
1910: if ( !g ) {
1911: if ( d )
1912: d->sugar = sugar;
1913: _dptodp(d,rp); _free_dp(d); return;
1914: }
1915: break;
1916: }
1917: }
1918: if ( !h_reducible ) {
1919: /* head term is not reducible */
1920: if ( !full ) {
1921: if ( g )
1922: g->sugar = sugar;
1923: _dptodp(g,rp); _free_dp(g); return;
1924: } else {
1925: m = BDY(g);
1926: if ( NEXT(m) ) {
1927: BDY(g) = NEXT(m); NEXT(m) = 0;
1928: } else {
1929: _FREEDP(g); g = 0;
1930: }
1931: if ( d ) {
1932: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
1933: NEXT(mrd) = m;
1934: } else {
1935: _MKDP(n,m,d);
1936: }
1937: }
1938: }
1.5 noro 1939: }
1.7 noro 1940: if ( d )
1941: d->sugar = sugar;
1942: _dptodp(d,rp); _free_dp(d);
1.5 noro 1943: }
1.13 noro 1944:
1945: /* reduction by linear base over a field */
1946:
1.20 noro 1947: void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p)
1.13 noro 1948: {
1949: DP r1,r2,b1,b2,t,s;
1950: Obj c,c1,c2;
1951: NODE l,b;
1952: int n;
1953:
1954: if ( !p1 ) {
1955: *r1p = p1; *r2p = p2; return;
1956: }
1957: n = p1->nv;
1958: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1959: if ( !r1 ) {
1960: *r1p = r1; *r2p = r2; return;
1961: }
1962: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1963: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1964: b2 = (DP)BDY(NEXT(b));
1965: divr(CO,(Obj)ONE,(Obj)BDY(b1)->c,&c1);
1966: mulr(CO,c1,(Obj)BDY(r1)->c,&c2); chsgnr(c2,&c);
1.66 ! noro 1967: muldc(CO,b1,(Obj)c,&t); addd(CO,r1,t,&s); r1 = s;
! 1968: muldc(CO,b2,(Obj)c,&t); addd(CO,r2,t,&s); r2 = s;
1.13 noro 1969: }
1970: }
1971: *r1p = r1; *r2p = r2;
1972: }
1973:
1974: /* reduction by linear base over GF(mod) */
1.5 noro 1975:
1.20 noro 1976: void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p)
1.5 noro 1977: {
1.7 noro 1978: DP r1,r2,b1,b2,t,s;
1979: P c;
1980: MQ c1,c2;
1981: NODE l,b;
1982: int n;
1983:
1984: if ( !p1 ) {
1985: *r1p = p1; *r2p = p2; return;
1986: }
1987: n = p1->nv;
1988: for ( l = g, r1 = p1, r2 = p2; l; l = NEXT(l) ) {
1989: if ( !r1 ) {
1990: *r1p = r1; *r2p = r2; return;
1991: }
1992: b = BDY((LIST)BDY(l)); b1 = (DP)BDY(b);
1993: if ( dl_equal(n,BDY(r1)->dl,BDY(b1)->dl) ) {
1994: b2 = (DP)BDY(NEXT(b));
1995: invmq(mod,(MQ)BDY(b1)->c,&c1);
1996: mulmq(mod,c1,(MQ)BDY(r1)->c,&c2); chsgnmp(mod,(P)c2,&c);
1997: mulmdc(CO,mod,b1,c,&t); addmd(CO,mod,r1,t,&s); r1 = s;
1998: mulmdc(CO,mod,b2,c,&t); addmd(CO,mod,r2,t,&s); r2 = s;
1999: }
2000: }
2001: *r1p = r1; *r2p = r2;
1.5 noro 2002: }
2003:
1.20 noro 2004: void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp)
1.5 noro 2005: {
1.7 noro 2006: DP s,t,u;
2007: MP m;
2008: DL h;
2009: int i,n;
2010:
2011: if ( !p ) {
2012: *rp = p; return;
2013: }
2014: n = p->nv;
2015: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
2016: h = m->dl;
2017: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
2018: i++;
1.66 ! noro 2019: mulmdc(CO,mod,(DP)BDY(NEXT(BDY(tab[i]))),(P)m->c,&t);
1.7 noro 2020: addmd(CO,mod,s,t,&u); s = u;
1.24 noro 2021: }
2022: *rp = s;
2023: }
2024:
2025: void dp_nf_tab_f(DP p,LIST *tab,DP *rp)
2026: {
2027: DP s,t,u;
2028: MP m;
2029: DL h;
2030: int i,n;
2031:
2032: if ( !p ) {
2033: *rp = p; return;
2034: }
2035: n = p->nv;
2036: for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) {
2037: h = m->dl;
2038: while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) )
2039: i++;
2040: muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t);
2041: addd(CO,s,t,&u); s = u;
1.7 noro 2042: }
2043: *rp = s;
1.5 noro 2044: }
2045:
1.7 noro 2046: /*
2047: * setting flags
1.30 noro 2048: * call create_order_spec with vl=0 to set old type order.
1.7 noro 2049: *
2050: */
2051:
1.27 noro 2052: int create_order_spec(VL vl,Obj obj,struct order_spec **specp)
1.5 noro 2053: {
1.62 noro 2054: int i,j,n,s,row,col,ret,wlen;
1.27 noro 2055: struct order_spec *spec;
1.7 noro 2056: struct order_pair *l;
1.62 noro 2057: Obj wp,wm;
2058: NODE node,t,tn,wpair;
1.7 noro 2059: MAT m;
1.49 noro 2060: VECT v;
2061: pointer **b,*bv;
1.7 noro 2062: int **w;
1.5 noro 2063:
1.37 noro 2064: if ( vl && obj && OID(obj) == O_LIST ) {
2065: ret = create_composite_order_spec(vl,(LIST)obj,specp);
2066: if ( show_orderspec )
2067: print_composite_order_spec(*specp);
2068: return ret;
2069: }
1.27 noro 2070:
2071: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 2072: if ( !obj || NUM(obj) ) {
2073: spec->id = 0; spec->obj = obj;
2074: spec->ord.simple = QTOS((Q)obj);
2075: return 1;
2076: } else if ( OID(obj) == O_LIST ) {
1.62 noro 2077: /* module order; obj = [0|1,w,ord] or [0|1,ord] */
1.7 noro 2078: node = BDY((LIST)obj);
1.53 noro 2079: if ( !BDY(node) || NUM(BDY(node)) ) {
1.62 noro 2080: switch ( length(node) ) {
2081: case 2:
2082: create_order_spec(0,(Obj)BDY(NEXT(node)),&spec);
2083: spec->id += 256; spec->obj = obj;
2084: spec->top_weight = 0;
2085: spec->module_rank = 0;
2086: spec->module_top_weight = 0;
2087: spec->ispot = (BDY(node)!=0);
2088: if ( spec->ispot ) {
2089: n = QTOS((Q)BDY(node));
2090: if ( n < 0 )
2091: spec->pot_nelim = -n;
2092: else
2093: spec->pot_nelim = 0;
2094: }
2095: break;
2096:
2097: case 3:
2098: create_order_spec(0,(Obj)BDY(NEXT(NEXT(node))),&spec);
2099: spec->id += 256; spec->obj = obj;
2100: spec->ispot = (BDY(node)!=0);
2101: node = NEXT(node);
2102: if ( !BDY(node) || OID(BDY(node)) != O_LIST )
2103: error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
2104: wpair = BDY((LIST)BDY(node));
2105: if ( length(wpair) != 2 )
2106: error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
2107:
2108: wp = BDY(wpair);
2109: wm = BDY(NEXT(wpair));
2110: if ( !wp || OID(wp) != O_LIST || !wm || OID(wm) != O_LIST )
2111: error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight");
2112: spec->nv = length(BDY((LIST)wp));
2113: spec->top_weight = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
2114: for ( i = 0, t = BDY((LIST)wp); i < spec->nv; t = NEXT(t), i++ )
2115: spec->top_weight[i] = QTOS((Q)BDY(t));
2116:
2117: spec->module_rank = length(BDY((LIST)wm));
2118: spec->module_top_weight = (int *)MALLOC_ATOMIC(spec->module_rank*sizeof(int));
2119: for ( i = 0, t = BDY((LIST)wm); i < spec->module_rank; t = NEXT(t), i++ )
2120: spec->module_top_weight[i] = QTOS((Q)BDY(t));
2121: break;
2122: default:
2123: error("create_order_spec : invalid arguments for module order");
2124: }
2125:
1.53 noro 2126: *specp = spec;
2127: return 1;
1.62 noro 2128: } else {
2129: /* block order in polynomial ring */
2130: for ( n = 0, t = node; t; t = NEXT(t), n++ );
2131: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
2132: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
2133: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
2134: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
2135: s += l[i].length;
2136: }
2137: spec->id = 1; spec->obj = obj;
2138: spec->ord.block.order_pair = l;
2139: spec->ord.block.length = n; spec->nv = s;
2140: return 1;
2141: }
1.7 noro 2142: } else if ( OID(obj) == O_MAT ) {
2143: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
2144: w = almat(row,col);
2145: for ( i = 0; i < row; i++ )
2146: for ( j = 0; j < col; j++ )
2147: w[i][j] = QTOS((Q)b[i][j]);
2148: spec->id = 2; spec->obj = obj;
2149: spec->nv = col; spec->ord.matrix.row = row;
2150: spec->ord.matrix.matrix = w;
2151: return 1;
2152: } else
1.5 noro 2153: return 0;
2154: }
2155:
1.28 noro 2156: void print_composite_order_spec(struct order_spec *spec)
2157: {
2158: int nv,n,len,i,j,k,start;
2159: struct weight_or_block *worb;
2160:
2161: nv = spec->nv;
2162: n = spec->ord.composite.length;
2163: worb = spec->ord.composite.w_or_b;
2164: for ( i = 0; i < n; i++, worb++ ) {
2165: len = worb->length;
2166: printf("[ ");
2167: switch ( worb->type ) {
2168: case IS_DENSE_WEIGHT:
2169: for ( j = 0; j < len; j++ )
2170: printf("%d ",worb->body.dense_weight[j]);
2171: for ( ; j < nv; j++ )
2172: printf("0 ");
2173: break;
2174: case IS_SPARSE_WEIGHT:
2175: for ( j = 0, k = 0; j < nv; j++ )
2176: if ( j == worb->body.sparse_weight[k].pos )
2177: printf("%d ",worb->body.sparse_weight[k++].value);
2178: else
2179: printf("0 ");
2180: break;
2181: case IS_BLOCK:
2182: start = worb->body.block.start;
2183: for ( j = 0; j < start; j++ ) printf("0 ");
2184: switch ( worb->body.block.order ) {
2185: case 0:
2186: for ( k = 0; k < len; k++, j++ ) printf("R ");
2187: break;
2188: case 1:
2189: for ( k = 0; k < len; k++, j++ ) printf("G ");
2190: break;
2191: case 2:
2192: for ( k = 0; k < len; k++, j++ ) printf("L ");
2193: break;
2194: }
2195: for ( ; j < nv; j++ ) printf("0 ");
2196: break;
2197: }
2198: printf("]\n");
2199: }
1.38 noro 2200: }
2201:
2202: struct order_spec *append_block(struct order_spec *spec,
2203: int nv,int nalg,int ord)
2204: {
2205: MAT m,mat;
2206: int i,j,row,col,n;
2207: Q **b,**wp;
2208: int **w;
2209: NODE t,s,s0;
2210: struct order_pair *l,*l0;
2211: int n0,nv0;
2212: LIST list0,list1,list;
2213: Q oq,nq;
2214: struct order_spec *r;
2215:
2216: r = (struct order_spec *)MALLOC(sizeof(struct order_spec));
2217: switch ( spec->id ) {
2218: case 0:
2219: STOQ(spec->ord.simple,oq); STOQ(nv,nq);
2220: t = mknode(2,oq,nq); MKLIST(list0,t);
2221: STOQ(ord,oq); STOQ(nalg,nq);
2222: t = mknode(2,oq,nq); MKLIST(list1,t);
2223: t = mknode(2,list0,list1); MKLIST(list,t);
2224: l = (struct order_pair *)MALLOC_ATOMIC(2*sizeof(struct order_pair));
2225: l[0].order = spec->ord.simple; l[0].length = nv;
2226: l[1].order = ord; l[1].length = nalg;
2227: r->id = 1; r->obj = (Obj)list;
2228: r->ord.block.order_pair = l;
2229: r->ord.block.length = 2;
2230: r->nv = nv+nalg;
2231: break;
2232: case 1:
2233: if ( spec->nv != nv )
2234: error("append_block : number of variables mismatch");
2235: l0 = spec->ord.block.order_pair;
2236: n0 = spec->ord.block.length;
2237: nv0 = spec->nv;
2238: list0 = (LIST)spec->obj;
2239: n = n0+1;
2240: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
2241: for ( i = 0; i < n0; i++ )
2242: l[i] = l0[i];
2243: l[i].order = ord; l[i].length = nalg;
2244: for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) {
2245: NEXTNODE(s0,s); BDY(s) = BDY(t);
2246: }
2247: STOQ(ord,oq); STOQ(nalg,nq);
2248: t = mknode(2,oq,nq); MKLIST(list,t);
2249: NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0;
2250: MKLIST(list,s0);
2251: r->id = 1; r->obj = (Obj)list;
2252: r->ord.block.order_pair = l;
2253: r->ord.block.length = n;
2254: r->nv = nv+nalg;
2255: break;
2256: case 2:
2257: if ( spec->nv != nv )
2258: error("append_block : number of variables mismatch");
2259: m = (MAT)spec->obj;
2260: row = m->row; col = m->col; b = (Q **)BDY(m);
2261: w = almat(row+nalg,col+nalg);
2262: MKMAT(mat,row+nalg,col+nalg); wp = (Q **)BDY(mat);
2263: for ( i = 0; i < row; i++ )
2264: for ( j = 0; j < col; j++ ) {
2265: w[i][j] = QTOS(b[i][j]);
2266: wp[i][j] = b[i][j];
2267: }
2268: for ( i = 0; i < nalg; i++ ) {
2269: w[i+row][i+col] = 1;
2270: wp[i+row][i+col] = ONE;
2271: }
2272: r->id = 2; r->obj = (Obj)mat;
2273: r->nv = col+nalg; r->ord.matrix.row = row+nalg;
2274: r->ord.matrix.matrix = w;
2275: break;
2276: case 3:
2277: default:
2278: /* XXX */
2279: error("append_block : not implemented yet");
2280: }
2281: return r;
1.28 noro 2282: }
2283:
1.37 noro 2284: int comp_sw(struct sparse_weight *a, struct sparse_weight *b)
2285: {
2286: if ( a->pos > b->pos ) return 1;
2287: else if ( a->pos < b->pos ) return -1;
2288: else return 0;
2289: }
2290:
1.27 noro 2291: /* order = [w_or_b, w_or_b, ... ] */
2292: /* w_or_b = w or b */
2293: /* w = [1,2,...] or [x,1,y,2,...] */
2294: /* b = [@lex,x,y,...,z] etc */
2295:
2296: int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp)
2297: {
2298: NODE wb,t,p;
2299: struct order_spec *spec;
2300: VL tvl;
1.29 noro 2301: int n,i,j,k,l,start,end,len,w;
1.27 noro 2302: int *dw;
2303: struct sparse_weight *sw;
2304: struct weight_or_block *w_or_b;
2305: Obj a0;
2306: NODE a;
1.29 noro 2307: V v,sv,ev;
2308: SYMBOL sym;
2309: int *top;
1.27 noro 2310:
2311: /* l = number of vars in vl */
2312: for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ );
2313: /* n = number of primitives in order */
2314: wb = BDY(order);
2315: n = length(wb);
2316: *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec));
2317: spec->id = 3;
2318: spec->obj = (Obj)order;
2319: spec->nv = l;
2320: spec->ord.composite.length = n;
1.28 noro 2321: w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *)
1.29 noro 2322: MALLOC(sizeof(struct weight_or_block)*(n+1));
2323:
2324: /* top : register the top variable in each w_or_b specification */
2325: top = (int *)ALLOCA(l*sizeof(int));
2326: for ( i = 0; i < l; i++ ) top[i] = 0;
2327:
1.28 noro 2328: for ( t = wb, i = 0; t; t = NEXT(t), i++ ) {
1.30 noro 2329: if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST )
2330: error("a list of lists must be specified for the key \"order\"");
1.28 noro 2331: a = BDY((LIST)BDY(t));
1.27 noro 2332: len = length(a);
2333: a0 = (Obj)BDY(a);
2334: if ( !a0 || OID(a0) == O_N ) {
1.28 noro 2335: /* a is a dense weight vector */
1.27 noro 2336: dw = (int *)MALLOC(sizeof(int)*len);
1.30 noro 2337: for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) {
2338: if ( !INT((Q)BDY(p)) )
2339: error("a dense weight vector must be specified as a list of integers");
1.27 noro 2340: dw[j] = QTOS((Q)BDY(p));
1.30 noro 2341: }
1.27 noro 2342: w_or_b[i].type = IS_DENSE_WEIGHT;
2343: w_or_b[i].length = len;
2344: w_or_b[i].body.dense_weight = dw;
1.29 noro 2345:
2346: /* find the top */
2347: for ( k = 0; k < len && !dw[k]; k++ );
2348: if ( k < len ) top[k] = 1;
2349:
1.27 noro 2350: } else if ( OID(a0) == O_P ) {
1.28 noro 2351: /* a is a sparse weight vector */
2352: len >>= 1;
1.27 noro 2353: sw = (struct sparse_weight *)
2354: MALLOC(sizeof(struct sparse_weight)*len);
2355: for ( j = 0, p = a; j < len; j++ ) {
1.30 noro 2356: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
2357: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 2358: v = VR((P)BDY(p)); p = NEXT(p);
1.27 noro 2359: for ( tvl = vl, k = 0; tvl && tvl->v != v;
2360: k++, tvl = NEXT(tvl) );
2361: if ( !tvl )
1.30 noro 2362: error("invalid variable name in a sparse weight vector");
1.27 noro 2363: sw[j].pos = k;
1.30 noro 2364: if ( !INT((Q)BDY(p)) )
2365: error("a sparse weight vector must be specified as [var1,weight1,...]");
1.28 noro 2366: sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p);
1.27 noro 2367: }
1.37 noro 2368: qsort(sw,len,sizeof(struct sparse_weight),
2369: (int (*)(const void *,const void *))comp_sw);
1.27 noro 2370: w_or_b[i].type = IS_SPARSE_WEIGHT;
2371: w_or_b[i].length = len;
2372: w_or_b[i].body.sparse_weight = sw;
1.29 noro 2373:
2374: /* find the top */
2375: for ( k = 0; k < len && !sw[k].value; k++ );
2376: if ( k < len ) top[sw[k].pos] = 1;
2377: } else if ( OID(a0) == O_RANGE ) {
2378: /* [range(v1,v2),w] */
2379: sv = VR((P)(((RANGE)a0)->start));
2380: ev = VR((P)(((RANGE)a0)->end));
2381: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
2382: if ( !tvl )
2383: error("invalid range");
2384: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
2385: if ( !tvl )
2386: error("invalid range");
2387: len = end-start+1;
2388: sw = (struct sparse_weight *)
2389: MALLOC(sizeof(struct sparse_weight)*len);
2390: w = QTOS((Q)BDY(NEXT(a)));
2391: for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) );
2392: for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) {
2393: sw[j].pos = k;
2394: sw[j].value = w;
2395: }
2396: w_or_b[i].type = IS_SPARSE_WEIGHT;
2397: w_or_b[i].length = len;
2398: w_or_b[i].body.sparse_weight = sw;
2399:
2400: /* register the top */
2401: if ( w ) top[start] = 1;
1.28 noro 2402: } else if ( OID(a0) == O_SYMBOL ) {
2403: /* a is a block */
1.29 noro 2404: sym = (SYMBOL)a0; a = NEXT(a); len--;
2405: if ( OID((Obj)BDY(a)) == O_RANGE ) {
2406: sv = VR((P)(((RANGE)BDY(a))->start));
2407: ev = VR((P)(((RANGE)BDY(a))->end));
2408: for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) );
2409: if ( !tvl )
2410: error("invalid range");
2411: for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) );
2412: if ( !tvl )
2413: error("invalid range");
2414: len = end-start+1;
2415: } else {
2416: for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a));
1.28 noro 2417: tvl = NEXT(tvl), start++ );
1.29 noro 2418: for ( p = NEXT(a), tvl = NEXT(tvl); p;
1.30 noro 2419: p = NEXT(p), tvl = NEXT(tvl) ) {
2420: if ( !BDY(p) || OID((P)BDY(p)) != O_P )
2421: error("a block must be specified as [ordsymbol,var1,var2,...]");
1.29 noro 2422: if ( tvl->v != VR((P)BDY(p)) ) break;
1.30 noro 2423: }
1.29 noro 2424: if ( p )
1.30 noro 2425: error("a block must be contiguous in the variable list");
1.29 noro 2426: }
1.28 noro 2427: w_or_b[i].type = IS_BLOCK;
2428: w_or_b[i].length = len;
2429: w_or_b[i].body.block.start = start;
2430: if ( !strcmp(sym->name,"@grlex") )
2431: w_or_b[i].body.block.order = 0;
2432: else if ( !strcmp(sym->name,"@glex") )
2433: w_or_b[i].body.block.order = 1;
2434: else if ( !strcmp(sym->name,"@lex") )
2435: w_or_b[i].body.block.order = 2;
2436: else
1.29 noro 2437: error("invalid ordername");
2438: /* register the tops */
2439: for ( j = 0, k = start; j < len; j++, k++ )
2440: top[k] = 1;
1.28 noro 2441: }
1.29 noro 2442: }
2443: for ( k = 0; k < l && top[k]; k++ );
2444: if ( k < l ) {
2445: /* incomplete order specification; add @grlex */
2446: w_or_b[n].type = IS_BLOCK;
2447: w_or_b[n].length = l;
2448: w_or_b[n].body.block.start = 0;
2449: w_or_b[n].body.block.order = 0;
2450: spec->ord.composite.length = n+1;
1.27 noro 2451: }
2452: }
2453:
1.35 noro 2454: /* module order spec */
2455:
2456: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s)
2457: {
2458: struct modorder_spec *spec;
2459: NODE n,t;
2460: LIST list;
2461: int *ds;
2462: int i,l;
2463: Q q;
2464:
2465: *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec));
2466: spec->id = id;
2467: if ( shift ) {
2468: n = BDY(shift);
2469: spec->len = l = length(n);
2470: spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int));
2471: for ( t = n, i = 0; t; t = NEXT(t), i++ )
2472: ds[i] = QTOS((Q)BDY(t));
2473: } else {
2474: spec->len = 0;
2475: spec->degree_shift = 0;
2476: }
2477: STOQ(id,q);
2478: n = mknode(2,q,shift);
2479: MKLIST(list,n);
2480: spec->obj = (Obj)list;
2481: }
2482:
1.7 noro 2483: /*
2484: * converters
2485: *
2486: */
2487:
1.20 noro 2488: void dp_homo(DP p,DP *rp)
1.5 noro 2489: {
1.7 noro 2490: MP m,mr,mr0;
2491: int i,n,nv,td;
2492: DL dl,dlh;
1.5 noro 2493:
1.7 noro 2494: if ( !p )
2495: *rp = 0;
2496: else {
2497: n = p->nv; nv = n + 1;
2498: m = BDY(p); td = sugard(m);
2499: for ( mr0 = 0; m; m = NEXT(m) ) {
2500: NEXTMP(mr0,mr); mr->c = m->c;
2501: dl = m->dl;
2502: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
2503: dlh->td = td;
2504: for ( i = 0; i < n; i++ )
2505: dlh->d[i] = dl->d[i];
2506: dlh->d[n] = td - dl->td;
2507: }
2508: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1.5 noro 2509: }
2510: }
2511:
1.20 noro 2512: void dp_dehomo(DP p,DP *rp)
1.5 noro 2513: {
1.7 noro 2514: MP m,mr,mr0;
2515: int i,n,nv;
2516: DL dl,dlh;
1.5 noro 2517:
1.7 noro 2518: if ( !p )
2519: *rp = 0;
2520: else {
2521: n = p->nv; nv = n - 1;
2522: m = BDY(p);
2523: for ( mr0 = 0; m; m = NEXT(m) ) {
2524: NEXTMP(mr0,mr); mr->c = m->c;
2525: dlh = m->dl;
2526: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
2527: dl->td = dlh->td - dlh->d[nv];
2528: for ( i = 0; i < nv; i++ )
2529: dl->d[i] = dlh->d[i];
2530: }
2531: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
2532: }
1.5 noro 2533: }
2534:
1.20 noro 2535: void dp_mod(DP p,int mod,NODE subst,DP *rp)
1.5 noro 2536: {
1.7 noro 2537: MP m,mr,mr0;
2538: P t,s,s1;
2539: V v;
2540: NODE tn;
1.5 noro 2541:
1.7 noro 2542: if ( !p )
2543: *rp = 0;
2544: else {
2545: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1.66 ! noro 2546: for ( tn = subst, s = (P)m->c; tn; tn = NEXT(tn) ) {
1.7 noro 2547: v = VR((P)BDY(tn)); tn = NEXT(tn);
2548: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
2549: }
2550: ptomp(mod,s,&t);
2551: if ( t ) {
1.66 ! noro 2552: NEXTMP(mr0,mr); mr->c = (Obj)t; mr->dl = m->dl;
1.7 noro 2553: }
2554: }
2555: if ( mr0 ) {
2556: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
2557: } else
2558: *rp = 0;
2559: }
1.5 noro 2560: }
2561:
1.20 noro 2562: void dp_rat(DP p,DP *rp)
1.5 noro 2563: {
1.7 noro 2564: MP m,mr,mr0;
1.5 noro 2565:
1.7 noro 2566: if ( !p )
2567: *rp = 0;
2568: else {
2569: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1.66 ! noro 2570: NEXTMP(mr0,mr); mptop((P)m->c,(P *)&mr->c); mr->dl = m->dl;
1.5 noro 2571: }
1.7 noro 2572: if ( mr0 ) {
2573: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
2574: } else
2575: *rp = 0;
1.5 noro 2576: }
2577: }
2578:
2579:
1.27 noro 2580: void homogenize_order(struct order_spec *old,int n,struct order_spec **newp)
1.5 noro 2581: {
1.7 noro 2582: struct order_pair *l;
2583: int length,nv,row,i,j;
2584: int **newm,**oldm;
1.27 noro 2585: struct order_spec *new;
1.31 noro 2586: int onv,nnv,nlen,olen,owlen;
2587: struct weight_or_block *owb,*nwb;
1.5 noro 2588:
1.27 noro 2589: *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec));
1.7 noro 2590: switch ( old->id ) {
2591: case 0:
2592: switch ( old->ord.simple ) {
2593: case 0:
2594: new->id = 0; new->ord.simple = 0; break;
2595: case 1:
2596: l = (struct order_pair *)
2597: MALLOC_ATOMIC(2*sizeof(struct order_pair));
2598: l[0].length = n; l[0].order = 1;
2599: l[1].length = 1; l[1].order = 2;
2600: new->id = 1;
2601: new->ord.block.order_pair = l;
2602: new->ord.block.length = 2; new->nv = n+1;
2603: break;
2604: case 2:
2605: new->id = 0; new->ord.simple = 1; break;
2606: case 3: case 4: case 5:
2607: new->id = 0; new->ord.simple = old->ord.simple+3;
2608: dp_nelim = n-1; break;
2609: case 6: case 7: case 8: case 9:
2610: new->id = 0; new->ord.simple = old->ord.simple; break;
2611: default:
2612: error("homogenize_order : invalid input");
2613: }
2614: break;
1.50 noro 2615: case 1: case 257:
1.7 noro 2616: length = old->ord.block.length;
2617: l = (struct order_pair *)
2618: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
2619: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
2620: l[length].order = 2; l[length].length = 1;
1.50 noro 2621: new->id = old->id; new->nv = n+1;
1.7 noro 2622: new->ord.block.order_pair = l;
2623: new->ord.block.length = length+1;
1.51 noro 2624: new->ispot = old->ispot;
1.7 noro 2625: break;
1.50 noro 2626: case 2: case 258:
1.7 noro 2627: nv = old->nv; row = old->ord.matrix.row;
2628: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
2629: for ( i = 0; i <= nv; i++ )
2630: newm[0][i] = 1;
2631: for ( i = 0; i < row; i++ ) {
2632: for ( j = 0; j < nv; j++ )
2633: newm[i+1][j] = oldm[i][j];
2634: newm[i+1][j] = 0;
2635: }
1.50 noro 2636: new->id = old->id; new->nv = nv+1;
1.7 noro 2637: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
1.51 noro 2638: new->ispot = old->ispot;
1.31 noro 2639: break;
1.50 noro 2640: case 3: case 259:
1.31 noro 2641: onv = old->nv;
2642: nnv = onv+1;
2643: olen = old->ord.composite.length;
2644: nlen = olen+1;
2645: owb = old->ord.composite.w_or_b;
2646: nwb = (struct weight_or_block *)
2647: MALLOC(nlen*sizeof(struct weight_or_block));
2648: for ( i = 0; i < olen; i++ ) {
2649: nwb[i].type = owb[i].type;
2650: switch ( owb[i].type ) {
2651: case IS_DENSE_WEIGHT:
2652: owlen = owb[i].length;
2653: nwb[i].length = owlen+1;
2654: nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int));
2655: for ( j = 0; j < owlen; j++ )
2656: nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j];
2657: nwb[i].body.dense_weight[owlen] = 0;
2658: break;
2659: case IS_SPARSE_WEIGHT:
2660: nwb[i].length = owb[i].length;
2661: nwb[i].body.sparse_weight = owb[i].body.sparse_weight;
2662: break;
2663: case IS_BLOCK:
2664: nwb[i].length = owb[i].length;
2665: nwb[i].body.block = owb[i].body.block;
2666: break;
2667: }
2668: }
2669: nwb[i].type = IS_SPARSE_WEIGHT;
2670: nwb[i].body.sparse_weight =
2671: (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
2672: nwb[i].body.sparse_weight[0].pos = onv;
2673: nwb[i].body.sparse_weight[0].value = 1;
1.50 noro 2674: new->id = old->id;
1.31 noro 2675: new->nv = nnv;
2676: new->ord.composite.length = nlen;
2677: new->ord.composite.w_or_b = nwb;
1.51 noro 2678: new->ispot = old->ispot;
1.31 noro 2679: print_composite_order_spec(new);
1.7 noro 2680: break;
1.50 noro 2681: case 256: /* simple module order */
2682: switch ( old->ord.simple ) {
2683: case 0:
2684: new->id = 256; new->ord.simple = 0; break;
2685: case 1:
2686: l = (struct order_pair *)
2687: MALLOC_ATOMIC(2*sizeof(struct order_pair));
2688: l[0].length = n; l[0].order = old->ord.simple;
2689: l[1].length = 1; l[1].order = 2;
2690: new->id = 257;
2691: new->ord.block.order_pair = l;
2692: new->ord.block.length = 2; new->nv = n+1;
2693: break;
2694: case 2:
2695: new->id = 256; new->ord.simple = 1; break;
2696: default:
2697: error("homogenize_order : invalid input");
2698: }
1.51 noro 2699: new->ispot = old->ispot;
1.50 noro 2700: break;
1.7 noro 2701: default:
2702: error("homogenize_order : invalid input");
1.5 noro 2703: }
1.7 noro 2704: }
2705:
1.20 noro 2706: void qltozl(Q *w,int n,Q *dvr)
1.7 noro 2707: {
2708: N nm,dn;
2709: N g,l1,l2,l3;
2710: Q c,d;
2711: int i;
2712: struct oVECT v;
1.5 noro 2713:
2714: for ( i = 0; i < n; i++ )
1.7 noro 2715: if ( w[i] && !INT(w[i]) )
2716: break;
2717: if ( i == n ) {
2718: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
2719: igcdv(&v,dvr); return;
2720: }
1.56 noro 2721: for ( i = 0; !w[i]; i++ );
2722: c = w[i]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
2723: for ( i++; i < n; i++ ) {
2724: c = w[i];
2725: if ( !c ) continue;
2726: l1 = INT(c) ? ONEN : DN(c);
1.7 noro 2727: gcdn(nm,NM(c),&g); nm = g;
2728: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1.5 noro 2729: }
1.7 noro 2730: if ( UNIN(dn) )
2731: NTOQ(nm,1,d);
2732: else
2733: NDTOQ(nm,dn,1,d);
2734: *dvr = d;
2735: }
1.5 noro 2736:
1.20 noro 2737: int comp_nm(Q *a,Q *b)
1.7 noro 2738: {
2739: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
2740: }
2741:
1.20 noro 2742: void sortbynm(Q *w,int n)
1.7 noro 2743: {
2744: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
2745: }
1.5 noro 2746:
2747:
1.7 noro 2748: /*
2749: * simple operations
2750: *
2751: */
1.5 noro 2752:
1.20 noro 2753: int dp_redble(DP p1,DP p2)
1.7 noro 2754: {
2755: int i,n;
2756: DL d1,d2;
1.5 noro 2757:
1.7 noro 2758: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
2759: if ( d1->td < d2->td )
2760: return 0;
2761: else {
2762: for ( i = 0, n = p1->nv; i < n; i++ )
2763: if ( d1->d[i] < d2->d[i] )
2764: return 0;
2765: return 1;
1.5 noro 2766: }
2767: }
2768:
1.66 ! noro 2769: int dpm_redble(DPM p1,DPM p2)
! 2770: {
! 2771: int i,n;
! 2772: DL d1,d2;
! 2773:
! 2774: if ( BDY(p1)->pos != BDY(p2)->pos ) return 0;
! 2775: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
! 2776: if ( d1->td < d2->td )
! 2777: return 0;
! 2778: else {
! 2779: for ( i = 0, n = p1->nv; i < n; i++ )
! 2780: if ( d1->d[i] < d2->d[i] )
! 2781: return 0;
! 2782: return 1;
! 2783: }
! 2784: }
! 2785:
! 2786:
1.20 noro 2787: void dp_subd(DP p1,DP p2,DP *rp)
1.5 noro 2788: {
1.7 noro 2789: int i,n;
1.5 noro 2790: DL d1,d2,d;
2791: MP m;
1.7 noro 2792: DP s;
1.5 noro 2793:
2794: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1.7 noro 2795: NEWDL(d,n); d->td = d1->td - d2->td;
1.5 noro 2796: for ( i = 0; i < n; i++ )
1.7 noro 2797: d->d[i] = d1->d[i]-d2->d[i];
1.66 ! noro 2798: NEWMP(m); m->dl = d; m->c = (Obj)ONE; NEXT(m) = 0;
1.7 noro 2799: MKDP(n,m,s); s->sugar = d->td;
2800: *rp = s;
2801: }
2802:
1.20 noro 2803: void dltod(DL d,int n,DP *rp)
1.7 noro 2804: {
2805: MP m;
2806: DP s;
2807:
1.66 ! noro 2808: NEWMP(m); m->dl = d; m->c = (Obj)ONE; NEXT(m) = 0;
1.7 noro 2809: MKDP(n,m,s); s->sugar = d->td;
2810: *rp = s;
1.5 noro 2811: }
2812:
1.20 noro 2813: void dp_hm(DP p,DP *rp)
1.5 noro 2814: {
2815: MP m,mr;
2816:
2817: if ( !p )
2818: *rp = 0;
2819: else {
2820: m = BDY(p);
2821: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
2822: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
2823: }
2824: }
2825:
1.35 noro 2826: void dp_ht(DP p,DP *rp)
2827: {
2828: MP m,mr;
2829:
2830: if ( !p )
2831: *rp = 0;
2832: else {
2833: m = BDY(p);
1.66 ! noro 2834: NEWMP(mr); mr->dl = m->dl; mr->c = (Obj)ONE; NEXT(mr) = 0;
1.35 noro 2835: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
2836: }
2837: }
2838:
1.66 ! noro 2839: void dpm_hm(DPM p,DPM *rp)
! 2840: {
! 2841: DMM m,mr;
! 2842:
! 2843: if ( !p )
! 2844: *rp = 0;
! 2845: else {
! 2846: m = BDY(p);
! 2847: NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; NEXT(mr) = 0;
! 2848: MKDPM(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
! 2849: }
! 2850: }
! 2851:
! 2852: void dpm_ht(DPM p,DPM *rp)
! 2853: {
! 2854: DMM m,mr;
! 2855:
! 2856: if ( !p )
! 2857: *rp = 0;
! 2858: else {
! 2859: m = BDY(p);
! 2860: NEWDMM(mr); mr->dl = m->dl; mr->pos = m->pos; mr->c = (Obj)ONE; NEXT(mr) = 0;
! 2861: MKDPM(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
! 2862: }
! 2863: }
! 2864:
! 2865:
1.20 noro 2866: void dp_rest(DP p,DP *rp)
1.5 noro 2867: {
2868: MP m;
2869:
2870: m = BDY(p);
2871: if ( !NEXT(m) )
2872: *rp = 0;
2873: else {
2874: MKDP(p->nv,NEXT(m),*rp);
2875: if ( *rp )
2876: (*rp)->sugar = p->sugar;
2877: }
2878: }
2879:
1.66 ! noro 2880: void dpm_rest(DPM p,DPM *rp)
! 2881: {
! 2882: DMM m;
! 2883:
! 2884: m = BDY(p);
! 2885: if ( !NEXT(m) )
! 2886: *rp = 0;
! 2887: else {
! 2888: MKDPM(p->nv,NEXT(m),*rp);
! 2889: if ( *rp )
! 2890: (*rp)->sugar = p->sugar;
! 2891: }
! 2892: }
! 2893:
1.20 noro 2894: DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl)
1.5 noro 2895: {
1.21 noro 2896: register int i, *d1, *d2, *d, td;
1.5 noro 2897:
2898: if ( !dl ) NEWDL(dl,nv);
2899: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1.21 noro 2900: for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) {
2901: *d = *d1 > *d2 ? *d1 : *d2;
2902: td += MUL_WEIGHT(*d,i);
2903: }
1.5 noro 2904: dl->td = td;
2905: return dl;
2906: }
2907:
1.20 noro 2908: int dl_equal(int nv,DL dl1,DL dl2)
1.5 noro 2909: {
2910: register int *d1, *d2, n;
2911:
2912: if ( dl1->td != dl2->td ) return 0;
2913: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
2914: if ( *d1 != *d2 ) return 0;
2915: return 1;
2916: }
2917:
1.20 noro 2918: int dp_nt(DP p)
1.5 noro 2919: {
2920: int i;
2921: MP m;
2922:
2923: if ( !p )
2924: return 0;
2925: else {
2926: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
2927: return i;
2928: }
2929: }
2930:
1.20 noro 2931: int dp_homogeneous(DP p)
1.15 noro 2932: {
2933: MP m;
2934: int d;
2935:
2936: if ( !p )
2937: return 1;
2938: else {
2939: m = BDY(p);
2940: d = m->dl->td;
2941: m = NEXT(m);
2942: for ( ; m; m = NEXT(m) ) {
2943: if ( m->dl->td != d )
2944: return 0;
2945: }
2946: return 1;
2947: }
1.16 noro 2948: }
2949:
1.20 noro 2950: void _print_mp(int nv,MP m)
1.16 noro 2951: {
2952: int i;
2953:
1.17 noro 2954: if ( !m )
1.16 noro 2955: return;
2956: for ( ; m; m = NEXT(m) ) {
2957: fprintf(stderr,"%d<",ITOS(C(m)));
2958: for ( i = 0; i < nv; i++ ) {
2959: fprintf(stderr,"%d",m->dl->d[i]);
2960: if ( i != nv-1 )
2961: fprintf(stderr," ");
2962: }
2963: fprintf(stderr,">",C(m));
2964: }
2965: fprintf(stderr,"\n");
1.15 noro 2966: }
1.26 noro 2967:
2968: static int cmp_mp_nvar;
2969:
2970: int comp_mp(MP *a,MP *b)
2971: {
2972: return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl);
2973: }
2974:
2975: void dp_sort(DP p,DP *rp)
2976: {
2977: MP t,mp,mp0;
2978: int i,n;
2979: DP r;
2980: MP *w;
2981:
2982: if ( !p ) {
2983: *rp = 0;
2984: return;
2985: }
2986: for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ );
2987: w = (MP *)ALLOCA(n*sizeof(MP));
2988: for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ )
2989: w[i] = t;
2990: cmp_mp_nvar = NV(p);
2991: qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp);
2992: mp0 = 0;
2993: for ( i = n-1; i >= 0; i-- ) {
2994: NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]);
2995: NEXT(mp) = mp0; mp0 = mp;
2996: }
2997: MKDP(p->nv,mp0,r);
2998: r->sugar = p->sugar;
2999: *rp = r;
3000: }
3001:
1.32 noro 3002: DP extract_initial_term_from_dp(DP p,int *weight,int n);
3003: LIST extract_initial_term(LIST f,int *weight,int n);
3004:
3005: DP extract_initial_term_from_dp(DP p,int *weight,int n)
3006: {
1.34 noro 3007: int w,t,i,top;
1.32 noro 3008: MP m,r0,r;
3009: DP dp;
3010:
3011: if ( !p ) return 0;
1.34 noro 3012: top = 1;
1.32 noro 3013: for ( m = BDY(p); m; m = NEXT(m) ) {
3014: for ( i = 0, t = 0; i < n; i++ )
3015: t += weight[i]*m->dl->d[i];
1.34 noro 3016: if ( top || t > w ) {
1.32 noro 3017: r0 = 0;
3018: w = t;
1.34 noro 3019: top = 0;
1.32 noro 3020: }
3021: if ( t == w ) {
3022: NEXTMP(r0,r);
3023: r->dl = m->dl;
3024: r->c = m->c;
3025: }
3026: }
3027: NEXT(r) = 0;
3028: MKDP(p->nv,r0,dp);
3029: return dp;
3030: }
3031:
3032: LIST extract_initial_term(LIST f,int *weight,int n)
3033: {
3034: NODE nd,r0,r;
3035: Obj p;
3036: LIST l;
3037:
3038: nd = BDY(f);
3039: for ( r0 = 0; nd; nd = NEXT(nd) ) {
3040: NEXTNODE(r0,r);
3041: p = (Obj)BDY(nd);
3042: BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n);
3043: }
3044: if ( r0 ) NEXT(r) = 0;
3045: MKLIST(l,r0);
3046: return l;
3047: }
3048:
3049: LIST dp_initial_term(LIST f,struct order_spec *ord)
3050: {
3051: int n,l,i;
3052: struct weight_or_block *worb;
3053: int *weight;
3054:
3055: switch ( ord->id ) {
3056: case 2: /* matrix order */
3057: /* extract the first row */
3058: n = ord->nv;
3059: weight = ord->ord.matrix.matrix[0];
3060: return extract_initial_term(f,weight,n);
3061: case 3: /* composite order */
3062: /* the first w_or_b */
3063: worb = ord->ord.composite.w_or_b;
3064: switch ( worb->type ) {
3065: case IS_DENSE_WEIGHT:
3066: n = worb->length;
3067: weight = worb->body.dense_weight;
3068: return extract_initial_term(f,weight,n);
3069: case IS_SPARSE_WEIGHT:
3070: n = ord->nv;
3071: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 3072: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 3073: l = worb->length;
3074: for ( i = 0; i < l; i++ )
3075: weight[worb->body.sparse_weight[i].pos]
3076: = worb->body.sparse_weight[i].value;
3077: return extract_initial_term(f,weight,n);
3078: default:
3079: error("dp_initial_term : unsupported order");
3080: }
3081: default:
3082: error("dp_initial_term : unsupported order");
3083: }
3084: }
3085:
3086: int highest_order_dp(DP p,int *weight,int n);
3087: LIST highest_order(LIST f,int *weight,int n);
3088:
3089: int highest_order_dp(DP p,int *weight,int n)
3090: {
1.34 noro 3091: int w,t,i,top;
1.32 noro 3092: MP m;
3093:
3094: if ( !p ) return -1;
1.34 noro 3095: top = 1;
1.32 noro 3096: for ( m = BDY(p); m; m = NEXT(m) ) {
3097: for ( i = 0, t = 0; i < n; i++ )
3098: t += weight[i]*m->dl->d[i];
1.34 noro 3099: if ( top || t > w ) {
1.32 noro 3100: w = t;
1.34 noro 3101: top = 0;
3102: }
1.32 noro 3103: }
3104: return w;
3105: }
3106:
3107: LIST highest_order(LIST f,int *weight,int n)
3108: {
3109: int h;
3110: NODE nd,r0,r;
3111: Obj p;
3112: LIST l;
3113: Q q;
3114:
3115: nd = BDY(f);
3116: for ( r0 = 0; nd; nd = NEXT(nd) ) {
3117: NEXTNODE(r0,r);
3118: p = (Obj)BDY(nd);
3119: h = highest_order_dp((DP)p,weight,n);
3120: STOQ(h,q);
3121: BDY(r) = (pointer)q;
3122: }
3123: if ( r0 ) NEXT(r) = 0;
3124: MKLIST(l,r0);
3125: return l;
3126: }
3127:
3128: LIST dp_order(LIST f,struct order_spec *ord)
3129: {
3130: int n,l,i;
3131: struct weight_or_block *worb;
3132: int *weight;
3133:
3134: switch ( ord->id ) {
3135: case 2: /* matrix order */
3136: /* extract the first row */
3137: n = ord->nv;
3138: weight = ord->ord.matrix.matrix[0];
3139: return highest_order(f,weight,n);
3140: case 3: /* composite order */
3141: /* the first w_or_b */
3142: worb = ord->ord.composite.w_or_b;
3143: switch ( worb->type ) {
3144: case IS_DENSE_WEIGHT:
3145: n = worb->length;
3146: weight = worb->body.dense_weight;
3147: return highest_order(f,weight,n);
3148: case IS_SPARSE_WEIGHT:
3149: n = ord->nv;
3150: weight = (int *)ALLOCA(n*sizeof(int));
1.33 noro 3151: for ( i = 0; i < n; i++ ) weight[i] = 0;
1.32 noro 3152: l = worb->length;
3153: for ( i = 0; i < l; i++ )
3154: weight[worb->body.sparse_weight[i].pos]
3155: = worb->body.sparse_weight[i].value;
3156: return highest_order(f,weight,n);
3157: default:
3158: error("dp_initial_term : unsupported order");
3159: }
3160: default:
3161: error("dp_initial_term : unsupported order");
1.35 noro 3162: }
3163: }
3164:
3165: int dpv_ht(DPV p,DP *h)
3166: {
3167: int len,max,maxi,i,t;
3168: DP *e;
3169: MP m,mr;
3170:
3171: len = p->len;
3172: e = p->body;
3173: max = -1;
3174: maxi = -1;
3175: for ( i = 0; i < len; i++ )
3176: if ( e[i] && (t = BDY(e[i])->dl->td) > max ) {
3177: max = t;
3178: maxi = i;
3179: }
3180: if ( max < 0 ) {
3181: *h = 0;
3182: return -1;
3183: } else {
3184: m = BDY(e[maxi]);
1.66 ! noro 3185: NEWMP(mr); mr->dl = m->dl; mr->c = (Obj)ONE; NEXT(mr) = 0;
1.35 noro 3186: MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td; /* XXX */
3187: return maxi;
1.32 noro 3188: }
3189: }
1.42 noro 3190:
3191: /* return 1 if 0 <_w1 v && v <_w2 0 */
3192:
3193: int in_c12(int n,int *v,int row1,int **w1,int row2, int **w2)
3194: {
3195: int t1,t2;
3196:
3197: t1 = compare_zero(n,v,row1,w1);
3198: t2 = compare_zero(n,v,row2,w2);
3199: if ( t1 > 0 && t2 < 0 ) return 1;
3200: else return 0;
3201: }
3202:
3203: /* 0 < u => 1, 0 > u => -1 */
3204:
3205: int compare_zero(int n,int *u,int row,int **w)
3206: {
3207: int i,j,t;
3208: int *wi;
3209:
3210: for ( i = 0; i < row; i++ ) {
3211: wi = w[i];
3212: for ( j = 0, t = 0; j < n; j++ ) t += u[j]*wi[j];
3213: if ( t > 0 ) return 1;
3214: else if ( t < 0 ) return -1;
3215: }
3216: return 0;
3217: }
3218:
3219: /* functions for generic groebner walk */
3220: /* u=0 means u=-infty */
3221:
3222: int compare_facet_preorder(int n,int *u,int *v,
3223: int row1,int **w1,int row2,int **w2)
3224: {
3225: int i,j,s,t,tu,tv;
3226: int *w2i,*uv;
3227:
3228: if ( !u ) return 1;
3229: uv = W_ALLOC(n);
3230: for ( i = 0; i < row2; i++ ) {
3231: w2i = w2[i];
3232: for ( j = 0, tu = tv = 0; j < n; j++ )
3233: if ( s = w2i[j] ) {
3234: tu += s*u[j]; tv += s*v[j];
3235: }
3236: for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu;
3237: t = compare_zero(n,uv,row1,w1);
3238: if ( t > 0 ) return 1;
3239: else if ( t < 0 ) return 0;
3240: }
3241: return 1;
3242: }
3243:
1.48 noro 3244: Q inner_product_with_small_vector(VECT w,int *v)
3245: {
3246: int n,i;
3247: Q q,s,t,u;
3248:
3249: n = w->len;
3250: s = 0;
3251: for ( i = 0; i < n; i++ ) {
3252: STOQ(v[i],q); mulq((Q)w->body[i],q,&t); addq(t,s,&u); s = u;
3253: }
3254: return s;
3255: }
3256:
3257: Q compute_last_t(NODE g,NODE gh,Q t,VECT w1,VECT w2,NODE *homo,VECT *wp)
3258: {
3259: int n,i;
3260: int *wt;
3261: Q last,d1,d2,dn,nm,s,t1;
3262: VECT wd,wt1,wt2,w;
3263: NODE tg,tgh;
3264: MP f;
3265: int *h;
3266: NODE r0,r;
3267: MP m0,m;
3268: DP d;
3269:
3270: n = w1->len;
3271: wt = W_ALLOC(n);
3272: last = ONE;
3273: /* t1 = 1-t */
3274: for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
3275: f = BDY((DP)BDY(tg));
3276: h = BDY((DP)BDY(tgh))->dl->d;
3277: for ( ; f; f = NEXT(f) ) {
3278: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
3279: for ( i = 0; i < n && !wt[i]; i++ );
3280: if ( i == n ) continue;
3281: d1 = inner_product_with_small_vector(w1,wt);
3282: d2 = inner_product_with_small_vector(w2,wt);
3283: nm = d1; subq(d1,d2,&dn);
3284: /* if d1=d2 then nothing happens */
3285: if ( !dn ) continue;
3286: /* s satisfies ds = 0*/
3287: divq(nm,dn,&s);
3288:
3289: if ( cmpq(s,t) > 0 && cmpq(s,last) < 0 )
3290: last = s;
3291: else if ( !cmpq(s,t) ) {
3292: if ( cmpq(d2,0) < 0 ) {
3293: last = t;
3294: break;
3295: }
3296: }
3297: }
3298: }
3299: if ( !last ) {
3300: dn = ONE; nm = 0;
3301: } else {
3302: NTOQ(NM(last),1,nm);
3303: if ( INT(last) ) dn = ONE;
3304: else {
3305: NTOQ(DN(last),1,dn);
3306: }
3307: }
3308: /* (1-n/d)*w1+n/d*w2 -> w=(d-n)*w1+n*w2 */
3309: subq(dn,nm,&t1); mulvect(CO,(Obj)w1,(Obj)t1,(Obj *)&wt1);
3310: mulvect(CO,(Obj)w2,(Obj)nm,(Obj *)&wt2); addvect(CO,wt1,wt2,&w);
3311:
3312: r0 = 0;
3313: for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
3314: f = BDY((DP)BDY(tg));
3315: h = BDY((DP)BDY(tgh))->dl->d;
3316: for ( m0 = 0; f; f = NEXT(f) ) {
3317: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
3318: for ( i = 0; i < n && !wt[i]; i++ );
3319: if ( !inner_product_with_small_vector(w,wt) ) {
3320: NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
3321: }
3322: }
3323: NEXT(m) = 0;
3324: MKDP(((DP)BDY(tg))->nv,m0,d); d->sugar = ((DP)BDY(tg))->sugar;
3325: NEXTNODE(r0,r); BDY(r) = (pointer)d;
3326: }
3327: NEXT(r) = 0;
3328: *homo = r0;
3329: *wp = w;
3330: return last;
3331: }
3332:
1.42 noro 3333: /* return 0 if last_w = infty */
3334:
3335: NODE compute_last_w(NODE g,NODE gh,int n,int **w,
3336: int row1,int **w1,int row2,int **w2)
3337: {
3338: DP d;
3339: MP f,m0,m;
3340: int *wt,*v,*h;
3341: NODE t,s,n0,tn,n1,r0,r;
3342: int i;
3343:
3344: wt = W_ALLOC(n);
3345: n0 = 0;
3346: for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
3347: f = BDY((DP)BDY(t));
3348: h = BDY((DP)BDY(s))->dl->d;
3349: for ( ; f; f = NEXT(f) ) {
3350: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
3351: for ( i = 0; i < n && !wt[i]; i++ );
3352: if ( i == n ) continue;
3353:
3354: if ( in_c12(n,wt,row1,w1,row2,w2) &&
3355: compare_facet_preorder(n,*w,wt,row1,w1,row2,w2) ) {
3356: v = (int *)MALLOC_ATOMIC(n*sizeof(int));
3357: for ( i = 0; i < n; i++ ) v[i] = wt[i];
3358: MKNODE(n1,v,n0); n0 = n1;
3359: }
3360: }
3361: }
3362: if ( !n0 ) return 0;
3363: for ( t = n0; t; t = NEXT(t) ) {
3364: v = (int *)BDY(t);
3365: for ( s = n0; s; s = NEXT(s) )
3366: if ( !compare_facet_preorder(n,v,(int *)BDY(s),row1,w1,row2,w2) )
3367: break;
3368: if ( !s ) {
3369: *w = v;
3370: break;
3371: }
3372: }
3373: if ( !t )
3374: error("compute_last_w : cannot happen");
3375: r0 = 0;
3376: for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) {
3377: f = BDY((DP)BDY(t));
3378: h = BDY((DP)BDY(s))->dl->d;
3379: for ( m0 = 0; f; f = NEXT(f) ) {
3380: for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
3381: for ( i = 0; i < n && !wt[i]; i++ );
3382: if ( i == n ||
3383: (compare_facet_preorder(n,wt,*w,row1,w1,row2,w2)
3384: && compare_facet_preorder(n,*w,wt,row1,w1,row2,w2)) ) {
3385: NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
3386: }
3387: }
1.43 noro 3388: NEXT(m) = 0;
1.42 noro 3389: MKDP(((DP)BDY(t))->nv,m0,d); d->sugar = ((DP)BDY(t))->sugar;
3390: NEXTNODE(r0,r); BDY(r) = (pointer)d;
3391: }
3392: NEXT(r) = 0;
3393: return r0;
3394: }
1.44 noro 3395:
3396: /* compute a sufficient set of d(f)=u-v */
3397:
3398: NODE compute_essential_df(DP *g,DP *gh,int ng)
3399: {
1.45 noro 3400: int nv,i,j,k,t,lj;
3401: NODE r,r1,ri,rt,r0;
3402: MP m;
3403: MP *mj;
3404: DL di,hj,dl,dlt;
3405: int *d,*dt;
3406: LIST l;
1.44 noro 3407: Q q;
1.45 noro 3408:
3409: nv = g[0]->nv;
3410: r = 0;
3411: for ( j = 0; j < ng; j++ ) {
3412: for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ );
3413: mj = (MP *)ALLOCA(lj*sizeof(MP));
3414: for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ )
3415: mj[k] = m;
3416: for ( i = 0; i < lj; i++ ) {
3417: for ( di = mj[i]->dl, k = i+1; k < lj; k++ )
3418: if ( _dl_redble(di,mj[k]->dl,nv) ) break;
3419: if ( k < lj ) mj[i] = 0;
3420: }
3421: hj = BDY(gh[j])->dl;
3422: _NEWDL(dl,nv); d = dl->d;
3423: r0 = r;
3424: for ( i = 0; i < lj; i++ ) {
3425: if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) {
3426: for ( k = 0, t = 0; k < nv; k++ ) {
3427: d[k] = hj->d[k]-di->d[k];
3428: t += d[k];
3429: }
3430: dl->td = t;
3431: #if 1
3432: for ( rt = r0; rt; rt = NEXT(rt) ) {
3433: dlt = (DL)BDY(rt);
3434: if ( dlt->td != dl->td ) continue;
3435: for ( dt = dlt->d, k = 0; k < nv; k++ )
3436: if ( d[k] != dt[k] ) break;
3437: if ( k == nv ) break;
3438: }
3439: #else
3440: rt = 0;
3441: #endif
3442: if ( !rt ) {
3443: MKNODE(r1,dl,r); r = r1;
3444: _NEWDL(dl,nv); d = dl->d;
3445: }
1.44 noro 3446: }
3447: }
3448: }
1.45 noro 3449: for ( rt = r; rt; rt = NEXT(rt) ) {
3450: dl = (DL)BDY(rt); d = dl->d;
3451: ri = 0;
3452: for ( k = nv-1; k >= 0; k-- ) {
3453: STOQ(d[k],q);
3454: MKNODE(r1,q,ri); ri = r1;
1.44 noro 3455: }
1.45 noro 3456: MKNODE(r1,0,ri); MKLIST(l,r1);
3457: BDY(rt) = (pointer)l;
1.44 noro 3458: }
3459: return r;
3460: }
1.57 noro 3461:
3462: int comp_bits_divisible(int *a,int *b,int n)
3463: {
3464: int bpi,i,wi,bi;
3465:
3466: bpi = (sizeof(int)/sizeof(char))*8;
3467: for ( i = 0; i < n; i++ ) {
3468: wi = i/bpi; bi = i%bpi;
3469: if ( !(a[wi]&(1<<bi)) && (b[wi]&(1<<bi)) ) return 0;
3470: }
3471: return 1;
3472: }
3473:
3474: int comp_bits_lex(int *a,int *b,int n)
3475: {
3476: int bpi,i,wi,ba,bb,bi;
3477:
3478: bpi = (sizeof(int)/sizeof(char))*8;
3479: for ( i = 0; i < n; i++ ) {
3480: wi = i/bpi; bi = i%bpi;
3481: ba = (a[wi]&(1<<bi))?1:0;
3482: bb = (b[wi]&(1<<bi))?1:0;
3483: if ( ba > bb ) return 1;
3484: else if ( ba < bb ) return -1;
3485: }
3486: return 0;
3487: }
3488:
3489: NODE mono_raddec(NODE ideal)
3490: {
3491: DP p;
3492: int nv,w,i,bpi,di,c,len;
3493: int *d,*s,*u,*new;
3494: NODE t,t1,v,r,rem,prev;
3495:
3496: if( !ideal ) return 0;
3497: p = (DP)BDY(ideal);
3498: nv = NV(p);
3499: bpi = (sizeof(int)/sizeof(char))*8;
3500: w = (nv+(bpi-1))/bpi;
3501: d = p->body->dl->d;
3502: if ( !NEXT(ideal) ) {
3503: for ( t = 0, i = nv-1; i >= 0; i-- ) {
3504: if ( d[i] ) {
3505: s = (int *)CALLOC(w,sizeof(int));
3506: s[i/bpi] |= 1<<(i%bpi);
3507: MKNODE(t1,s,t);
3508: t = t1;
3509: }
3510: }
3511: return t;
3512: }
3513: rem = mono_raddec(NEXT(ideal));
3514: r = 0;
3515: len = w*sizeof(int);
3516: u = (int *)CALLOC(w,sizeof(int));
3517: for ( i = nv-1; i >= 0; i-- ) {
3518: if ( d[i] ) {
3519: for ( t = rem; t; t = NEXT(t) ) {
3520: bcopy((char *)BDY(t),(char *)u,len);
3521: u[i/bpi] |= 1<<(i%bpi);
3522: for ( v = r; v; v = NEXT(v) ) {
3523: if ( comp_bits_divisible(u,(int *)BDY(v),nv) ) break;
3524: }
3525: if ( v ) continue;
3526: for ( v = r, prev = 0; v; v = NEXT(v) ) {
3527: if ( comp_bits_divisible((int *)BDY(v),u,nv) ) {
3528: if ( prev ) NEXT(prev) = NEXT(v);
3529: else r = NEXT(r);
3530: } else prev =v;
3531: }
3532: for ( v = r, prev = 0; v; prev = v, v = NEXT(v) ) {
3533: if ( comp_bits_lex(u,(int *)BDY(v),nv) < 0 ) break;
3534: }
3535: new = (int *)CALLOC(w,sizeof(int));
3536: bcopy((char *)u,(char *)new,len);
3537: MKNODE(t1,new,v);
3538: if ( prev ) NEXT(prev) = t1;
3539: else r = t1;
3540: }
3541: }
3542: }
3543: return r;
3544: }