hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

cleanup some powering code

Browse source
This commit is contained in:
Fredrik Johansson 2013-09-19 00:30:04 +01:00
parent a51ed9cf30
commit e6d8bb2e16
5 changed files with 88 additions and 92 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

62
fmpcb/pow.c
View file

@ -129,49 +129,43 @@ _fmpcb_pow_fmprb_exp(fmpcb_t z, const fmpcb_t x, const fmprb_t y, long prec)
fmpcb_clear(t);
}
#define BINEXP_LIMIT 64
void
fmpcb_pow_fmprb(fmpcb_t z, const fmpcb_t x, const fmprb_t y, long prec)
{
if (fmprb_is_exact(y))
const fmpr_struct * ymid = fmprb_midref(y);
const fmpr_struct * yrad = fmprb_radref(y);
if (fmprb_is_zero(y))
{
if (fmpr_is_zero(fmprb_midref(y)))
{
fmpcb_one(z);
return;
}
fmpcb_one(z);
return;
}
if (!fmpr_is_special(fmprb_midref(y)))
if (fmpr_is_zero(yrad))
{
/* small half-integer or integer */
if (fmpr_cmpabs_2exp_si(ymid, BINEXP_LIMIT) < 0 &&
fmpr_is_int_2exp_si(ymid, -1))
{
const fmpz * exp_exp = fmpr_expref(fmprb_midref(y));
fmpz_t e;
fmpz_init(e);
/* smallish integer powers and square roots */
if (!COEFF_IS_MPZ(*exp_exp) && (*exp_exp >= -1L))
if (fmpr_is_int(ymid))
{
long exp_bits;
exp_bits = *exp_exp + fmpz_bits(fmpr_manref(fmprb_midref(y)));
if (exp_bits < 64)
{
fmpz_t e;
fmpz_init(e);
if (*exp_exp == -1L)
{
fmpz_set(e, fmpr_manref(fmprb_midref(y)));
fmpcb_sqrt(z, x, prec + exp_bits);
fmpcb_pow_fmpz(z, z, e, prec);
}
else
{
fmpz_mul_2exp(e, fmpr_manref(fmprb_midref(y)), *exp_exp);
fmpcb_pow_fmpz(z, x, e, prec);
}
fmpz_clear(e);
return;
}
fmpr_get_fmpz_fixed_si(e, ymid, 0);
fmpcb_pow_fmpz_binexp(z, x, e, prec);
}
else
{
fmpr_get_fmpz_fixed_si(e, ymid, -1);
fmpcb_sqrt(z, x, prec + fmpz_bits(e));
fmpcb_pow_fmpz_binexp(z, z, e, prec);
}
fmpz_clear(e);
return;
}
}

11
fmpcb/test/t-pow.c
View file

@ -36,7 +36,7 @@ int main()
flint_randinit(state);
/* check large arguments */
for (iter = 0; iter < 10000; iter++)
for (iter = 0; iter < 20000; iter++)
{
fmpcb_t a, b, c, d, e, f;
long prec1, prec2;
@ -51,8 +51,13 @@ int main()
fmpcb_init(e);
fmpcb_init(f);
fmpcb_randtest(a, state, 1 + n_randint(state, 1000), 200);
fmpcb_randtest(b, state, 1 + n_randint(state, 1000), 200);
fmpcb_randtest(a, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 200));
if (n_randint(state, 4) == 0)
fmprb_zero(fmpcb_imagref(a));
fmpcb_randtest(b, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 200));
if (n_randint(state, 4) == 0)
fmprb_zero(fmpcb_imagref(b));
fmpcb_pow(c, a, b, prec1);
fmpcb_pow(d, a, b, prec2);

103
fmprb/pow.c
View file

@ -25,12 +25,22 @@
#include "fmprb.h"
#define BINEXP_LIMIT 64
void
_fmprb_pow_exp(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
_fmprb_pow_exp(fmprb_t z, const fmprb_t x, int negx, const fmprb_t y, long prec)
{
fmprb_t t;
fmprb_init(t);
fmprb_log(t, x, prec);
if (negx)
{
fmprb_neg(t, x);
fmprb_log(t, t, prec);
}
else
fmprb_log(t, x, prec);
fmprb_mul(t, t, y, prec);
fmprb_exp(z, t, prec);
fmprb_clear(t);
@ -39,64 +49,51 @@ _fmprb_pow_exp(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
void
fmprb_pow(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
{
if (fmprb_is_exact(y))
const fmpr_struct * ymid = fmprb_midref(y);
const fmpr_struct * yrad = fmprb_radref(y);
if (fmprb_is_zero(y))
{
if (fmpr_is_zero(fmprb_midref(y)))
fmprb_one(z);
return;
}
if (fmpr_is_zero(yrad) && !fmpr_is_special(fmprb_midref(x)))
{
/* small half-integer or integer */
if (fmpr_cmpabs_2exp_si(ymid, BINEXP_LIMIT) < 0 &&
fmpr_is_int_2exp_si(ymid, -1))
{
fmprb_one(z);
fmpz_t e;
fmpz_init(e);
if (fmpr_is_int(ymid))
{
fmpr_get_fmpz_fixed_si(e, ymid, 0);
fmprb_pow_fmpz_binexp(z, x, e, prec);
}
else
{
fmpr_get_fmpz_fixed_si(e, ymid, -1);
fmprb_sqrt(z, x, prec + fmpz_bits(e));
fmprb_pow_fmpz_binexp(z, z, e, prec);
}
fmpz_clear(e);
return;
}
if (!fmpr_is_special(fmprb_midref(y)) && !fmpr_is_special(fmprb_midref(x)))
else if (fmpr_is_int(ymid) && fmpr_sgn(fmprb_midref(x)) < 0)
{
const fmpz * exp_exp = fmpr_expref(fmprb_midref(y));
/* smallish integer powers and square roots */
if (!COEFF_IS_MPZ(*exp_exp) && (*exp_exp >= -1L))
{
long exp_bits;
exp_bits = *exp_exp + fmpz_bits(fmpr_manref(fmprb_midref(y)));
if (exp_bits < 64)
{
fmpz_t e;
fmpz_init(e);
if (*exp_exp == -1L)
{
fmpz_set(e, fmpr_manref(fmprb_midref(y)));
fmprb_sqrt(z, x, prec + exp_bits);
fmprb_pow_fmpz_binexp(z, z, e, prec);
}
else
{
fmpz_mul_2exp(e, fmpr_manref(fmprb_midref(y)), *exp_exp);
fmprb_pow_fmpz_binexp(z, x, e, prec);
}
fmpz_clear(e);
return;
}
}
/* (-x)^n = (-1)^n * x^n */
if (fmpz_sgn(exp_exp) >= 0 && fmpr_sgn(fmprb_midref(x)) < 0)
{
fmprb_t t;
int odd;
fmprb_init(t);
fmprb_neg(t, x);
odd = fmpz_is_zero(exp_exp);
_fmprb_pow_exp(z, t, y, prec);
if (odd)
fmprb_neg(z, z);
fmprb_clear(t);
return;
}
/* use (-x)^n = (-1)^n * x^n to avoid NaNs
at least at high enough precision */
int odd = !fmpr_is_int_2exp_si(ymid, 1);
_fmprb_pow_exp(z, x, 1, y, prec);
if (odd)
fmprb_neg(z, z);
return;
}
}
_fmprb_pow_exp(z, x, y, prec);
_fmprb_pow_exp(z, x, 0, y, prec);
}

2
fmprb/test/t-pow.c
View file

@ -36,7 +36,7 @@ int main()
flint_randinit(state);
/* check large arguments */
for (iter = 0; iter < 10000; iter++)
for (iter = 0; iter < 20000; iter++)
{
fmprb_t a, b, c, d, e, f;
long prec1, prec2;

2
fmprb/test/t-trim.c
View file

@ -37,7 +37,7 @@ int main()
for (iter = 0; iter < 100000; iter++)
{
fmprb_t x, y;
long prec, acc1, acc2;
long acc1, acc2;
int accuracy_ok;
fmprb_init(x);