OpenXM_contrib/pari-2.2/src/basemath/trans2.c - diff

Return to trans2.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/trans2.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:06

version 1.2, 2002/09/11 07:26:53

Line 15 Foundation, Inc., 59 Temple Place - Suite 330, Boston,

/********************************************************************/

/** **/

/** TRANSCENDENTAL FONCTIONS **/

/** TRANSCENDENTAL FUNCTIONS **/

/** (part 2) **/

/** **/

/********************************************************************/

Line 33 static GEN mpach(GEN x);

static GEN

mpatan(GEN x)

{

long l,l1,l2,n,m,u,i,av0,av,lp,e,sx,s;

long l, l1, l2, n, m, u, i, lp, e, sx, s;

gpmem_t av0, av;

double alpha,beta,gama=1.0,delta,fi;

GEN y,p1,p2,p3,p4,p5,unr;

sx=signe(x);

if (!sx)

if (!sx) return realzero_bit(expo(x));

{

y=cgetr(3); y[1]=x[1]; y[2]=0; return y;

}

l = lp = lg(x);

if (sx<0) setsigne(x,1);

u = cmprs(x,1);

Line 126 mpatan(GEN x)

Line 124 mpatan(GEN x)

GEN

gatan(GEN x, long prec)

{

long av,tetpil;

gpmem_t av, tetpil;

GEN y,p1;

switch(typ(x))

Line 161 gatan(GEN x, long prec)

Line 159 gatan(GEN x, long prec)

void

gatanz(GEN x, GEN y)

{

long av=avma,prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gatanz");

gaffect(gatan(x,prec),y); avma=av;

Line 177 gatanz(GEN x, GEN y)

Line 176 gatanz(GEN x, GEN y)

static GEN

mpasin(GEN x)

{

long l,u,v,av;

long l, u, v;

gpmem_t av;

GEN y,p1;

u=cmprs(x,1); v=cmpsr(-1,x);

Line 199 mpasin(GEN x)

GEN

gasin(GEN x, long prec)

{

long av,tetpil,l,sx;

long l, sx;

gpmem_t av, tetpil;

GEN y,p1;

switch(typ(x))

{

case t_REAL: sx=signe(x);

if (!sx) { y=cgetr(3); y[1]=x[1]; y[2]=0; return y; }

if (!sx) return realzero_bit(expo(x));

if (sx<0) setsigne(x,1);

if (cmpsr(1,x)>=0) { setsigne(x,sx); return mpasin(x); }

Line 248 gasin(GEN x, long prec)

Line 249 gasin(GEN x, long prec)

void

gasinz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gasinz");

gaffect(gasin(x,prec),y); avma=av;

Line 264 gasinz(GEN x, GEN y)

Line 266 gasinz(GEN x, GEN y)

static GEN

mpacos(GEN x)

{

long l,u,v,sx,av;

long l, u, v, sx;

gpmem_t av;

GEN y,p1,p2;

u=cmprs(x,1); v=cmpsr(-1,x); sx = signe(x);

Line 274 mpacos(GEN x)

Line 277 mpacos(GEN x)

y=mppi(2-l); setexpo(y,0); return y;

}

l=lg(x);

if (!u)

if (!u) return realzero_bit( -(bit_accuracy(l)>>1) );

{

y = cgetr(3);

y[1] = evalexpo(-(bit_accuracy(l)>>1));

y[2] = 0; return y;

}

if (!v) return mppi(l);

y=cgetr(l); av=avma;

Line 309 mpacos(GEN x)

Line 307 mpacos(GEN x)

GEN

gacos(GEN x, long prec)

{

long av,tetpil,l,sx;

long l, sx;

gpmem_t av, tetpil;

GEN y,p1;

switch(typ(x))

Line 355 gacos(GEN x, long prec)

Line 354 gacos(GEN x, long prec)

void

gacosz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gacosz");

gaffect(gacos(x,prec),y); avma=av;

Line 377 mparg(GEN x, GEN y)

sx=signe(x); sy=signe(y);

if (!sy)

{

if (sx>0)

if (sx>0) return realzero_bit(expo(y) - expo(x));

{

theta=cgetr(3); theta[1]=y[1]-expo(x);

theta[2]=0; return theta;

}

return mppi(lg(x));

}

prec = lg(y); if (prec<lg(x)) prec = lg(x);

Line 423 rfix(GEN x,long prec)

Line 419 rfix(GEN x,long prec)

static GEN

sarg(GEN x, GEN y, long prec)

{

long av=avma;

gpmem_t av=avma;

x = rfix(x,prec); y = rfix(y,prec);

return gerepileupto(av,mparg(x,y));

}

Line 432 GEN

Line 428 GEN

garg(GEN x, long prec)

{

GEN p1;

long av,tx=typ(x),tetpil;

long tx=typ(x);

gpmem_t av, tetpil;

if (gcmp0(x)) err(talker,"zero argument in garg");

switch(tx)

Line 458 garg(GEN x, long prec)

Line 455 garg(GEN x, long prec)

/********************************************************************/

/** **/

/** FONCTION COSINUS HYPERBOLIQUE **/

/** HYPERBOLIC COSINE **/

/** **/

/********************************************************************/

static GEN

mpch(GEN x)

{

long l,av;

gpmem_t av;

GEN y,p1;

if (gcmp0(x)) return gaddsg(1,x);

l=lg(x); y=cgetr(l); av=avma;

y = cgetr(lg(x)); av = avma;

p1=mpexp(x); p1 = addrr(p1, divsr(1,p1));

p1 = mpexp(x); p1 = addrr(p1, ginv(p1));

setexpo(p1, expo(p1)-1);

affrr(p1,y); avma=av; return y;

affrr(p1, y); avma = av; return y;

}

GEN

gch(GEN x, long prec)

{

long av,tetpil;

gpmem_t av;

GEN p1;

switch(typ(x))

Line 487 gch(GEN x, long prec)

Line 484 gch(GEN x, long prec)

case t_REAL:

return mpch(x);

case t_SER:

if (gcmp0(x) && valp(x) == 0) return gcopy(x);

case t_COMPLEX:

av=avma; p1=gexp(x,prec);

av = avma; p1 = gexp(x,prec); p1 = gadd(p1, ginv(p1));

p1=gadd(p1,ginv(p1)); tetpil=avma;

return gerepileupto(av, gmul2n(p1,-1));

return gerepile(av,tetpil,gmul2n(p1,-1));

case t_SER:

av=avma; p1=gexp(x,prec); p1=gadd(p1,gdivsg(1,p1));

tetpil=avma; return gerepile(av,tetpil,gmul2n(p1,-1));

case t_INTMOD: case t_PADIC:

err(typeer,"gch");

}

Line 505 gch(GEN x, long prec)

Line 499 gch(GEN x, long prec)

void

gchz(GEN x, GEN y)

{

long av=avma,prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gchz");

gaffect(gch(x,prec),y); avma=av;

Line 513 gchz(GEN x, GEN y)

Line 508 gchz(GEN x, GEN y)

/********************************************************************/

/** **/

/** FONCTION SINUS HYPERBOLIQUE **/

/** HYPERBOLIC SINE **/

/** **/

/********************************************************************/

static GEN

mpsh(GEN x)

{

long l,av;

gpmem_t av;

GEN y,p1;

if (!signe(x))

if (!signe(x)) return realzero_bit(expo(x));

{

y = cgetr(lg(x)); av = avma;

y=cgetr(3); y[1]=x[1]; y[2]=0; return y;

p1 = mpexp(x); p1 = addrr(p1, divsr(-1,p1));

}

l=lg(x); y=cgetr(l); av=avma;

p1=mpexp(x); p1 = addrr(p1, divsr(-1,p1));

setexpo(p1, expo(p1)-1);

affrr(p1,y); avma=av; return y;

affrr(p1,y); avma = av; return y;

}

GEN

gsh(GEN x, long prec)

{

long av,tetpil;

gpmem_t av;

GEN p1;

switch(typ(x))

Line 544 gsh(GEN x, long prec)

Line 536 gsh(GEN x, long prec)

case t_REAL:

return mpsh(x);

case t_SER:

if (gcmp0(x) && valp(x) == 0) return gcopy(x);

case t_COMPLEX:

av=avma; p1=gexp(x,prec);

av = avma; p1 = gexp(x,prec); p1 = gsub(p1, ginv(p1));

p1=gsub(p1,ginv(p1)); tetpil=avma;

return gerepileupto(av, gmul2n(p1,-1));

return gerepile(av,tetpil,gmul2n(p1,-1));

case t_SER:

if (gcmp0(x)) return gcopy(x);

av=avma; p1=gexp(x,prec); p1=gsub(p1,gdivsg(1,p1));

tetpil=avma; return gerepile(av,tetpil,gmul2n(p1,-1));

case t_INTMOD: case t_PADIC:

err(typeer,"gsh");

}

Line 564 gsh(GEN x, long prec)

Line 551 gsh(GEN x, long prec)

void

gshz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gshz");

gaffect(gsh(x,prec),y); avma=av;

Line 579 gshz(GEN x, GEN y)

Line 567 gshz(GEN x, GEN y)

static GEN

mpth(GEN x)

{

long l,av;

long l;

gpmem_t av;

GEN y,p1,p2;

if (!signe(x))

if (!signe(x)) return realzero_bit(expo(x));

{

y=cgetr(3); y[1]=x[1]; y[2]=0;

return y;

}

l=lg(x); y=cgetr(l); av=avma;

p1=cgetr(l+1); affrr(x,p1);

Line 600 mpth(GEN x)

Line 585 mpth(GEN x)

GEN

gth(GEN x, long prec)

{

long av,tetpil;

gpmem_t av, tetpil;

GEN p1,p2,p3;

switch(typ(x))

Line 629 gth(GEN x, long prec)

Line 614 gth(GEN x, long prec)

void

gthz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gthz");

gaffect(gth(x,prec),y); avma=av;

Line 645 gthz(GEN x, GEN y)

Line 631 gthz(GEN x, GEN y)

static GEN

mpash(GEN x)

{

long s=signe(x),av;

long s=signe(x);

gpmem_t av;

GEN y,p1;

y=cgetr(lg(x)); av=avma;

Line 658 mpash(GEN x)

Line 645 mpash(GEN x)

GEN

gash(GEN x, long prec)

{

long av,tetpil,sx,sy,sz;

long sx, sy, sz;

gpmem_t av, tetpil;

GEN y,p1;

if (gcmp0(x)) return gcopy(x);

Line 699 gash(GEN x, long prec)

Line 687 gash(GEN x, long prec)

void

gashz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gashz");

gaffect(gash(x,prec),y); avma=av;

Line 714 gashz(GEN x, GEN y)

Line 703 gashz(GEN x, GEN y)

static GEN

mpach(GEN x)

{

long l,av;

long l;

gpmem_t av;

GEN y,p1;

l=lg(x); y=cgetr(l); av=avma;

Line 730 mpach(GEN x)

Line 720 mpach(GEN x)

GEN

gach(GEN x, long prec)

{

long av,tetpil;

gpmem_t av;

GEN y,p1;

long v;

switch(typ(x))

{

case t_REAL:

if (gcmpgs(x,1)>=0) return mpach(x);

if (gcmpgs(x,1) >= 0) return mpach(x);

y=cgetg(3,t_COMPLEX);

y = cgetg(3,t_COMPLEX);

if (gcmpgs(x,-1)>=0)

if (gcmpgs(x,-1) >= 0)

{

y[2]=lmpacos(x); y[1]=zero;

y[1] = zero;

return y;

y[2] = lmpacos(x); return y;

}

av=avma; p1=mpach(gneg_i(x)); tetpil=avma;

av = avma;

y[1]=lpile(av,tetpil,gneg(p1));

y[1] = lpileupto(av, gneg(mpach(gneg_i(x))));

y[2]=lmppi(lg(x));

y[2] = lmppi(lg(x)); return y;

return y;

case t_COMPLEX:

av=avma; p1=gaddsg(-1,gsqr(x));

av = avma; p1 = gaddsg(-1,gsqr(x));

p1=gadd(x,gsqrt(p1,prec)); tetpil=avma;

p1 = gadd(x, gsqrt(p1,prec)); /* x + sqrt(x^2-1) */

y=glog(p1,prec);

y = glog(p1,prec);

if (signe(y[2])<0) { tetpil=avma; y=gneg(y); }

if (signe(y[2]) < 0) y = gneg(y);

return gerepile(av,tetpil,y);

return gerepileupto(av, y);

case t_SER:

av=avma; if (valp(x)<0) err(negexper,"gach");

av = avma; v = valp(x);

p1=gdiv(derivser(x),gsqrt(gsubgs(gsqr(x),1),prec));

if (v < 0) err(negexper,"gach");

y=integ(p1,varn(x));

if (gcmp0(x))

if (!valp(x) && gcmp1((GEN)x[2])) return gerepileupto(av,y);

if (valp(x))

{

p1=cgetg(3,t_COMPLEX); p1[1]=zero; p1[2]=lmppi(prec);

if (!v) return gcopy(x);

setexpo(p1[2],0);

return gerepileupto(av, gadd(x, PiI2n(prec,-1)));

}

else p1=gach((GEN)x[2],prec);

p1 = gdiv(derivser(x), gsqrt(gsubgs(gsqr(x),1),prec));

tetpil=avma; return gerepile(av,tetpil,gadd(p1,y));

y = integ(p1, varn(x));

if (v)

p1 = PiI2n(prec, -1); /* I Pi/2 */

else

{

p1 = (GEN)x[2];

if (gcmp1(p1)) return gerepileupto(av,y);

p1 = gach(p1, prec);

}

return gerepileupto(av, gadd(p1,y));

case t_INTMOD: case t_PADIC:

err(typeer,"gach");

Line 778 gach(GEN x, long prec)

Line 775 gach(GEN x, long prec)

void

gachz(GEN x, GEN y)

{

long av=avma,prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gachz");

gaffect(gach(x,prec),y); avma=av;

Line 790 gachz(GEN x, GEN y)

Line 788 gachz(GEN x, GEN y)

/** **/

/********************************************************************/

/* |x| < 1 */

static GEN

mpath(GEN x)

{

long av;

gpmem_t av;

GEN y,p1;

if (!signe(x))

if (!signe(x)) return realzero_bit(expo(x));

{

y=cgetr(3); y[1]=x[1]; y[2]=0;

return y;

}

y=cgetr(lg(x)); av=avma;

p1 = addrs(divsr(2,subsr(1,x)),-1);

affrr(mplog(p1), y); avma=av;

Line 810 mpath(GEN x)

Line 805 mpath(GEN x)

GEN

gath(GEN x, long prec)

{

long av,tetpil;

gpmem_t av, tetpil;

GEN y,p1;

switch(typ(x))

Line 847 gath(GEN x, long prec)

Line 842 gath(GEN x, long prec)

void

gathz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gathz");

gaffect(gath(x,prec),y); avma=av;

Line 863 gathz(GEN x, GEN y)

Line 859 gathz(GEN x, GEN y)

void

mpbern(long nb, long prec)

{

long n,m,i,j,d,d1,d2,l,av,av2,code0;

long n, m, i, j, d, d1, d2, l, code0;

gpmem_t av, av2;

GEN p1,p2, B;

if (nb < 0) nb = 0;

Line 888 mpbern(long nb, long prec)

Line 885 mpbern(long nb, long prec)

{

fprintferr("caching Bernoulli numbers 2*%ld to 2*%ld, prec = %ld\n",

i,nb,prec);

timer2();

(void)timer2();

}

p2 = p1; av2=avma;

Line 926 bernreal(long n, long prec)

Line 923 bernreal(long n, long prec)

{

GEN B;

if (n==1) { B=cgetr(prec); affsr(-1,B); setexpo(B,-1); return B; }

if (n==1) { B = stor(-1, prec); setexpo(B,-1); return B; }

if (n<0 || n&1) return gzero;

n >>= 1; mpbern(n+1,prec); B=cgetr(prec);

affrr(bern(n),B); return B;

Line 936 bernreal(long n, long prec)

Line 933 bernreal(long n, long prec)

static GEN

bernfracspec(long k)

{

long n,av,lim, K = k+1;

ulong n, K = k+1;

GEN s,c,N,h;

gpmem_t av, lim;

GEN s, c, N, b;

c = N = stoi(K); s = gun; h = gzero;

c = N = utoi(K); s = gun; b = gzero;

av = avma; lim = stack_lim(av,2);

for (n=2; ; n++)

for (n=2; ; n++) /* n <= k+1 */

{

c = gdivgs(gmulgs(c,n-k-2),n);

c = diviiexact(muliu(c,k+2-n), utoi(n));

h = gadd(h, gdivgs(gmul(c,s), n));

if (n & 1) setsigne(c, 1); else setsigne(c, -1);

if (n==K) return gerepileupto(av,h);

/* c = (-1)^(n-1) binomial(k+1, n), s = 1^k + ... + (n-1)^k */

N[2] = n; s = addii(s, gpuigs(N,k));

b = gadd(b, gdivgs(mulii(c,s), n));

if (n == K) return gerepileupto(av, b);

N[2] = n; s = addii(s, gpowgs(N,k));

if (low_stack(lim, stack_lim(av,2)))

{

GEN *gptr[3]; gptr[0]=&c; gptr[1]=&h; gptr[2]=&s;

if (DEBUGMEM>1) err(warnmem,"bernfrac");

gerepilemany(av,gptr,3);

gerepileall(av,3, &c,&b,&s);

}

Line 960 GEN

Line 961 GEN

bernfrac(long k)

{

if (!k) return gun;

if (k==1) return gneg(ghalf);

if (k == 1) return gneg(ghalf);

if (k<0 || k&1) return gzero;

if (k < 0 || k & 1) return gzero;

return bernfracspec(k);

}

static GEN

bernvec2(long k)

{

GEN B = cgetg(k+2,t_VEC);

long i;

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

B[i+1]=(long)bernfracspec(i<<1);

B[1]=un; return B;

}

/* mpbern as exact fractions */

GEN

bernvec(long nb)

{

long n,m,i,j,d1,d2,av,tetpil;

GEN y = cgetg(nb+2,t_VEC);

GEN p1,y;

long n, i;

if (nb < 300) return bernvec2(nb);

if (nb > 46340 && BITS_IN_LONG == 32) err(impl, "bernvec for n > 46340");

y=cgetg(nb+2,t_VEC); y[1]=un;

y[1] = un;

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

for (n = 1; n <= nb; n++)

{

{ /* compute y[n+1] = B_{2n} */

av=avma; p1=gzero;

gpmem_t av = avma;

n=8; m=5; d1=i-1; d2=2*i-3;

GEN b = gmul2n(stoi(1-2*n), -1); /* 1 + (2n+1)B_1 = (1-2n) /2 */

for (j=d1; j>0; j--)

GEN c = gun;

{

ulong u1 = 2*n + 1, u2 = n, d1 = 1, d2 = 1;

p1 = gmulsg(n*m, gadd(p1,(GEN)y[j+1]));

p1 = gdivgs(p1, d1*d2);

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

n+=4; m+=2; d1--; d2-=2;

{

c = diviiexact(muliu(c, u1*u2), utoi(d1*d2)); /* = binomial(2n+1, 2*i) */

b = gadd(b, gmul(c, (GEN)y[i+1]));

u1 -= 2; u2--; d1++; d2 += 2;

}

p1 = gsubsg(1,gdivgs(gaddsg(1,p1),2*i+1));

y[n+1] = lpileupto(av, gdivgs(b, -(1+2*n)));

tetpil=avma; p1 = gmul2n(p1,-2*i);

y[i+1] = lpile(av,tetpil,p1);

}

return y;

}

Line 1013 static GEN

Line 1004 static GEN

mpgamma(GEN x)

{

GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;

long l,l1,l2,u,i,k,e,s,sx,n,p,av,av1;

long l, l1, l2, u, i, k, e, s, sx, n, p;

gpmem_t av, av1;

double alpha,beta,dk;

sx=signe(x); if (!sx) err(gamer2);

Line 1087 static GEN

Line 1079 static GEN

cxgamma(GEN x, long prec)

{

GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;

long l,l1,l2,u,i,k,e,s,n,p,av,av1;

long l, l1, l2, u, i, k, e, s, n, p;

gpmem_t av, av1;

double alpha,beta,dk;

if (gcmp0((GEN)x[2])) return ggamma((GEN)x[1],prec);

Line 1200 double

Line 1193 double

dnorm(double s, double t) { return s*s + t*t; }

GEN

trans_fix_arg(long *prec, GEN *s0, GEN *sig, long *av, GEN *res)

trans_fix_arg(long *prec, GEN *s0, GEN *sig, gpmem_t *av, GEN *res)

{

GEN s, p1;

long l;

Line 1232 GEN

Line 1225 GEN

gammanew(GEN s0, long la, long prec)

{

GEN s, u, a, y, res, tes, sig, invn2, p1, unr, nnx, pitemp;

long i, nn, lim, av, av2, avlim;

long i, lim, nn;

gpmem_t av, av2, avlim;

int funeq = 0;

if (DEBUGLEVEL) timer2();

if (DEBUGLEVEL) (void)timer2();

s = trans_fix_arg(&prec,&s0,&sig,&av,&res);

if (signe(sig) <= 0 || expo(sig) < -1)

Line 1278 gammanew(GEN s0, long la, long prec)

Line 1272 gammanew(GEN s0, long la, long prec)

nn = 1;

if (DEBUGLEVEL) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn);

if ((ulong)nn >= maxprime()) err(primer1);

}

prec++; unr = realun(prec);

Line 1338 gammanew(GEN s0, long la, long prec)

Line 1331 gammanew(GEN s0, long la, long prec)

GEN

ggamma(GEN x, long prec)

{

gpmem_t av=avma, tetpil;

long m, ma;

GEN p1,p2,p3;

switch(typ(x))

{

case t_INT:

if (signe(x)<=0) err(gamer2);

return transc(ggamma,x,prec);

if (cmpis(x,481177) > 0) err(talker,"argument too large in ggamma");

/* heuristic */

if (cmpis(x,350 + 70*(prec-2)) > 0) return transc(ggamma,x,prec);

p2 = cgetr(prec); av = avma;

p1 = mpfact(itos(x) - 1);

affir(p1,p2); avma = av;

return p2;

case t_REAL:

return mpgamma(x);

case t_FRAC:

if (cmpis((GEN)x[2],2) == 0)

{

p2 = cgetr(prec); av = avma;

if (cmpis(mpabs((GEN)x[1]),962354) > 0)

err(talker,"argument too large in ggamma");

/* heuristic */

if (cmpis((GEN)x[1],200 + 50*(prec-2)) > 0)

return transc(ggamma,x,prec);

p1 = gsqrt(mppi(prec),prec);

m = itos((GEN)x[1]) - 1; ma = labs(m);

p3 = gmul2n(divii(mpfact(ma),mpfact(ma/2)),-ma);

if (m >= 0) affrr(gmul(p3,p1),p2);

else

{

affrr(gdiv(p1,p3),p2);

if ((m&3) == 2) setsigne(p2,-1);

}

avma = av;

return p2;

}

else return transc(ggamma,x,prec);

case t_FRACN:

p1 = gred(x);

tetpil = avma;

return gerepile(av,tetpil,ggamma(p1,prec));

case t_COMPLEX:

return cxgamma(x,prec);

Line 1364 ggamma(GEN x, long prec)

Line 1395 ggamma(GEN x, long prec)

void

ggammaz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"ggammaz");

gaffect(ggamma(x,prec),y); avma=av;

Line 1374 static GEN

Line 1406 static GEN

mplngamma(GEN x)

{

GEN z,y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;

long l,l1,l2,u,i,k,e,f,s,sx,n,p,av,av1;

long l, l1, l2, u, i, k, e, f, s, sx, n, p;

gpmem_t av, av1;

double alpha,beta,dk;

sx=signe(x); if (!sx) err(talker,"zero argument in mplngamma");

Line 1466 mplngamma(GEN x)

Line 1499 mplngamma(GEN x)

}

/* t_REAL result */

y[3] = y[0]; y += 3;

affrr(p4,y); avma = (long)y; return y;

affrr(p4,y); avma = (gpmem_t)y; return y;

}

static GEN

cxlngamma(GEN x, long prec)

{

GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;

long l,l1,l2,flag,i,k,e,s,n,p,av,av1;

long l, l1, l2, flag, i, k, e, s, n, p;

gpmem_t av, av1;

double alpha,beta,dk;

if (gcmp0((GEN)x[2])) return glngamma((GEN)x[1],prec);

Line 1581 cxlngamma(GEN x, long prec)

Line 1615 cxlngamma(GEN x, long prec)

GEN

glngamma(GEN x, long prec)

{

long i,av,n;

long i, n;

GEN y,p1;

gpmem_t av;

GEN y,p1,p2;

switch(typ(x))

{

case t_INT:

if (signe(x)<=0) err(gamer2);

return transc(glngamma,x,prec);

p2 = cgetr(prec); av = avma;

/* heuristic */

if (cmpis(x,200 + 50*(prec-2)) > 0)

return transc(glngamma,x,prec);

p1 = glog(mpfact(itos(x) - 1),prec);

affrr(p1,p2); avma = av;

return p2;

case t_REAL:

return mplngamma(x);

Line 1618 glngamma(GEN x, long prec)

Line 1659 glngamma(GEN x, long prec)

void

glngammaz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"glngammaz");

gaffect(glngamma(x,prec),y); avma=av;

Line 1633 glngammaz(GEN x, GEN y)

Line 1675 glngammaz(GEN x, GEN y)

static GEN

mpgamd(long x, long prec)

{

long i,j,a,l,av;

long i, j, a, l;

gpmem_t av;

GEN y,p1,p2;

a = labs(x);

Line 1660 mpgamd(long x, long prec)

Line 1703 mpgamd(long x, long prec)

void

mpgamdz(long s, GEN y)

{

long av=avma;

gpmem_t av=avma;

affrr(mpgamd(s,lg(y)),y); avma=av;

}

Line 1668 mpgamdz(long s, GEN y)

Line 1711 mpgamdz(long s, GEN y)

GEN

ggamd(GEN x, long prec)

{

long av,tetpil;

gpmem_t av, tetpil;

switch(typ(x))

{

Line 1690 ggamd(GEN x, long prec)

Line 1733 ggamd(GEN x, long prec)

void

ggamdz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"ggamdz");

gaffect(ggamd(x,prec),y); avma=av;

Line 1702 ggamdz(GEN x, GEN y)

Line 1746 ggamdz(GEN x, GEN y)

/** **/

/********************************************************************/

#if 1

GEN

psinew(GEN s0, long prec)

{

gpmem_t av;

GEN sum,z,a,res,tes,in2,sig,s,unr;

long lim,nn,k;

const long la = 3;

if (DEBUGLEVEL>2) (void)timer2();

s = trans_fix_arg(&prec,&s0,&sig,&av,&res);

if (signe(sig) <= 0)

{

GEN pi = mppi(prec);

z = gsub(psinew(gsubsg(1,s), prec), gmul(pi, gcotan(gmul(pi,s), prec)));

gaffect(z, res); avma = av; return res;

}

{

double ssig = rtodbl(sig);

double st = rtodbl(gimag(s));

double l;

{

double rlog, ilog; /* log (s-gamma) */

dcxlog(ssig - rtodbl(mpeuler(3)), st, &rlog,&ilog);

l = dnorm(rlog,ilog);

}

if (l < 0.000001) l = 0.000001;

l = log(l) / 2.;

lim = 2 + (long)ceil((pariC2*(prec-2) - l) / (2*(1+log((double)la))));

if (lim < 2) lim = 2;

l = (2*lim-1)*la / (2.*PI);

l = l*l - st*st;

if (l < 0.) l = 0.;

nn = (long)ceil( sqrt(l) - ssig );

if (nn < 1) nn = 1;

if (DEBUGLEVEL>2) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn);

}

prec++; unr = realun(prec); /* one extra word of precision */

a = gdiv(unr, gaddgs(s, nn)); /* 1 / (s+n) */

sum = gmul2n(a,-1);

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

sum = gadd(sum, gdiv(unr, gaddgs(s, k)));

z = gsub(glog(gaddgs(s, nn), prec), sum);

if (DEBUGLEVEL>2) msgtimer("sum from 0 to N-1");

tes = divrs(bernreal(2*lim, prec), 2*lim);

in2 = gsqr(a);

for (k=2*lim-2; k>=2; k-=2)

{

tes = gadd(gmul(in2,tes), divrs(bernreal(k, prec), k));

}

if (DEBUGLEVEL>2) msgtimer("Bernoulli sum");

z = gsub(z, gmul(in2,tes));

gaffect(z, res); avma = av; return res;

}

#if 0

static GEN

mppsi(GEN z) /* version p=2 */

{

long l,n,k,x,xx,av,av1,tetpil;

long l, n, k, x, xx;

gpmem_t av, av1, tetpil;

GEN zk,u,v,a,b;

av=avma; l=lg(z);

x=(long)(1 + (bit_accuracy(l)>>1)*LOG2 + 1.58*rtodbl(absr(z)));

if (expo(z)>=15 || x>46340) err(impl,"psi(x) for x>=29000");

xx=x*x; n=(long)(1+3.591*x);

affsr(x,a=cgetr(l));

a = stor(x, l);

a = mplog(a);

gaffect(a,u=cgetr(l));

gaffsg(1,b=cgetr(l));

Line 1727 mppsi(GEN z) /* version p=2 */

Line 1831 mppsi(GEN z) /* version p=2 */

}

tetpil=avma; return gerepile(av,tetpil,divrr(u,v));

}

#else

static GEN /* by Manfred Radimersky */

mppsi(GEN z)

{

Line 1750 mppsi(GEN z)

Line 1853 mppsi(GEN z)

gerepile(head, tail, subrr(s, a));

}

a = cgetr(len);

a = stor(0, len);

affsr(0, a);

x = cgetr(len);

affrr(z, x);

tail = avma;

while(cmprs(x, num >> 2) < 0) {

gaddz(a, divsr(1, x), a);

gaddz(a, ginv(x), a);

gaddgsz(x, 1, x);

avma = tail;

}

Line 1764 mppsi(GEN z)

Line 1866 mppsi(GEN z)

s = mplog(x);

gsubz(s, a, s);

b = gmul2n(x, 1);

gsubz(s, divsr(1, b), s);

gsubz(s, ginv(b), s);

mpbern(num, len);

affsr(-1, a);

gdivgsz(a, 2, a);

f = mulrr(x, x);

f = divsr(1, f);

f = ginv(f);

k = 1;

do {

gmulz(a, f, a);

Line 1783 mppsi(GEN z)

Line 1885 mppsi(GEN z)

return gerepilecopy(head, s);

}

#endif

#if 1

static GEN

cxpsi(GEN z, long prec) /* version p=2 */

{

long l,n,k,x,xx,av,av1,tetpil;

long l, n, k, x, xx;

gpmem_t av, av1, tetpil;

GEN zk,u,v,a,b,p1;

if (gcmp0((GEN)z[2])) return gpsi((GEN)z[1],prec);

Line 1810 cxpsi(GEN z, long prec) /* version p=2 */

Line 1910 cxpsi(GEN z, long prec) /* version p=2 */

}

tetpil=avma; return gerepile(av,tetpil,gdiv(u,v));

}

#else

GEN

cxpsi(GEN z, long prec) /* by Manfred Radimersky */

{

Line 1877 gpsi(GEN x, long prec)

Line 1976 gpsi(GEN x, long prec)

{

switch(typ(x))

{

case t_REAL:

case t_REAL: case t_COMPLEX:

return mppsi(x);

return psinew(x,prec);

case t_COMPLEX:

return cxpsi(x,prec);

case t_INTMOD: case t_PADIC:

err(typeer,"gpsi");

case t_SER:

Line 1894 gpsi(GEN x, long prec)

Line 1990 gpsi(GEN x, long prec)

void

gpsiz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gpsiz");

gaffect(gpsi(x,prec),y); avma=av;

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