OpenXM_contrib2/asir2000/engine/Hgfs.c - diff

Return to Hgfs.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/Hgfs.c between version 1.12 and 1.26

version 1.12, 2001/06/27 04:07:57

version 1.26, 2002/10/25 02:43:40

Line 1

/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.11 2001/06/26 09:47:05 noro Exp $ */

/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.25 2002/10/23 07:54:58 noro Exp $ */

#include "ca.h"

#include "inline.h"

struct p_pair {

void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1);

UM p0;

void extractcoefbm(BM f,int dx,UM r);

UM p1;

struct p_pair *next;

};

void canzassf(UM,int,UM *);

int comp_dum(DUM a,DUM b)

void lnfsf(int,UM,UM,struct p_pair *,UM,UM);

void minipolysf(UM,UM,UM);

void czsfum(UM,UM *);

void gensqfrsfum(UM,DUM);

void sfbmtop(int,BM,V,V,P *);

void pp_sfp(VL,P,P *);

void sfcsump(VL,P,P *);

void mulsfbmarray(int,BM,ML,int,int *,V,V,P *);

void const_term(P,UM);

int comp_dum(a,b)

DUM a,b;

{

if ( DEG(a->f) > DEG(b->f) )

return -1;

Line 30 DUM a,b;

Line 16 DUM a,b;

return 0;

}

void fctrsf(p,dcp)

void ufctrsf(P p,DCP *dcp)

P p;

DCP *dcp;

{

int n,i,j,k;

DCP dc,dc0;

P lc;

P zp;

UM mp;

UM *tl;

Obj obj;

struct oDUM *udc,*udc1;

simp_ff(p,&zp); p = zp;

simp_ff((Obj)p,&obj); p = (P)obj;

if ( !p ) {

*dcp = 0; return;

NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE;

NEXT(dc) = 0; *dcp = dc;

return;

}

mp = W_UMALLOC(UDEG(p));

ptosfum(p,mp);

Line 89 DCP *dcp;

Line 75 DCP *dcp;

NEXT(dc) = 0; *dcp = dc0;

}

void gensqfrsfum(p,dc)

void gensqfrsfum(UM p,struct oDUM *dc)

UM p;

struct oDUM *dc;

{

int n,i,j,d,mod;

UM t,s,g,f,f1,b;

Line 144 struct oDUM *dc;

Line 128 struct oDUM *dc;

}

void randsfum(d,p)

void randsfum(int d,UM p)

int d;

UM p;

{

int i;

Line 156 UM p;

Line 138 UM p;

p->d = i;

}

void pwrmodsfum(p,e,f,pr)

void pwrmodsfum(UM p,int e,UM f,UM pr)

int e;

UM p,f,pr;

{

UM wt,ws,q;

Line 186 UM p,f,pr;

Line 166 UM p,f,pr;

}

void spwrsfum(m,f,e,r)

void spwrsfum(UM m,UM f,N e,UM r)

UM f,m,r;

N e;

{

UM t,s,q;

N e1;

Line 210 N e;

Line 188 N e;

}

void tracemodsfum(m,f,e,r)

void tracemodsfum(UM m,UM f,int e,UM r)

UM f,m,r;

int e;

{

UM t,s,q,u;

int i;

Line 232 int e;

Line 208 int e;

cpyum(s,r);

}

void make_qmatsf(p,tab,mp)

void make_qmatsf(UM p,UM *tab,int ***mp)

UM p;

UM *tab;

int ***mp;

{

int n,i,j;

int *c;

Line 256 int ***mp;

Line 229 int ***mp;

mat[i][i] = _subsf(mat[i][i],one);

}

void nullsf(mat,n,ind)

void nullsf(int **mat,int n,int *ind)

int **mat;

int *ind;

int n;

{

int i,j,l,s,h,inv;

int *t,*u;

Line 291 int n;

Line 261 int n;

}

void null_to_solsf(mat,ind,n,r)

void null_to_solsf(int **mat,int *ind,int n,UM *r)

int **mat;

int *ind;

int n;

UM *r;

{

int i,j,k,l;

int *c;

Line 324 nullsf(mat,n,ind)

Line 290 nullsf(mat,n,ind)

null_to_solsf(ind,n,r)

void czsfum(f,r)

void czsfum(UM f,UM *r)

UM f,*r;

{

int i,j;

int d,n,ord;

Line 379 UM f,*r;

Line 344 UM f,*r;

r[j] = 0;

}

int berlekampsf(p,df,tab,r)

int berlekampsf(UM p,int df,UM *tab,UM *r)

UM p;

int df;

UM *tab,*r;

{

int n,i,j,k,nf,d,nr;

int **mat;

Line 428 UM *tab,*r;

Line 390 UM *tab,*r;

}

/* NOT REACHED */

error("berlekampsf : cannot happen");

return 0;

}

void minipolysf(f,p,mp)

void minipolysf(UM f,UM p,UM mp)

UM f,p,mp;

{

struct p_pair *list,*l,*l1,*lprev;

int n,d;

Line 466 UM f,p,mp;

Line 430 UM f,p,mp;

}

void lnfsf(n,p0,p1,list,np0,np1)

void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1)

int n;

UM p0,p1;

struct p_pair *list;

UM np0,np1;

{

int inv,h,d1;

int h,d1;

UM t0,t1,s0,s1;

struct p_pair *l;

Line 489 UM np0,np1;

Line 449 UM np0,np1;

}

int find_rootsf(p,root)

int find_rootsf(UM p,int *root)

UM p;

int *root;

{

UM *r;

int i,j,n;

int i,n;

n = DEG(p);

r = ALLOCA((DEG(p))*sizeof(UM));

Line 504 int *root;

Line 462 int *root;

return n;

}

void canzassf(f,d,r)

void canzassf(UM f,int d,UM *r)

UM f,*r;

int d;

{

UM t,s,u,w,g,o;

N n1,n2,n3,n4,n5;

UM *b;

int n,m,i,q,ed;

int n,q,ed;

if ( DEG(f) == d ) {

r[0] = UMALLOC(d); cpyum(f,r[0]);

Line 552 int d;

Line 508 int d;

int sfberle(VL,P,int,GFS *,DCP *);

void sfgcdgen(P,ML,ML *);

void sfhenmain(LUM,ML,ML,ML *);

void ptosflum(int,P,LUM);

void sfhenmain2(BM,UM,UM,int,BM *);

void ptosfbm(int,P,BM);

/* f = f(x,y) */

void sfhensel(count,f,x,evp,sfp,listp)

void sfhensel(int count,P f,V x,int degbound,GFS *evp,P *sfp,ML *listp)

int count;

P f;

V x;

GFS *evp;

P *sfp;

ML *listp;

{

int i,j;

int i;

int fn,n,bound;

int fn;

ML rlist;

BM fl;

VL vl,nvl;

V y;

int dx,dy,mev;

int dx,dy,bound;

GFS ev;

P f1,t,yev,c,sf;

P f1,t,c,sf;

DCP dc;

DCP dc,dct,dc0;

UM w,w1,q,fm,hm;

UM q,fm,hm;

UM *gm;

struct oEGT tmp0,tmp1,eg_hensel,eg_hensel_t;

Line 599 ML *listp;

Line 547 ML *listp;

*listp = rlist = MLALLOC(1); rlist->n = fn; rlist->c[0] = 0;

return;

}

/* pass the the leading coeff. to the first element */

if ( degbound >= 0 ) {

c = dc->c; dc = NEXT(dc);

mulp(vl,dc->c,c,&t); dc->c = t;

* reconstruct dc so that

* dc[1],... : factors satisfy degree bound

* dc[0] : product of others

c = dc->c; dc = NEXT(dc);

dc0 = 0;

fn = 0;

while ( dc ) {

if ( getdeg(x,COEF(dc)) <= degbound ) {

dct = NEXT(dc); NEXT(dc) = dc0; dc0 = dc; dc = dct;

fn++;

} else {

mulp(vl,COEF(dc),c,&t); c = t;

dc = NEXT(dc);

}

if ( OID(c) == O_P ) {

NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = dc0;

fn++;

} else {

mulp(vl,dc0->c,c,&t); dc0->c = t; dc = dc0;

}

} else {

/* pass the the leading coeff. to the first element */

c = dc->c; dc = NEXT(dc);

mulp(vl,dc->c,c,&t); dc->c = t;

}

/* convert mod y-a factors into UM */

gm = (UM *)ALLOCA(fn*sizeof(UM));

Line 610 ML *listp;

Line 584 ML *listp;

ptosfum(dc->c,gm[i]);

}

bound = dy+1;

/* set bound */

/* g | f, lc_y(g) = lc_y(f) => deg_y(g) <= deg_y(f) */

/* so, bound = dy is sufficient, but we use slightly large value */

bound = dy+2;

/* f(x,y) -> f(x,y+ev) */

fl = BMALLOC(dx,bound);

ptosfbm(bound,f,fl);

shiftsfbm(bound,fl,FTOIF(CONT(ev)));

if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));

/* sf = f(x+ev) */

sfbmtop(bound,fl,x,y,&sf);

sfbmtop(fl,x,y,&sf);

/* fm = fl mod y */

fm = W_UMALLOC(dx);

Line 653 ML *listp;

Line 631 ML *listp;

/* y -> y-a */

mev = _chsgnsf(FTOIF(CONT(ev)));

for ( i = 0; i < fn; i++ )

shiftsfbm(bound,(BM)(rlist->c[i]),mev);

shiftsfbm((BM)(rlist->c[i]),mev);

#endif

*evp = ev;

*sfp = sf;

Line 662 ML *listp;

Line 640 ML *listp;

/* main variable of f = x */

int sfberle(vl,f,count,ev,dcp)

int sfberle(VL vl,P f,int count,GFS *ev,DCP *dcp)

VL vl;

P f;

int count;

GFS *ev;

DCP *dcp;

{

UM wf,wf1,wf2,wfs,gcd;

ML flist;

int fn,n;

int fn,fn1,n;

GFS m,fm;

DCP dc,dct,dc0;

VL nvl;

V x,y;

P g,lc,lc0,f0;

P lc,lc0,f0;

int j,q1,index,i;

Obj obj;

int j,q,index,i;

clctv(vl,f,&nvl); vl = nvl;

x = vl->v; y = vl->next->v;

simp_ff(f,&g); g = f;

simp_ff((Obj)f,&obj); f = (P)obj;

n = QTOS(DEG(DC(f)));

wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n);

wfs = W_UMALLOC(n); gcd = W_UMALLOC(n);

q1 = field_order_sf()-1;

q = field_order_sf();

lc = DC(f)->c;

for ( j = 0, fn = n + 1, index = 0;

index < q1 && j < count && fn > 1; index++ ) {

index < q && j < count && fn > 1; index++ ) {

MKGFS(index,m);

indextogfs(index,&m);

substp(vl,lc,y,(P)m,&lc0);

if ( lc0 ) {

substp(vl,f,y,(P)m,&f0);

ptosfum(f0,wf); cpyum(wf,wf1);

diffsfum(wf1,wf2); gcdsfum(wf1,wf2,gcd);

if ( DEG(gcd) == 0 ) {

fctrsf(f0,&dc);

ufctrsf(f0,&dc);

for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ );

if ( i < fn ) {

dc0 = dc; fn = i; fm = m;

Line 705 DCP *dcp;

Line 678 DCP *dcp;

}

if ( index == q1 )

if ( index == q )

return 0;

else if ( fn == 1 )

return 1;

Line 716 DCP *dcp;

Line 689 DCP *dcp;

}

void sfgcdgen(f,blist,clistp)

void sfgcdgen(P f,ML blist,ML *clistp)

P f;

ML blist,*clistp;

{

int i;

int n,d,np;

Line 745 ML blist,*clistp;

Line 716 ML blist,*clistp;

}

/* f = g0*h0 mod y -> f = gk*hk mod y^(dy+1), f is replaced by hk */

sfhenmain(fl,bqlist,cqlist,listp)

void sfhenmain(f,bqlist,cqlist,listp)

void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp)

LUM f;

ML bqlist,cqlist,*listp;

{

int i,j,k;

int i,k;

int *px,*py;

int dx;

int **pp,**pp1;

UM wt,wa,wb,q,w1,w2,wh1,wg1,ws;

int n,np,bound,dr,tmp;

UM wt,wq0,wq,wr,wm,wm0,wa,q;

LUM wb0,wb1,tlum;

UM *b,*c;

LUM *l;

ML list;

n = DEG(f); np = bqlist->n; bound = bqlist->bound;

*listp = list = MLALLOC(n);

list->n = np; list->bound = bound;

W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);

wt = W_UMALLOC(n); wq0 = W_UMALLOC(n); wq = W_UMALLOC(n);

wr = W_UMALLOC(n); wm = W_UMALLOC(2*n); wm0 = W_UMALLOC(2*n);

wa = W_UMALLOC(2*n); q = W_UMALLOC(2*n);

b = (UM *)bqlist->c; c = (UM *)cqlist->c; l = (LUM *)list->c;

for ( i = 0; i < np; i++ ) {

l[i] = LUMALLOC(DEG(b[i]),bound);

for ( j = DEG(b[i]), pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- )

pp[j][0] = px[j];

}

for ( i = 1; i < bound; i++ ) {

fprintf(stderr,".");

/* at this point, f = l[0]*l[1]*...*l[np-1] mod y^i */

mulsflum(i+1,l[0],l[1],wb0);

for ( j = 2; j < np; j++ ) {

mulsflum(i+1,l[j],wb0,wb1);

tlum = wb0; wb0 = wb1; wb1 = tlum;

}

#if 0

/* check */

for ( j = 0, pp = COEF(f), pp1 = COEF(wb0); j <= n; j++ )

for ( k = 0; k < i; k++ )

if ( pp[j][k] != pp1[j][k] )

error("henmain : cannot happen");

#endif

for ( j = n, px = COEF(wt); j >= 0; j-- )

px[j] = 0;

for ( j = n, pp = COEF(f), pp1 = COEF(wb0); j >= 0; j-- )

COEF(wt)[j] = _subsf(pp[j][i],pp1[j][i]);

degum(wt,n);

for ( j = n, px = COEF(wq0); j >= 0; j-- )

px[j] = 0;

for ( j = 1; j < np; j++ ) {

mulsfum(wt,c[j],wm); dr = divsfum(wm,b[j],q);

for ( k = DEG(q), px = COEF(wq0), py = COEF(q); k >= 0; k-- )

px[k] = _addsf(px[k],py[k]);

for ( k = dr, pp = COEF(l[j]), px = COEF(wm); k >= 0; k-- )

pp[k][i] = px[k];

}

degum(wq0,n); mulsfum(wq0,b[0],wm);

mulsfum(wt,c[0],wm0); addsfum(wm,wm0,wa);

for ( j = DEG(wa), pp = COEF(l[0]), px = COEF(wa); j >= 0; j-- )

pp[j][i] = px[j];

for ( j = n, px = COEF(wq0); j >= 0; j-- )

px[j] = 0;

}

fprintf(stderr,"\n");

}

/* f = g0*h0 mod y -> f = gk*hk mod y^bound, f is replaced by hk */

void sfhenmain2(f,g0,h0,bound,gp)

BM f;

UM g0,h0;

int bound;

BM *gp;

{

int i,j,k,l;

int *px,*py;

int **pp,**pp1;

int n,np,dr,tmp;

UM wt,wa,wb,wq,wm,q,w1,w2,wh1,wg1,ws;

UM wc,wd,we,wz;

BM wb0,wb1;

int ng,nh;

int dg,dh;

BM fk,gk,hk;

n = degsfbm(bound,f);

if ( DEG(f) < dy )

ng = g0->d;

error("sfhenmain2 : invalid input");

nh = h0->d;

W_BMALLOC(n,bound,wb0); W_BMALLOC(n,bound,wb1);

dx = degbm(f);

wt = W_UMALLOC(n); ws = W_UMALLOC(n);

dg = DEG(g0);

wq = W_UMALLOC(n); q = W_UMALLOC(2*n);

dh = DEG(h0);

wg1 = W_UMALLOC(2*n); wh1 = W_UMALLOC(2*n);

W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1);

wt = W_UMALLOC(dx); ws = W_UMALLOC(dx); q = W_UMALLOC(2*dx);

wg1 = W_UMALLOC(2*dx); wh1 = W_UMALLOC(2*dx);

/* fk = gk*hk mod y^k */

W_BMALLOC(n,bound,fk);

W_BMALLOC(dx,dy,fk);

cpyum(COEF(f)[0],COEF(fk)[0]);

gk = BMALLOC(ng,bound);

gk = BMALLOC(dg,dy);

cpyum(g0,COEF(gk)[0]);

W_BMALLOC(nh,bound,hk);

W_BMALLOC(dh,dy,hk);

cpyum(h0,COEF(hk)[0]);

wc = W_UMALLOC(2*n); wd = W_UMALLOC(2*n);

wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx);

we = W_UMALLOC(2*n); wz = W_UMALLOC(2*n);

we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx);

/* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */

w1 = W_UMALLOC(ng); cpyum(g0,w1);

w1 = W_UMALLOC(dg); cpyum(g0,w1);

w2 = W_UMALLOC(nh); cpyum(h0,w2);

w2 = W_UMALLOC(dh); cpyum(h0,w2);

wa = W_UMALLOC(2*n); wb = W_UMALLOC(2*n); /* XXX */

wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */

eucsfum(w1,w2,wa,wb);

#if 0

fprintf(stderr,"dy=%d\n",dy);

mulsfum(wa,g0,wc); mulsfum(wb,h0,wd); addsfum(wc,wd,we);

for ( k = 1; k <= dy; k++ ) {

if ( DEG(we) != 0 || COEF(we)[0] != _onesf() )

error("henmain2 : cannot happen(euc)");

#endif

fprintf(stderr,"bound=%d\n",bound);

for ( k = 1; k < bound; k++ ) {

fprintf(stderr,".");

/* at this point, f = gk*hk mod y^k */

#if 0

for ( j = 0; j < k; j++ )

if ( !isequalum(COEF(f)[j],COEF(fk)[j]) )

error("henmain2 : cannot happen(history)");

#endif

/* clear wt */

bzero(COEF(wt),(n+1)*sizeof(int));

clearum(wt,dx);

/* wt = (f-gk*hk)/y^k */

subsfum(COEF(f)[k],COEF(fk)[k],wt);

/* clear wq */

bzero(COEF(wq),(n+1)*sizeof(int));

/* compute wf1,wg1 s.t. wh1*g0+wg1*h0 = wt */

mulsfum(wa,wt,wh1); DEG(wh1) = divsfum(wh1,h0,q);

mulsfum(wh1,g0,wc); subsfum(wt,wc,wd); DEG(wd) = divsfum(wd,h0,wg1);

Line 900 BM *gp;

Line 783 BM *gp;

error("henmain2 : cannot happen");

#endif

/* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod y^bound */

/* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod y^(dy+1) */

/* wb0 = wh1*y^k */

clearsfbm(bound,n,wb0);

clearbm(dx,wb0);

DEG(wb0) = bound;

cpyum(wh1,COEF(wb0)[k]);

/* wb1 = gk*wb0 mod y^bound */

/* wb1 = gk*wb0 mod y^(dy+1) */

clearsfbm(bound,n,wb1);

clearbm(dx,wb1);

mulsfbm(bound,gk,wb0,wb1);

mulsfbm(gk,wb0,wb1);

/* fk += wb1 */

addtosfbm(bound,wb1,fk);

addtosfbm(wb1,fk);

/* wb0 = wg1*y^k */

clearsfbm(bound,n,wb0);

clearbm(dx,wb0);

DEG(wb0) = bound;

cpyum(wg1,COEF(wb0)[k]);

/* wb1 = hk*wb0 mod y^bound */

/* wb1 = hk*wb0 mod y^(dy+1) */

clearsfbm(bound,n,wb1);

clearbm(dx,wb1);

mulsfbm(bound,hk,wb0,wb1);

mulsfbm(hk,wb0,wb1);

/* fk += wb1 */

addtosfbm(bound,wb1,fk);

addtosfbm(wb1,fk);

/* fk += wg1*wh1*y^(2*k) mod y^bound */

/* fk += wg1*wh1*y^(2*k) mod y^(dy+1) */

if ( 2*k < bound ) {

if ( 2*k <= dy ) {

mulsfum(wg1,wh1,wt); addsfum(COEF(fk)[2*k],wt,ws);

cpyum(ws,COEF(fk)[2*k]);

}

Line 935 BM *gp;

Line 816 BM *gp;

}

fprintf(stderr,"\n");

*gp = gk;

DEG(f) = bound;

DEG(f) = dy;

for ( i = 0; i < bound; i++ )

for ( i = 0; i <= dy; i++ )

cpyum(COEF(hk)[i],COEF(f)[i]);

}

void ptosflum(bound,f,fl)

int bound;

P f;

LUM fl;

{

DCP dc;

int **pp;

int d;

UM t;

t = UMALLOC(bound);

for ( dc = DC(f), pp = COEF(fl); dc; dc = NEXT(dc) ) {

d = QTOS(DEG(dc));

ptosfum(COEF(dc),t);

bcopy(t->c,pp[d],(t->d+1)*sizeof(int));

}

/* fl->c[i] = coef_y(f,i) */

void ptosfbm(bound,f,fl)

void ptosfbm(int dy,P f,BM fl)

int bound;

P f;

BM fl;

{

DCP dc;

int d,i,n;

int d,i,dx;

UM t;

n = QTOS(DEG(DC(f)));

dx = QTOS(DEG(DC(f)));

clearsfbm(bound,n,fl);

if ( DEG(fl) < dy )

DEG(fl) = bound;

error("ptosfbm : invalid input");

t = UMALLOC(bound);

DEG(fl) = dy;

clearbm(dx,fl);

t = UMALLOC(dy);

for ( dc = DC(f); dc; dc = NEXT(dc) ) {

d = QTOS(DEG(dc));

ptosfum(COEF(dc),t);

for ( i = 0; i <= DEG(t); i++ )

COEF(COEF(fl)[i])[d] = COEF(t)[i];

}

for ( i = 0; i < bound; i++ )

for ( i = 0; i <= dy; i++ )

degum(COEF(fl)[i],n);

degum(COEF(fl)[i],dx);

}

/* x : main variable */

void sflumtop(bound,fl,x,y,fp)

void sfbmtop(BM f,V x,V y,P *fp)

int bound;

LUM fl;

V x,y;

P *fp;

{

int i,j,n;

int **c;

UM w;

int *coef;

DCP dc,dct;

n = fl->d;

c = fl->c;

w = W_UMALLOC(bound);

for ( i = 0, dc = 0; i <= n; i++ ) {

coef = c[i];

for ( j = bound-1; j >= 0 && coef[j] == 0; j-- );

if ( j < 0 )

continue;

DEG(w) = j; bcopy(coef,COEF(w),(j+1)*sizeof(int));

NEWDC(dct); STOQ(i,DEG(dct));

sfumtop(y,w,&COEF(dct)); NEXT(dct) = dc; dc = dct;

}

MKP(x,dc,*fp);

}

/* x : main variable */

void sfbmtop(bound,f,x,y,fp)

int bound;

BM f;

V x,y;

P *fp;

{

UM *c;

int i,j,d,a;

int i,j,d,a,dy;

GFS b;

DCP dc0,dc,dct;

dy = DEG(f);

c = COEF(f);

d = DEG(c[0]);

d = degbm(f);

for ( i = 1; i < bound; i++ )

d = MAX(DEG(c[i]),d);

dc0 = 0;

for ( i = 0; i <= d; i++ ) {

dc = 0;

for ( j = 0; j < bound; j++ ) {

for ( j = 0; j <= dy; j++ ) {

if ( DEG(c[j]) >= i && (a = COEF(c[j])[i]) ) {

NEWDC(dct);

STOQ(j,DEG(dct));

MKGFS(IFTOF(a),b);

iftogfs(a,&b);

COEF(dct) = (P)b;

NEXT(dct) = dc;

dc = dct;

Line 1057 P *fp;

Line 885 P *fp;

*fp = 0;

}

void sfsqfr(P f,DCP *dcp)

{

Obj obj;

DCP dc;

VL vl;

simp_ff((Obj)f,&obj); f = (P)obj;

clctv(CO,f,&vl);

if ( !vl ) {

/* f is a const */

NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; *dcp = dc;

} else if ( !NEXT(vl) )

sfusqfr(f,dcp);

else

sqfrsf(f,dcp);

}

void sfusqfr(P f,DCP *dcp)

{

DCP dc,dct;

struct oDUM *udc;

V x;

P lc;

int n,i;

UM mf;

x = VR(f);

n = getdeg(x,f);

mf = W_UMALLOC(n);

ptosfum(f,mf);

lc = COEF(DC(f));

if ( !_isonesf(COEF(mf)[n]) ) {

monicsfum(mf);

}

W_CALLOC(n+1,struct oDUM,udc);

gensqfrsfum(mf,udc);

for ( i = 0, dc = 0; udc[i].f; i++ ) {

NEWDC(dct); STOQ(udc[i].n,DEG(dct));

sfumtop(x,udc[i].f,&COEF(dct));

NEXT(dct) = dc; dc = dct;

}

NEWDC(dct); DEG(dct) = ONE; COEF(dct) = (P)lc; NEXT(dct) = dc;

*dcp = dct;

}

#if 0

void sfbsqfrmain(P f,V x,V y,DCP *dcp)

{

/* XXX*/

}

/* f is bivariate */

void sfbsqfr(P f,V x,V y,DCP *dcp)

{

P t,rf,cx,cy;

VL vl,rvl;

DCP dcx,dcy,dct,dc;

struct oVL vl0,vl1;

/* cy(y) = cont(f,x), f /= cy */

cont_pp_sfp(vl,f,&cy,&t); f = t;

/* rvl = [y,x] */

reordvar(vl,y,&rvl); reorderp(rvl,vl,f,&rf);

/* cx(x) = cont(rf,y), Rf /= cy */

cont_pp_sfp(rvl,rf,&cx,&t); rf = t;

reorderp(vl,rvl,rf,&f);

/* f -> cx*cy*f */

sfsqfr(cx,&dcx); dcx = NEXT(dcx);

sfsqfr(cy,&dcy); dcy = NEXT(dcy);

if ( dcx ) {

for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );

NEXT(dct) = dcy;

} else

dcx = dcy;

if ( OID(f) == O_N )

*dcp = dcx;

else {

/* f must be bivariate */

sfbsqfrmain(f,x,y,&dc);

if ( dcx ) {

for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );

NEXT(dct) = dc;

} else

dcx = dc;

*dcp = dcx;

}

#endif

void sfdtest(P,ML,V,V,DCP *);

void sfbfctr(f,x,y,dcp)

/* if degbound >= 0 find factor s.t. deg_x(factor) <= degbound */

P f;

V x,y;

void sfbfctr(P f,V x,V y,int degbound,DCP *dcp)

DCP *dcp;

{

ML list;

P sf;

GFS ev;

DCP dc,dct;

BM fl;

int n,bound;

int dx,dy;

/* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */

sfhensel(5,f,x,&ev,&sf,&list);

sfhensel(5,f,x,degbound,&ev,&sf,&list);

if ( list->n == 0 )

error("sfbfctr : short of evaluation points");

else if ( list->n == 1 ) {

/* f is irreducible */

NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;

*dcp = dc;

return;

}

sfdtest(sf,list,x,y,&dc);

n = getdeg(x,sf);

if ( ev ) {

bound = list->bound;

dx = getdeg(x,sf);

W_BMALLOC(n,bound,fl);

dy = getdeg(y,sf);

for ( dct = dc; dct; dct = NEXT(dct) ) {

W_BMALLOC(dx,dy,fl);

ptosfbm(bound,COEF(dct),fl);

for ( dct = dc; dct; dct = NEXT(dct) ) {

shiftsfbm(bound,fl,_chsgnsf(FTOIF(CONT(ev))));

ptosfbm(dy,COEF(dct),fl);

sfbmtop(bound,fl,x,y,&COEF(dct));

shiftsfbm(fl,_chsgnsf(FTOIF(CONT(ev))));

sfbmtop(fl,x,y,&COEF(dct));

}

*dcp = dc;

}

/* f = f(x,y) = list->c[0]*list->c[1]*... mod y^list->bound */

/* returns shifted f, shifted factors and the eval pt */

void sfdtest(f,list,x,y,dcp)

void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P *sfp,DCP *dcp)

P f;

ML list;

V x,y;

DCP *dcp;

{

int n,np,bound;

ML list;

int i,j,k;

P sf;

GFS ev;

DCP dc,dct;

int dx,dy;

/* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */

sfhensel(5,f,x,degbound,&ev,&sf,&list);

if ( list->n == 0 )

error("sfbfctr_shift : short of evaluation points");

else if ( list->n == 1 ) {

/* f is irreducible */

NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;

*evp = 0;

*sfp = f;

*dcp = dc;

} else {

sfdtest(sf,list,x,y,dcp);

*evp = ev;

*sfp = sf;

}

/* f = f(x,y) = list->c[0]*list->c[1]*... mod y^(list->bound+1) */

void sfdtest(P f,ML list,V x,V y,DCP *dcp)

{

int np,dx,dy;

int i,j,k,bound;

int *win;

P g,lcg,factor,cofactor,lcyx;

P t,csum;

P csum;

DCP dcf,dcf0,dc;

BM *c;

BM lcy;

UM lcg0;

UM lcg0,lcy0,w;

UM *d1c;

ML wlist;

struct oVL vl1,vl0;

VL vl;

int z;

int z,dt,dtok;

/* vl = [x,y] */

vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0;

/* setup various structures and arrays */

n = UDEG(f); np = list->n; bound = list->bound; win = W_ALLOC(np+1);

dx = getdeg(x,f);

wlist = W_MLALLOC(np); wlist->n = list->n;

dy = getdeg(y,f);

wlist->bound = list->bound;

np = list->n;

win = W_ALLOC(np+1);

wlist = W_MLALLOC(np);

wlist->n = list->n;

bound = wlist->bound = list->bound;

c = (BM *)COEF(wlist);

bcopy((char *)COEF(list),(char *)c,(int)(sizeof(BM)*np));

lcg0 = W_UMALLOC(2*bound);

lcg0 = W_UMALLOC(2*dy);

/* initialize g by f */

g = f;

Line 1132 DCP *dcp;

Line 1088 DCP *dcp;

sfcsump(vl,lcg,&csum);

/* initialize lcy by LC(f) */

W_BMALLOC(0,bound,lcy);

W_BMALLOC(0,dy,lcy);

NEWDC(dc); COEF(dc) = COEF(DC(g)); DEG(dc) = 0;

NEWP(lcyx); VR(lcyx) = x; DC(lcyx) = dc;

ptosfbm(bound,lcyx,lcy);

ptosfbm(dy,lcyx,lcy);

/* initialize lcy0 by LC(f) */

lcy0 = W_UMALLOC(bound);

ptosfum(COEF(DC(g)),lcy0);

/* ((d-1 coefs)*lcy0 */

d1c = (UM *)W_ALLOC(np*sizeof(UM));

w = W_UMALLOC(2*bound);

for ( i = 1; i < np; i++ ) {

extractcoefbm(c[i],degbm(c[i])-1,w);

d1c[i] = W_UMALLOC(2*bound);

mulsfum(w,lcy0,d1c[i]);

/* d1c[i] = d1c[i] mod y^(bound+1) */

if ( DEG(d1c[i]) > bound ) {

for ( j = DEG(d1c[i]); j > bound; j-- )

COEF(d1c[i])[j] = 0;

degum(d1c[i],bound);

}

fprintf(stderr,"np = %d\n",np);

dtok = 0;

for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; z++ ) {

if ( !(z % 1000) ) fprintf(stderr,".");

if ( sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,k,win,&factor,&cofactor) ) {

dt = sfdegtest(dy,bound,d1c,k,win);

if ( dt )

dtok++;

if ( dt && sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,

k,win,&factor,&cofactor) ) {

NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;

g = cofactor;

Line 1153 DCP *dcp;

Line 1133 DCP *dcp;

/* update csum */

sfcsump(vl,lcg,&csum);

/* update dy */

dy = getdeg(y,g);

/* update lcy */

clearsfbm(bound,0,lcy);

clearbm(0,lcy);

COEF(dc) = COEF(DC(g)); ptosfbm(bound,lcyx,lcy);

COEF(dc) = COEF(DC(g));

ptosfbm(dy,lcyx,lcy);

for ( i = 0; i < k - 1; i++ )

for ( j = win[i] + 1; j < win[i + 1]; j++ )

Line 1173 DCP *dcp;

Line 1157 DCP *dcp;

else

for ( i = 1; i < k; i++ )

win[i] = win[0] + i;

/* update lcy0 */

ptosfum(COEF(DC(g)),lcy0);

/* update d-1 coeffs */

for ( i = 1; i <= np; i++ ) {

extractcoefbm(c[i],degbm(c[i])-1,w);

mulsfum(w,lcy0,d1c[i]);

/* d1c[i] = d1c[1] mod y^(bound+1) */

if ( DEG(d1c[i]) > bound ) {

for ( j = DEG(d1c[i]); j > bound; j-- )

COEF(d1c[i])[j] = 0;

degum(d1c[i],bound);

}

} else if ( !ncombi(1,np,k,win) )

if ( k == np )

break;

Line 1180 DCP *dcp;

for ( i = 0, ++k; i < k; i++ )

win[i] = i + 1;

}

fprintf(stderr,"\n");

fprintf(stderr,"total %d, omitted by degtest %d\n",z,z-dtok);

NEXTDC(dcf0,dcf); COEF(dcf) = g;

DEG(dcf) = ONE; NEXT(dcf) = 0; *dcp = dcf0;

}

void extractcoefbm(BM f,int dx,UM r)

{

int j;

UM fj;

for ( j = DEG(f); j >= 0; j-- ) {

fj = COEF(f)[j];

if ( fj && DEG(fj) >= dx ) {

COEF(r)[j] = COEF(fj)[dx];

} else

COEF(r)[j] = 0;

}

degum(r,DEG(f));

}

/* deg_y(prod mod y^(bound+1)) <= dy ? */

int sfdegtest(int dy,int bound,UM *d1c,int k,int *in)

{

int i,j;

UM w,w1,wt;

BM t;

w = W_UMALLOC(bound);

w1 = W_UMALLOC(bound);

clearum(w,bound);

for ( i = 0; i < k; i++ ) {

addsfum(w,d1c[in[i]],w1); wt = w; w = w1; w1 = wt;

}

return DEG(w) <= dy ? 1 : 0;

}

/* lcy = LC(g), lcg = lcy*g, lcg0 = const part of lcg */

int sfdtestmain(vl,lcg,lcg0,lcy,csum,list,k,in,fp,cofp)

VL vl;

int sfdtestmain(VL vl,P lcg,UM lcg0,BM lcy,P csum,ML list,

P lcg;

int k,int *in,P *fp,P *cofp)

UM lcg0;

BM lcy;

P csum;

ML list;

int k;

int *in;

P *fp,*cofp;

{

P fmul,csumg,q;

P fmul,csumg,q,cont;

V x,y;

x = vl->v;

Line 1213 P *fp,*cofp;

Line 1238 P *fp,*cofp;

}

if ( divtp_by_sfbm(vl,lcg,fmul,&q) ) {

pp_sfp(vl,fmul,fp);

cont_pp_sfp(vl,fmul,&cont,fp);

pp_sfp(vl,q,cofp);

cont_pp_sfp(vl,q,&cont,cofp);

return 1;

} else

return 0;

}

void const_term(f,c)

void const_term(P f,UM c)

P f;

UM c;

{

DCP dc;

Line 1233 UM c;

Line 1256 UM c;

DEG(c) = -1;

}

void const_term_sfbm(f,bound,c)

void const_term_sfbm(BM f,UM c)

BM f;

int bound;

UM c;

{

int i;

int i,dy;

for ( i = 0; i < bound; i++ )

dy = DEG(f);

for ( i = 0; i <= dy; i++ )

if ( DEG(COEF(f)[i]) >= 0 )

COEF(c)[i] = COEF(COEF(f)[i])[0];

else

COEF(c)[i] = 0;

degum(c,bound-1);

degum(c,dy);

}

/* lcy*(product of const part) | lcg0 ? */

int sfctest(lcg0,lcy,list,k,in)

int sfctest(UM lcg0,BM lcy,ML list,int k,int *in)

UM lcg0;

BM lcy;

ML list;

int k;

int *in;

{

DCP dc;

int dy,i,dr;

int bound,i,dr;

UM t,s,u,w;

BM *l;

bound = list->bound;

dy = list->bound;

t = W_UMALLOC(2*bound);

t = W_UMALLOC(2*dy);

s = W_UMALLOC(2*bound);

s = W_UMALLOC(2*dy);

u = W_UMALLOC(2*bound);

u = W_UMALLOC(2*dy);

const_term_sfbm(lcy,bound,t);

const_term_sfbm(lcy,t);

if ( DEG(t) < 0 )

return 1;

l = (BM *)list->c;

for ( i = 0; i < k; i++ ) {

const_term_sfbm(l[in[i]],bound,s);

const_term_sfbm(l[in[i]],s);

mulsfum(t,s,u);

if ( DEG(u) >= bound )

if ( DEG(u) > dy )

degum(u,bound-1);

degum(u,dy);

w = t; t = u; u = w;

}

cpyum(lcg0,s);

Line 1288 int *in;

Line 1303 int *in;

/* main var of f is x */

void mulsfbmarray(n,lcy,list,k,in,x,y,g)

void mulsfbmarray(int dx,BM lcy,ML list,int k,int *in,V x,V y,P *g)

int n;

BM lcy;

ML list;

int k;

int *in;

V x,y;

P *g;

{

int bound,i;

int dy,i;

BM wb0,wb1,t,lcbm;

BM wb0,wb1,t;

BM *l;

bound = list->bound;

dy = list->bound;

W_BMALLOC(n,bound,wb0); W_BMALLOC(n,bound,wb1);

W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1);

l = (BM *)list->c;

clearsfbm(bound,n,wb0);

clearbm(dx,wb0);

mulsfbm(bound,lcy,l[in[0]],wb0);

mulsfbm(lcy,l[in[0]],wb0);

for ( i = 1; i < k; i++ ) {

clearsfbm(bound,n,wb1);

clearbm(dx,wb1);

mulsfbm(bound,l[in[i]],wb0,wb1);

mulsfbm(l[in[i]],wb0,wb1);

t = wb0; wb0 = wb1; wb1 = t;

}

sfbmtop(bound,wb0,x,y,g);

sfbmtop(wb0,x,y,g);

}

void sfcsump(vl,f,s)

void sfcsump(VL vl,P f,P *s)

VL vl;

P f;

P *s;

{

P t,u;

DCP dc;

Line 1330 P *s;

Line 1335 P *s;

/* *fp = primitive part of f w.r.t. x */

void pp_sfp(vl,f,fp)

void cont_pp_sfp(VL vl,P f,P *cp,P *fp)

VL vl;

P f;

P *fp;

{

V x,y;

int d;

UM t,s,gcd;

DCP dc;

P dvr;

GFS g;

x = vl->v;

y = vl->next->v;

d = getdeg(y,f);

if ( d == 0 )

if ( d == 0 ) {

itogfs(1,&g);

*cp = (P)g;

*fp = f; /* XXX */

else {

} else {

t = W_UMALLOC(2*d);

s = W_UMALLOC(2*d);

gcd = W_UMALLOC(2*d);

Line 1357 P *fp;

Line 1361 P *fp;

gcdsfum(gcd,t,s);

cpyum(s,gcd);

}

sfumtop(y,gcd,&dvr);

sfumtop(y,gcd,cp);

divsp(vl,f,dvr,fp);

divsp(vl,f,*cp,fp);

}

int divtp_by_sfbm(vl,f,g,qp)

int divtp_by_sfbm(VL vl,P f,P g,P *qp)

VL vl;

P f,g;

P *qp;

{

V x,y;

int fx,fy,gx,gy;

Line 1380 P *qp;

Line 1381 P *qp;

if ( fx < gx || fy < gy )

return 0;

W_BMALLOC(fx,fy+1,fl); ptosfbm(fy+1,f,fl); cf = COEF(fl);

W_BMALLOC(fx,fy,fl); ptosfbm(fy,f,fl); cf = COEF(fl);

W_BMALLOC(gx,gy+1,gl); ptosfbm(gy+1,g,gl); cg = COEF(gl);

W_BMALLOC(gx,gy,gl); ptosfbm(gy,g,gl); cg = COEF(gl);

W_BMALLOC(fx-gx,fy-gy+1,ql); cq = COEF(ql);

W_BMALLOC(fx-gx,fy-gy,ql); cq = COEF(ql);

hg = cg[gy];

q = W_UMALLOC(fx); t = W_UMALLOC(fx); s = W_UMALLOC(fx);

Line 1405 P *qp;

Line 1406 P *qp;

for ( j = gy-1; j >= 0 && DEG(cf[j]) < 0; j-- );

if ( j >= 0 )

return 0;

sfbmtop(DEG(ql),ql,x,y,qp);

sfbmtop(ql,x,y,qp);

return 1;

}

/* XXX generate an irreducible poly of degree n */

extern int current_gfs_q1;

extern int *current_gfs_ntoi;

void generate_defpoly_sfum(int n,UM *dp)

{

UM r,dr,t,g;

UM *f;

int *c,*w;

int max,i,j;

*dp = r = UMALLOC(n);

DEG(r) = n;

c = COEF(r);

c[n] = _onesf();

max = current_gfs_q1;

w = (int *)ALLOCA(n*sizeof(int));

bzero(w,n*sizeof(int));

dr = W_UMALLOC(n); t = W_UMALLOC(n); g = W_UMALLOC(n);

f = (UM *)ALLOCA((n+1)*sizeof(UM));

while ( 1 ) {

for ( i = 0; i < n && w[i] == max; i++ );

if ( i == n ) {

/* XXX cannot happen */

error("generate_defpoly_sfum : cannot happen");

}

for ( j = 0; j < i; j++ )

w[j] = 0;

w[i]++;

if ( !current_gfs_ntoi )

for ( i = 0; i < n; i++ )

c[i] = w[i]?FTOIF(w[i]):0;

else

for ( i = 0; i < n; i++ )

c[i] = w[i]?FTOIF(w[i]-1):0;

if ( !c[0] )

continue;

diffsfum(r,dr); cpyum(r,t); gcdsfum(t,dr,g);

if ( DEG(g) > 0 )

continue;

czsfum(r,f);

for ( i = 0; f[i]; i++ );

if ( i == 1 )

return;

}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>