version 1.11, 2002/11/22 08:44:57 |
version 1.19, 2003/03/12 02:06:11 |
|
|
/* $OpenXM: OpenXM_contrib2/asir2000/engine/Fgfs.c,v 1.10 2002/11/22 07:32:10 noro Exp $ */ |
/* $OpenXM: OpenXM_contrib2/asir2000/engine/Fgfs.c,v 1.18 2003/01/17 00:17:10 noro Exp $ */ |
|
|
#include "ca.h" |
#include "ca.h" |
|
|
Line 12 void partial_sqfrsf(VL vl,V v,P f,P *r,DCP *dcp); |
|
Line 12 void partial_sqfrsf(VL vl,V v,P f,P *r,DCP *dcp); |
|
void gcdsf(VL vl,P *pa,int k,P *r); |
void gcdsf(VL vl,P *pa,int k,P *r); |
void mfctrsfmain(VL vl, P f, DCP *dcp); |
void mfctrsfmain(VL vl, P f, DCP *dcp); |
void next_evaluation_point(int *mev,int n); |
void next_evaluation_point(int *mev,int n); |
void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,P *lcp); |
void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,int *mev,P *lcp); |
void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P lcu,P lcv,P *up); |
void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P lcu,P lcv,int *mev,P *up); |
void substvp_sf(VL vl,VL rvl,P f,int *mev,P *r); |
void substvp_sf(VL vl,VL rvl,P f,int *mev,P *r); |
void shift_sf(VL vl, VL rvl, P f, int *mev, int sgn, P *r); |
void shift_sf(VL vl, VL rvl, P f, int *mev, int sgn, P *r); |
void adjust_coef_sf(VL vl,VL rvl,P lcu,P u0,P *r); |
void adjust_coef_sf(VL vl,VL rvl,P lcu,P u0,int *mev,P *r); |
void extended_gcd_modyk(P u0,P v0,V x,V y,int dy,P *cu,P *cv); |
void extended_gcd_modyk(P u0,P v0,V x,V y,int dy,P *cu,P *cv); |
void poly_to_gfsn_poly(VL vl,P f,V v,P *r); |
void poly_to_gfsn_poly(VL vl,P f,V v,P *r); |
void gfsn_poly_to_poly(VL vl,P f,V v,P *r); |
void gfsn_poly_to_poly(VL vl,P f,V v,P *r); |
Line 25 void gfsn_poly_to_poly_main(P f,V v,P *r); |
|
Line 25 void gfsn_poly_to_poly_main(P f,V v,P *r); |
|
void gfsn_univariate_to_sfbm(P f,int dy,BM *r); |
void gfsn_univariate_to_sfbm(P f,int dy,BM *r); |
void sfbm_to_gfsn_univariate(BM f,V x,V y,P *r); |
void sfbm_to_gfsn_univariate(BM f,V x,V y,P *r); |
|
|
|
void monomialfctr_sf(VL vl,P p,P *pr,DCP *dcp) |
|
{ |
|
VL nvl,avl; |
|
Q d; |
|
P f,t,s; |
|
DCP dc0,dc; |
|
Obj obj; |
|
|
|
clctv(vl,p,&nvl); |
|
for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) { |
|
getmindeg(avl->v,f,&d); |
|
if ( d ) { |
|
MKV(avl->v,t); |
|
simp_ff((Obj)t,&obj); t = (P)obj; |
|
NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t; |
|
pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t; |
|
} |
|
} |
|
if ( dc0 ) |
|
NEXT(dc) = 0; |
|
*pr = f; *dcp = dc0; |
|
} |
|
|
void lex_lc(P f,P *c) |
void lex_lc(P f,P *c) |
{ |
{ |
if ( !f || NUM(f) ) |
if ( !f || NUM(f) ) |
Line 50 void sqfrsf(VL vl, P f, DCP *dcp) |
|
Line 73 void sqfrsf(VL vl, P f, DCP *dcp) |
|
{ |
{ |
DCP dc,dct; |
DCP dc,dct; |
Obj obj; |
Obj obj; |
P t,s,c; |
P t,s,c,cont; |
VL tvl,nvl; |
VL tvl,onevl; |
|
|
simp_ff((Obj)f,&obj); f = (P)obj; |
simp_ff((Obj)f,&obj); f = (P)obj; |
lex_lc(f,&c); divsp(vl,f,c,&t); f = t; |
lex_lc(f,&c); divsp(vl,f,c,&t); f = t; |
monomialfctr(vl,f,&t,&dc); f = t; |
monomialfctr_sf(vl,f,&t,&dc); f = t; |
clctv(vl,f,&tvl); vl = tvl; |
clctv(vl,f,&tvl); vl = tvl; |
|
NEWVL(onevl); NEXT(onevl)=0; |
if ( !vl ) |
if ( !vl ) |
; |
; |
else if ( !NEXT(vl) ) { |
else if ( !NEXT(vl) ) { |
Line 65 void sqfrsf(VL vl, P f, DCP *dcp) |
|
Line 89 void sqfrsf(VL vl, P f, DCP *dcp) |
|
} else { |
} else { |
t = f; |
t = f; |
for ( tvl = vl; tvl; tvl = NEXT(tvl) ) { |
for ( tvl = vl; tvl; tvl = NEXT(tvl) ) { |
reordvar(vl,tvl->v,&nvl); |
onevl->v = tvl->v; |
cont_pp_mv_sf(vl,NEXT(nvl),t,&c,&s); t = s; |
cont_pp_mv_sf(vl,onevl,t,&cont,&s); t = s; |
sqfrsf(vl,c,&dct); |
sqfrsf(vl,cont,&dct); |
dc = append_dc(dc,NEXT(dct)); |
dc = append_dc(dc,NEXT(dct)); |
} |
} |
sqfrsfmain(vl,t,&dct); |
sqfrsfmain(vl,t,&dct); |
Line 473 void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp) |
|
Line 497 void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp) |
|
MP t; |
MP t; |
int i,m; |
int i,m; |
P *ps; |
P *ps; |
|
struct order_spec spec, currentspec; |
|
extern struct order_spec dp_current_spec; |
|
|
|
currentspec = dp_current_spec; |
|
create_order_spec(0,&spec); |
|
initd(&spec); |
ptod(vl,rvl,p,&dp); |
ptod(vl,rvl,p,&dp); |
for ( t = BDY(dp), m = 0; t; t = NEXT(t), m++ ); |
for ( t = BDY(dp), m = 0; t; t = NEXT(t), m++ ); |
ps = (P *)ALLOCA(m*sizeof(P)); |
ps = (P *)ALLOCA(m*sizeof(P)); |
Line 481 void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp) |
|
Line 510 void cont_pp_mv_sf(VL vl,VL rvl,P p,P *c,P *pp) |
|
ps[i] = C(t); |
ps[i] = C(t); |
gcdsf(vl,ps,m,c); |
gcdsf(vl,ps,m,c); |
divsp(vl,p,*c,pp); |
divsp(vl,p,*c,pp); |
|
initd(¤tspec); |
} |
} |
|
|
void mfctrsf(VL vl, P f, DCP *dcp) |
void mfctrsf(VL vl, P f, DCP *dcp) |
Line 515 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 545 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
P gcd,g,df,dfmin; |
P gcd,g,df,dfmin; |
P pa[2]; |
P pa[2]; |
P f0,pp0,spp0,c,c0,x,y,u,v,lcf,lcu,lcv,u0,v0,t,s; |
P f0,pp0,spp0,c,c0,x,y,u,v,lcf,lcu,lcv,u0,v0,t,s; |
P ype,yme; |
P ype,yme,fin; |
GFS ev,evy; |
GFS ev,evy; |
P *fp0; |
P *fp0; |
int *mev,*win; |
int *mev,*win; |
Line 531 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 561 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
return; |
return; |
} |
} |
for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ ); |
for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ ); |
va = (V *)ALLOCA(n*sizeof(int)); |
va = (V *)ALLOCA(n*sizeof(V)); |
da = (int *)ALLOCA(n*sizeof(int)); |
da = (int *)ALLOCA(n*sizeof(int)); |
/* find v s.t. diff(f,v) is nonzero and deg(f,v) is minimal */ |
/* find v s.t. diff(f,v) is nonzero and deg(f,v) is minimal */ |
imin = -1; |
imin = -1; |
Line 578 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 608 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
} |
} |
NEXT(tvl) = 0; |
NEXT(tvl) = 0; |
|
|
|
fin = f; |
reorderp(nvl,vl,f,&g); f = g; |
reorderp(nvl,vl,f,&g); f = g; |
vx = nvl->v; |
vx = nvl->v; |
vy = NEXT(nvl)->v; |
vy = NEXT(nvl)->v; |
Line 626 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 657 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
|
|
if ( !NEXT(dc) ) { |
if ( !NEXT(dc) ) { |
/* f is irreducible */ |
/* f is irreducible */ |
NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; |
NEWDC(dc); DEG(dc) = ONE; COEF(dc) = fin; NEXT(dc) = 0; |
*dcp = dc; |
*dcp = dc; |
return; |
return; |
} |
} |
Line 645 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 676 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
lcf = COEF(DC(f)); |
lcf = COEF(DC(f)); |
mfctrsf(nvl,lcf,&dct); |
mfctrsf(nvl,lcf,&dct); |
/* skip the first element (= a number) */ |
/* skip the first element (= a number) */ |
dct = NEXT(dct); |
lcfdc = NEXT(dct); |
|
|
/* shift lcfdc; c <- c(X+mev) */ |
|
for ( lcfdc = 0; dct; dct = NEXT(dct) ) { |
|
NEXTDC(lcfdc,dcs); |
|
DEG(dcs) = DEG(dct); |
|
shift_sf(nvl,rvl,COEF(dct),mev,1,&COEF(dcs)); |
|
} |
|
NEXT(dcs) = 0; |
|
|
|
/* np = number of bivariate factors */ |
/* np = number of bivariate factors */ |
for ( np = 0, dct = dc; dct; dct = NEXT(dct), np++ ); |
for ( np = 0, dct = dc; dct; dct = NEXT(dct), np++ ); |
fp0 = (P *)ALLOCA((np+1)*sizeof(P)); |
fp0 = (P *)ALLOCA((np+1)*sizeof(P)); |
Line 664 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 687 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
l = tl = (P *)ALLOCA((np+1)*sizeof(P)); |
l = tl = (P *)ALLOCA((np+1)*sizeof(P)); |
win = W_ALLOC(np+1); |
win = W_ALLOC(np+1); |
|
|
/* f <- f(X+mev) */ |
|
shift_sf(nvl,rvl,f,mev,1,&s); f = s; |
|
|
|
for ( k = 1, win[0] = 1, --np; ; ) { |
for ( k = 1, win[0] = 1, --np; ; ) { |
itogfs(1,&u0); |
itogfs(1,&u0); |
/* u0 = product of selected factors */ |
/* u0 = product of selected factors */ |
Line 675 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 695 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
} |
} |
/* we have to consider the content */ |
/* we have to consider the content */ |
/* f0 = c0*u0*v0 */ |
/* f0 = c0*u0*v0 */ |
mulp(nvl,LC(u0),c0,&c); estimatelc_sf(nvl,rvl,c,lcfdc,&lcu); |
mulp(nvl,LC(u0),c0,&c); estimatelc_sf(nvl,rvl,c,lcfdc,mev,&lcu); |
divsp(nvl,pp0,u0,&v0); |
divsp(nvl,pp0,u0,&v0); |
mulp(nvl,LC(v0),c0,&c); estimatelc_sf(nvl,rvl,c,lcfdc,&lcv); |
mulp(nvl,LC(v0),c0,&c); estimatelc_sf(nvl,rvl,c,lcfdc,mev,&lcv); |
mfctrsf_hensel(nvl,rvl,f,pp0,u0,v0,lcu,lcv,&u); |
mfctrsf_hensel(nvl,rvl,f,pp0,u0,v0,lcu,lcv,mev,&u); |
if ( u ) { |
if ( u ) { |
/* save the factor */ |
/* save the factor */ |
reorderp(vl,nvl,u,&t); |
reorderp(vl,nvl,u,&t); |
/* x -> x-mev, y -> y-evy */ |
/* y -> y-evy */ |
shift_sf(vl,rvl,t,mev,-1,&s); substp(vl,s,vy,yme,tl++); |
substp(vl,t,vy,yme,tl++); |
|
|
/* update f,pp0 */ |
/* update f,pp0 */ |
divsp(nvl,f,u,&t); f = t; |
divsp(nvl,f,u,&t); f = t; |
Line 714 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
Line 734 void mfctrsfmain(VL vl, P f, DCP *dcp) |
|
} |
} |
} |
} |
reorderp(vl,nvl,f,&t); |
reorderp(vl,nvl,f,&t); |
/* x -> x-mev, y -> y-evy */ |
/* y -> y-evy */ |
shift_sf(vl,rvl,t,mev,-1,&s); substp(vl,s,vy,yme,tl++); |
substp(vl,t,vy,yme,tl++); |
*tl = 0; |
*tl = 0; |
for ( dc0 = 0, i = 0; l[i]; i++ ) { |
for ( dc0 = 0, i = 0; l[i]; i++ ) { |
NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = l[i]; |
NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = l[i]; |
Line 742 void next_evaluation_point(int *e,int n) |
|
Line 762 void next_evaluation_point(int *e,int n) |
|
|
|
/* |
/* |
* dc : f1^E1*...*fk^Ek |
* dc : f1^E1*...*fk^Ek |
* find e1,...,ek s.t. fi(0)^ei | c |
* find e1,...,ek s.t. fi(mev)^ei | c |
* and return f1^e1*...*fk^ek |
* and return f1^e1*...*fk^ek |
* vl = (vx,vy,rvl) |
* vl = (vx,vy,rvl) |
*/ |
*/ |
|
|
void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,P *lcp) |
void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,int *mev,P *lcp) |
{ |
{ |
DCP dct; |
DCP dct; |
P r,c1,c2,t,s,f; |
P r,c1,c2,t,s,f; |
Line 758 void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,P *lcp) |
|
Line 778 void estimatelc_sf(VL vl,VL rvl,P c,DCP dc,P *lcp) |
|
if ( NUM(COEF(dct)) ) |
if ( NUM(COEF(dct)) ) |
continue; |
continue; |
/* constant part */ |
/* constant part */ |
substvp_sf(vl,rvl,COEF(dct),0,&f); |
substvp_sf(vl,rvl,COEF(dct),mev,&f); |
d = QTOS(DEG(dct)); |
d = QTOS(DEG(dct)); |
for ( i = 0, c1 = c; i < d; i++ ) |
for ( i = 0, c1 = c; i < d; i++ ) |
if ( !divtp(vl,c1,f,&c2) ) |
if ( !divtp(vl,c1,f,&c2) ) |
Line 819 void shift_sf(VL vl, VL rvl, P f, int *mev, int sgn, P |
|
Line 839 void shift_sf(VL vl, VL rvl, P f, int *mev, int sgn, P |
|
* pp(f(0)) = u0*v0 |
* pp(f(0)) = u0*v0 |
*/ |
*/ |
|
|
void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P lcu,P lcv,P *up) |
void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P lcu,P lcv,int *mev,P *up) |
{ |
{ |
VL tvl,onevl; |
VL tvl,onevl; |
P t,s,w,u,v,ff,si,wu,wv,fj,cont; |
P t,s,w,u,v,ff,si,wu,wv,fj,cont; |
Line 828 void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P |
|
Line 848 void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P |
|
int dy,n,i,dbd,nv,j; |
int dy,n,i,dbd,nv,j; |
int *md; |
int *md; |
P *uh,*vh; |
P *uh,*vh; |
P x,du0,dv0,m,q,r; |
P x,du0,dv0,m,q,r,fin; |
P *cu,*cv; |
P *cu,*cv; |
GFSN inv; |
GFSN inv; |
|
|
/* adjust coeffs */ |
/* check the validity of lc's and adjust coeffs */ |
/* u0 = am x^m+ ... -> lcu*x^m + a(m-1)*(lcu(0)/am)*x^(m-1)+... */ |
|
/* v0 = bm x^l+ ... -> lcv*x^l + b(l-1)*(lcv(0)/bl)*x^(l-1)+... */ |
|
/* f -> lcu*lcv*x^(m+l)+... */ |
/* f -> lcu*lcv*x^(m+l)+... */ |
adjust_coef_sf(vl,rvl,lcu,u0,&u); |
mulp(vl,lcu,lcv,&t); |
adjust_coef_sf(vl,rvl,lcv,v0,&v); |
if ( !divtp(vl,t,LC(f),&m) ) { |
mulp(vl,lcu,lcv,&t); divsp(vl,t,LC(f),&m); mulp(vl,m,f,&t); f = t; |
*up = 0; return; |
|
} |
|
mulp(vl,m,f,&t); f = t; |
|
/* u0 = am x^m+ ... -> lcu*x^m + a(m-1)*(lcu(mev)/am)*x^(m-1)+... */ |
|
/* v0 = bm x^l+ ... -> lcv*x^l + b(l-1)*(lcv(mev)/bl)*x^(l-1)+... */ |
|
adjust_coef_sf(vl,rvl,lcu,u0,mev,&u); |
|
adjust_coef_sf(vl,rvl,lcv,v0,mev,&v); |
|
|
|
/* f <- f(X+mev), u <- u(X+mev), v <- v(X+mev) */ |
|
fin = f; |
|
shift_sf(vl,rvl,f,mev,1,&s); f = s; |
|
shift_sf(vl,rvl,u,mev,1,&s); u = s; |
|
shift_sf(vl,rvl,v,mev,1,&s); v = s; |
|
|
vx = vl->v; vy = NEXT(vl)->v; |
vx = vl->v; vy = NEXT(vl)->v; |
n = getdeg(vx,f); |
n = getdeg(vx,f); |
dy = getdeg(vy,f)+1; |
dy = getdeg(vy,f)+1; |
Line 930 void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P |
|
Line 960 void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P |
|
} |
} |
if ( !wu ) { |
if ( !wu ) { |
gfsn_poly_to_poly(vl,u,vy,&t); |
gfsn_poly_to_poly(vl,u,vy,&t); |
if ( divtp(vl,f,t,&q) ) { |
shift_sf(vl,rvl,t,mev,-1,&s); |
cont_pp_mv_sf(vl,onevl,t,&cont,up); |
if ( divtp(vl,fin,s,&q) ) { |
|
cont_pp_mv_sf(vl,onevl,s,&cont,up); |
return; |
return; |
} |
} |
} |
} |
if ( !wv ) { |
if ( !wv ) { |
gfsn_poly_to_poly(vl,v,vy,&t); |
gfsn_poly_to_poly(vl,v,vy,&t); |
if ( divtp(vl,f,t,&q) ) { |
shift_sf(vl,rvl,t,mev,-1,&s); |
|
if ( divtp(vl,fin,s,&q) ) { |
cont_pp_mv_sf(vl,onevl,q,&cont,up); |
cont_pp_mv_sf(vl,onevl,q,&cont,up); |
return; |
return; |
} |
} |
Line 947 void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P |
|
Line 979 void mfctrsf_hensel(VL vl,VL rvl,P f,P pp0,P u0,P v0,P |
|
addp(vl,v,wv,&t); v = t; |
addp(vl,v,wv,&t); v = t; |
addp(vl,vh[j],wv,&t); vh[j] = t; |
addp(vl,vh[j],wv,&t); vh[j] = t; |
} |
} |
|
gfsn_poly_to_poly(vl,u,vy,&t); |
|
shift_sf(vl,rvl,t,mev,-1,&s); |
|
if ( divtp(vl,fin,s,&q) ) |
|
cont_pp_mv_sf(vl,onevl,s,&cont,up); |
|
else |
|
*up = 0; |
} |
} |
|
|
void adjust_coef_sf(VL vl,VL rvl,P lcu,P u0,P *r) |
void adjust_coef_sf(VL vl,VL rvl,P lcu,P u0,int *mev,P *r) |
{ |
{ |
P lcu0,cu; |
P lcu0,cu; |
DCP dc0,dcu,dc; |
DCP dc0,dcu,dc; |
|
|
substvp_sf(vl,rvl,lcu,0,&lcu0); |
substvp_sf(vl,rvl,lcu,mev,&lcu0); |
divsp(vl,lcu0,LC(u0),&cu); |
divsp(vl,lcu0,LC(u0),&cu); |
for ( dc0 = 0, dcu = DC(u0); dcu; dcu = NEXT(dcu) ) { |
for ( dc0 = 0, dcu = DC(u0); dcu; dcu = NEXT(dcu) ) { |
if ( !dc0 ) { |
if ( !dc0 ) { |