![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Fgfs.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to Fgfs.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM$ */
! 2:
! 3: #include "ca.h"
! 4:
! 5:
! 6: void gcdsf(VL vl,P *pa,int k,P *r)
! 7: {
! 8: P *ps,*pl,*pm,*cp;
! 9: int *cn;
! 10: DCP *ml;
! 11: Obj obj;
! 12: int i,j,l,m;
! 13: P mg,mgsf,t;
! 14: VL avl,nvl,tvl,svl;
! 15:
! 16: ps = (P *)ALLOCA(k*sizeof(P));
! 17: for ( i = 0, m = 0; i < k; i++ ) {
! 18: simp_ff((Obj)pa[i],&obj);
! 19: if ( obj )
! 20: ps[m++] = obj;
! 21: }
! 22: if ( !m ) {
! 23: *r = 0;
! 24: return;
! 25: }
! 26: if ( m == 1 ) {
! 27: *r = BDY(n0);
! 28: return;
! 29: }
! 30: pl = (P *)ALLOCA(m*sizeof(P));
! 31: ml = (DCP *)ALLOCA(m*sizeof(DCP));
! 32: for ( i = 0; i < m; i++ )
! 33: monomialfctr(vl,ps[i],&pl[i],&ml[i]);
! 34: gcdmonomial(vl,ml,m,&mg); simp_ff(mg,&mgsf);
! 35: for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) {
! 36: clctv(vl,pl[i],&tvl);
! 37: intersectv(nvl,tvl,&svl); nvl = svl;
! 38: mergev(vl,avl,tvl,&svl); avl = svl;
! 39: }
! 40: if ( !nvl ) {
! 41: *r = mgsf;
! 42: return;
! 43: }
! 44: if ( !NEXT(avl) ) {
! 45: ugcdsf(pl,m,&t);
! 46: mulp(vl,mgsf,t,r);
! 47: return;
! 48: }
! 49: for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ );
! 50: for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ );
! 51: if ( i == j ) {
! 52: /* all the pl[i]'s have the same variables */
! 53: gcdsf_main(avl,pl,m,&t);
! 54: } else {
! 55: cp = (P **)ALLOCA(m*sizeof(P *));
! 56: cn = (int *)ALLOCA(m*sizeof(int));
! 57: for ( i = 0; i < m; i++ ) {
! 58: cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P));
! 59: cn[i] = pcoef(vl,nvl,pl[i],cp[i]);
! 60: }
! 61: for ( i = j = 0; i < m; i++ )
! 62: j += cn[i];
! 63: pm = (P *)ALLOCA(j*sizeof(P));
! 64: for ( i = l = 0; i < m; i++ )
! 65: for ( j = 0; j < cn[i]; j++ )
! 66: pm[l++] = cp[i][j];
! 67: gcdsf(vl,pm,l,&t);
! 68: }
! 69: mulp(vl,mgsf,t,r);
! 70: }
! 71:
! 72: /* univariate gcd */
! 73:
! 74: void ugcdsf(P *pa,int m,P *r)
! 75: {
! 76: P *pa;
! 77: int i;
! 78: UM w1,w2,w3,w;
! 79: int d;
! 80: V v;
! 81:
! 82: if ( m == 1 ) {
! 83: *r = pa[0];
! 84: return;
! 85: }
! 86: ps = (P *)ALLOCA(m*sizeof(P));
! 87: v = VR(pa[0]);
! 88: sort_by_deg(m,pa,ps);
! 89: d = getdeg(ps[m-1],v);
! 90: w1 = W_UMALLOC(d);
! 91: w2 = W_UMALLOC(d);
! 92: w3 = W_UMALLOC(d);
! 93: ptosfum(ps[0],w1);
! 94: for ( i = 1; i < m; i++ ) {
! 95: ptosfum(ps[i],w2);
! 96: gcdsfum(w1,w2,w3);
! 97: w = w1; w1 = w3; w3 = w;
! 98: if ( !DEG(w1) ) {
! 99: MKGFS(0,*r);
! 100: return;
! 101: }
! 102: }
! 103: sfumtop(v,w1,r);
! 104: }
! 105:
! 106:
! 107: /* all the pa[i]'s have the same variables (=vl) */
! 108:
! 109: void gcdsf_main(VL vl,P *pa,int m,P *r)
! 110: {
! 111: for ( nv = 0, tvl = vl; tvl; tvl = NEXT(tvl), nv++);
! 112: if ( nv == 1 ) {
! 113: ugcdsf(pa,m,r);
! 114: return;
! 115: }
! 116: /* find v s.t. min(deg(pa[i],v)) is minimal */
! 117: tvl = vl;
! 118: do {
! 119: v = tvl->v;
! 120: i = 0;
! 121: do {
! 122: d = getdeg(pa[i],v);
! 123: if ( i == 0 || (d < d0) ) {
! 124: d0 = d; i0 = i; v0 = v;
! 125: }
! 126: } while ( ++i < m );
! 127: if ( tvl == vl || (d0 < dmin) ) {
! 128: dmin = d0; imin = i0; vmin = v0;
! 129: }
! 130: } while ( tvl = NEXT(tvl) );
! 131:
! 132: /* reorder variables so that vmin is the last variable */
! 133: for ( nvl0 = 0, rvl0 = 0, tvl = vl; tvl; tvl = NEXT(tvl) )
! 134: if ( tvl->v != vmin ) {
! 135: NEXTVL(nvl0,nvl); nvl->v = tvl->v;
! 136: NEXTVL(rvl0,rvl); rvl->v = tvl->v;
! 137: }
! 138: /* rvl = remaining variables */
! 139: NEXT(rvl) = 0;
! 140: /* nvl = ...,vmin */
! 141: NEXTVL(nvl0,nvl); nvl->v = vmin; NEXT(nvl) = 0;
! 142:
! 143: /* for content and primitive part */
! 144: pc = (P *)ALLOCA(m*sizeof(P));
! 145: ps = (P *)ALLOCA(m*sizeof(P));
! 146: ph = (P *)ALLOCA(m*sizeof(P));
! 147: /* separate the contents */
! 148: for ( i = 0; i < m; i++ ) {
! 149: reorderp(vl,nvl,pa[i],&t);
! 150: cont_pp_mv_sf(nvl,rvl,t,&pc[i],&ps[i]);
! 151: head_monomial(vmin,ps[i],&ph[i],&t);
! 152: }
! 153: ugcdsf(pc,m,&cont);
! 154: ugcdsf(ph,m,&hg);
! 155:
! 156: /* for hg*pp (used in check phase) */
! 157: lps = (P *)ALLOCA(m*sizeof(P));
! 158: for ( i = 0; i < m; i++ )
! 159: mulp(nvl,hg,ps[i],&lps[i]);
! 160:
! 161: while ( 1 ) {
! 162: g = 0;
! 163: cofmin = 0;
! 164: hm = 0;
! 165: MKGFS(0,mod);
! 166: index = 0;
! 167: while ( getdeg(mod,vmin) <= d+1 ) {
! 168: /* evaluation pt */
! 169: MKGFS(index,s);
! 170: substp(nvl,hg,vmin,s,&hge);
! 171: if ( !hge )
! 172: continue;
! 173: for ( i = 0; i < m; i++ )
! 174: substp(nvl,ps[i],vmin,s,&ph[i]);
! 175: /* ge = GCD(ps[0]|x=s,...,ps[m-1]|x=s) */
! 176: gcdsf(nvl,ph,m,&ge);
! 177: head_monomial(vmin,ge,&ce,&he);
! 178: if ( NUM(hge) ) {
! 179: *r = cont;
! 180: return;
! 181: }
! 182: divgfs(hge,ce,&t); mulp(nvl,t,ge,&u); ge = u;
! 183: divsp(nvl,ph[imin],ge,&t); mulp(nvl,hge,t,&cof1e);
! 184: if ( !hm || !cmpp(he,hm) ) {
! 185: }
! 186: }
! 187: }
! 188: }
! 189:
! 190: void head_monomial(V v,P p,P *coef,P *term)
! 191: {
! 192: P t,s,u;
! 193: DCP dc;
! 194: GFS one;
! 195:
! 196: MKGFS(0,one); t = one;
! 197: while ( 1 ) {
! 198: if ( NUM(p) || VR(p) == v ) {
! 199: *coef = p;
! 200: *term = t;
! 201: return;
! 202: } else {
! 203: NEWDC(dc); MKGFS(0,one);
! 204: COEF(dc) = one; DEG(dc) = DEG(DC(p));
! 205: MKP(VR(p),dc,s);
! 206: mulp(vl,t,s,&u); t = u;
! 207: p = COEF(DC(p));
! 208: }
! 209: }
! 210: }
! 211:
! 212: void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp)
! 213: {
! 214: DP dp;
! 215: MP t;
! 216: int i,m;
! 217: P *ps;
! 218:
! 219: ptod(vl,rvl,p,&dp);
! 220: for ( t = BDY(dp), m = 0; t; t = NEXT(t), m++ );
! 221: ps = (P *)ALLOCA(m*sizeof(P));
! 222: for ( t = BDY(dp), i = 0; t; t = NEXT(t), i++ );
! 223: ps[i] = C(t);
! 224: ugcdsf(ps,m,c);
! 225: divsp(vl,p,c,pp);
! 226: }