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.26 and 1.28

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

version 1.28, 2002/11/22 07:32:10

Line 1

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

/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.27 2002/11/01 05:43:35 noro Exp $ */

#include "ca.h"

#include "inline.h"

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

/* Hensel related functions */

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

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

void sfgcdgen(P,ML,ML *);

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

void ptosfbm(int,P,BM);

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

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

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

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

{

int i;

int fn;

ML rlist;

BM fl;

VL vl,nvl;

V y;

int dx,dy,bound;

GFS ev;

P f1,t,c,sf;

Line 534 void sfhensel(int count,P f,V x,int degbound,GFS *evp,

reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&f1);

vl = nvl; f = f1;

}

y = vl->next->v;

if ( vl->next )

y = vl->next->v;

dx = getdeg(x,f);

dy = getdeg(y,f);

if ( dx == 1 ) {

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

return;

}

fn = sfberle(vl,f,count,&ev,&dc);

fn = sfberle(x,y,f,count,&ev,&dc);

if ( fn <= 1 ) {

/* fn == 0 => short of evaluation points */

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

Line 640 void sfhensel(int count,P f,V x,int degbound,GFS *evp,

Line 641 void sfhensel(int count,P f,V x,int degbound,GFS *evp,

/* main variable of f = x */

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

int sfberle(V x,V y,P f,int count,GFS *ev,DCP *dcp)

{

UM wf,wf1,wf2,wfs,gcd;

int fn,n;

GFS m,fm;

DCP dc,dct,dc0;

VL nvl;

VL vl;

V x,y;

P lc,lc0,f0;

Obj obj;

int j,q,index,i;

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

NEWVL(vl); vl->v = x;

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

NEWVL(NEXT(vl)); NEXT(vl)->v = y;

NEXT(NEXT(vl)) =0;

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);

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

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

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

}

/* a0*g+b0*h = 1 mod y -> a*g+b*h = 1 mod y^(dy+1) */

void sfexgcd_by_hensel(BM g,BM h,int dy,BM *ap,BM *bp)

{

int i,k;

int dx;

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

UM wc,wd,we,wz,wa1,wb1;

BM wz0,wz1;

int dg,dh;

BM a,b,c;

dg = degbm(g);

dh = degbm(h);

dx = dg+dh;

a = BMALLOC(dh,dy);

b = BMALLOC(dg,dy);

/* c holds a*g+b*h-1 */

c = BMALLOC(dg+dh,dy);

W_BMALLOC(dx,dy,wz0); W_BMALLOC(dx,dy,wz1);

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

wa1 = W_UMALLOC(2*dx); wb1 = W_UMALLOC(2*dx);

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

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(dg); cpyum(COEF(g)[0],w1);

w2 = W_UMALLOC(dh); cpyum(COEF(h)[0],w2);

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

eucsfum(w1,w2,wa,wb);

cpyum(wa,COEF(a)[0]); cpyum(wb,COEF(b)[0]);

/* initialize c to a*g+b*h-1 */

mulsfbm(a,g,c); mulsfbm(b,h,wz0); addtosfbm(wz0,c);

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

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

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

fprintf(stderr,".");

/* at this point, a*g+b*h = 1 mod y^k, c = a*g+b*h-1 */

/* wt = -((a*g+b*h-1)/y^k) */

cpyum(COEF(c)[k],wt);

for ( i = DEG(wt); i >= 0; i-- )

COEF(wt)[i] = _chsgnsf(COEF(wt)[i]);

/* compute wa1,wb1 s.t. wa1*g0+wb1*h0 = wt */

mulsfum(wa,wt,wa1); DEG(wa1) = divsfum(wa1,COEF(h)[0],q);

mulsfum(wa1,COEF(g)[0],wc); subsfum(wt,wc,wd);

DEG(wd) = divsfum(wd,COEF(h)[0],wb1);

/* c += ((wa1*g+wb1*h)*y^k mod y^(dy+1) */

/* wz0 = wa1*y^k */

clearbm(dx,wz0);

cpyum(wa1,COEF(wz0)[k]);

/* wz1 = wz0*g mod y^(dy+1) */

clearbm(dx,wz1);

mulsfbm(g,wz0,wz1);

/* c += wz1 */

addtosfbm(wz1,c);

/* wz0 = wb1*y^k */

clearbm(dx,wz0);

cpyum(wb1,COEF(wz0)[k]);

/* wz1 = wz0*h mod y^(dy+1) */

clearbm(dx,wz1);

mulsfbm(h,wz0,wz1);

/* c += wz1 */

addtosfbm(wz1,c);

/* a += wa1*y^k, b += wb1*y^k */

cpyum(wa1,COEF(a)[k]);

cpyum(wb1,COEF(b)[k]);

}

fprintf(stderr,"\n");

DEG(a) = dy;

DEG(b) = dy;

*ap = a;

*bp = b;

}

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

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

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

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

int dx,dy;

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

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

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

if ( list->n == 0 )

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

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

Line 1024 void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P

Line 1112 void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P