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

fix fmprb_exp to support huge exponents; add fmprb_expm1

Browse source
This commit is contained in:
Fredrik Johansson 2013-03-24 18:46:04 +01:00
parent fb6e85669c
commit 4f102eb83d
6 changed files with 450 additions and 26 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
doc/source/fmprb.rst
View file

@ -214,6 +214,8 @@ Precision and comparisons
.. function:: void fmprb_add_error_2exp_si(fmprb_t x, long e)
.. function:: void fmprb_add_error_2exp_fmpz(fmprb_t x, const fmpz_t e)
Adds `2^e` to the radius of *x*.
.. function:: void fmprb_add_error(fmprb_t x, const fmprb_t err)
@ -361,6 +363,8 @@ Arithmetic
.. function:: void fmprb_mul_2exp_si(fmprb_t y, const fmprb_t x, long e)
.. function:: void fmprb_mul_2exp_fmpz(fmprb_t y, const fmprb_t x, const fmpz_t e)
Sets *y* to *x* multiplied by `2^e`.
.. function:: void fmprb_div(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
@ -479,10 +483,11 @@ Exponentials and logarithms
Sets `z = \exp(x)`. Error propagation is done using the following rule:
the error is largest at `m + r`, and we have
`\exp(m+r) - \exp(m) = \exp(m) (\exp(r)-1)`.
The last expression is calculated accurately for small radii
via *fmpr_expm1*.
`\exp(m+r) - \exp(m) = \exp(m) (\exp(r)-1) \le r \exp(m+r)`.
.. function:: void fmprb_expm1(fmprb_t z, const fmprb_t x, long prec)
Sets `z = \exp(x)-1`, computed accurately when `x \approx 0`.
Trigonometric functions
-------------------------------------------------------------------------------

9
fmprb.h
View file

@ -272,6 +272,7 @@ void fmprb_log_ui(fmprb_t z, ulong x, long prec);
void fmprb_log_fmpz(fmprb_t z, const fmpz_t x, long prec);
void fmprb_exp(fmprb_t z, const fmprb_t x, long prec);
void fmprb_expm1(fmprb_t z, const fmprb_t x, long prec);
void fmprb_sin(fmprb_t s, const fmprb_t x, long prec);
void fmprb_cos(fmprb_t c, const fmprb_t x, long prec);
@ -351,6 +352,13 @@ fmprb_mul_2exp_si(fmprb_t y, const fmprb_t x, long e)
fmpr_mul_2exp_si(fmprb_radref(y), fmprb_radref(x), e);
}
static __inline__ void
fmprb_mul_2exp_fmpz(fmprb_t y, const fmprb_t x, const fmpz_t e)
{
fmpr_mul_2exp_fmpz(fmprb_midref(y), fmprb_midref(x), e);
fmpr_mul_2exp_fmpz(fmprb_radref(y), fmprb_radref(x), e);
}
static __inline__ void
fmprb_set_fmpq(fmprb_t y, const fmpq_t x, long prec)
{
@ -507,6 +515,7 @@ fmprb_bits(const fmprb_t x)
void fmprb_add_error_fmpr(fmprb_t x, const fmpr_t err);
void fmprb_add_error_2exp_si(fmprb_t x, long err);
void fmprb_add_error_2exp_fmpz(fmprb_t x, const fmpz_t err);
void fmprb_add_error(fmprb_t x, const fmprb_t error);
void fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits);

11
fmprb/add_error.c
View file

@ -41,6 +41,17 @@ fmprb_add_error_2exp_si(fmprb_t x, long err)
fmpr_clear(t);
}
void
fmprb_add_error_2exp_fmpz(fmprb_t x, const fmpz_t err)
{
fmpr_t t;
fmpr_init(t);
fmpz_one(fmpr_manref(t));
fmpz_set(fmpr_expref(t), err);
fmprb_add_error_fmpr(x, t);
fmpr_clear(t);
}
void
fmprb_add_error(fmprb_t x, const fmprb_t error)
{

252
fmprb/exp.c
View file

@ -19,43 +19,251 @@
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
#include "fmprb.h"
void
fmprb_exp(fmprb_t z, const fmprb_t x, long prec)
fmpr_exp_ubound(fmpr_t y, const fmpr_t x, long prec, long maglim)
{
long r;
fmpz_t mag;
if (fmpr_is_special(x))
{
if (fmpr_is_pos_inf(x))
fmpr_pos_inf(y);
else if (fmpr_is_zero(x))
fmpr_one(y);
else if (fmpr_is_neg_inf(x))
fmpr_zero(y);
else
fmpr_pos_inf(y);
return;
}
fmpz_init(mag);
/* 2^(mag-1) <= |x| < 2^mag */
fmpz_add_ui(mag, fmpr_expref(x), fmpz_bits(fmpr_manref(x)));
/* normal range -- ok to just call mpfr */
if (!COEFF_IS_MPZ(*mag) && (FLINT_ABS(*mag) <= 24))
{
fmpr_exp(y, x, prec, FMPR_RND_UP);
}
/* close to zero */
else if (fmpz_sgn(mag) < 0)
{
/* x < 0 ==> exp(x) < 1 */
if (fmpz_sgn(fmpr_manref(x)) < 0)
{
fmpr_one(y);
}
/* x > 0 ==> exp(x) < 1 + 2x */
else
{
fmpr_mul_2exp_si(y, x, 1);
fmpr_add_ui(y, y, 1, prec, FMPR_RND_UP);
}
}
/* overflow */
else if (COEFF_IS_MPZ(*mag) || *mag > maglim)
{
if (fmpz_sgn(fmpr_manref(x)) > 0)
{
fmpr_pos_inf(y);
}
else
{
/*
mag > maglim
mag - 1 >= maglim
x <= -2^(mag-1)
x <= -2^maglim
exp(x) < 2^(-2^(mag-1)) <= 2^(-2^maglim)
*/
fmpz_set_si(fmpr_expref(y), -1);
fmpz_mul_2exp(fmpr_expref(y), fmpr_expref(y), maglim);
fmpz_one(fmpr_manref(y));
}
}
/* +x or -x is huge */
else
{
fmprb_t ln2, t, u;
fmpz_t q;
long wp;
fmprb_init(ln2);
fmprb_init(t);
fmprb_init(u);
fmpz_init(q);
wp = prec + *mag + 10;
fmprb_const_log2(ln2, wp);
fmprb_set_fmpr(t, x);
fmprb_div(u, t, ln2, wp);
fmpr_get_fmpz(q, fmprb_midref(u), FMPR_RND_DOWN);
fmprb_submul_fmpz(t, ln2, q, wp);
fmpr_add(y, fmprb_midref(t), fmprb_radref(t), prec, FMPR_RND_CEIL);
fmpr_exp(y, y, prec, FMPR_RND_UP);
fmpr_mul_2exp_fmpz(y, y, q);
fmprb_clear(ln2);
fmprb_clear(t);
fmprb_clear(u);
fmpz_clear(q);
}
fmpz_clear(mag);
}
void
fmprb_exp_fmpr(fmprb_t z, const fmpr_t x, long prec, long maglim, int m1)
{
long r, small_cutoff;
fmpz_t mag;
if (fmpr_is_special(x))
{
if (m1)
fmpr_expm1(fmprb_midref(z), x, prec, FMPR_RND_DOWN);
else
fmpr_exp(fmprb_midref(z), x, prec, FMPR_RND_DOWN);
fmpr_zero(fmprb_radref(z));
return;
}
fmpz_init(mag);
/* 2^(mag-1) <= |x| < 2^mag */
fmpz_add_ui(mag, fmpr_expref(x), fmpz_bits(fmpr_manref(x)));
/* magnitude for approximation by 1 + x */
if (m1)
small_cutoff = -FLINT_MAX(24, 4 + prec);
else
small_cutoff = -FLINT_MAX(24, 4 + prec / 2);
/* standard case */
if (!COEFF_IS_MPZ(*mag) && (*mag > small_cutoff) && (*mag <= 24))
{
if (m1)
r = fmpr_expm1(fmprb_midref(z), x, prec, FMPR_RND_DOWN);
else
r = fmpr_exp(fmprb_midref(z), x, prec, FMPR_RND_DOWN);
fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
}
/* close to zero */
else if (fmpz_sgn(mag) < 0)
{
/* exp(x) = 1 + x + eps, eps < x^2 < 2^(2mag) */
r = fmpr_add_ui(fmprb_midref(z), x, m1 ? 0 : 1, prec, FMPR_RND_DOWN);
fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
fmpz_mul_2exp(mag, mag, 1);
fmprb_add_error_2exp_fmpz(z, mag);
}
/* overflow */
else if (COEFF_IS_MPZ(*mag) || *mag > maglim)
{
if (fmpz_sgn(fmpr_manref(x)) > 0)
{
fmpr_pos_inf(fmprb_midref(z));
fmpr_pos_inf(fmprb_radref(z));
}
else
{
/*
mag > maglim
mag - 1 >= maglim
x <= -2^(mag-1)
x <= -2^maglim
0 < exp(x) < 2^(-2^(mag-1)) <= 2^(-2^maglim)
*/
fmpz_set_si(fmpr_expref(fmprb_midref(z)), -1);
fmpz_mul_2exp(fmpr_expref(fmprb_midref(z)), fmpr_expref(fmprb_midref(z)), maglim);
fmpz_one(fmpr_manref(fmprb_midref(z)));
fmpr_set(fmprb_radref(z), fmprb_midref(z));
if (m1)
fmprb_sub_ui(z, z, 1, prec);
}
}
/* huge */
else
{
fmprb_t ln2, t, u;
fmpz_t q;
long wp;
fmprb_init(ln2);
fmprb_init(t);
fmprb_init(u);
fmpz_init(q);
wp = prec + *mag + 10;
fmprb_const_log2(ln2, wp);
fmprb_set_fmpr(t, x);
fmprb_div(u, t, ln2, wp);
fmpr_get_fmpz(q, fmprb_midref(u), FMPR_RND_DOWN);
fmprb_submul_fmpz(t, ln2, q, wp);
fmprb_exp(z, t, prec);
fmprb_mul_2exp_fmpz(z, z, q);
if (m1)
fmprb_sub_ui(z, z, 1, prec);
fmprb_clear(ln2);
fmprb_clear(t);
fmprb_clear(u);
fmpz_clear(q);
}
fmpz_clear(mag);
}
void
_fmprb_exp(fmprb_t z, const fmprb_t x, long prec, int m1)
{
long maglim = FLINT_MAX(128, 2 * prec);
if (fmprb_is_exact(x))
{
r = fmpr_exp(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN);
fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
fmprb_exp_fmpr(z, fmprb_midref(x), prec, maglim, m1);
}
else
{
/*
exp(a+b) - exp(a) = exp(a) * (exp(b)-1)
*/
fmpr_t t, u;
/* exp(a+b) - exp(a) = exp(a) * (exp(b)-1) <= b * exp(a+b) */
fmpr_t t;
fmpr_init(t);
fmpr_init(u);
fmpr_exp(t, fmprb_midref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_expm1(u, fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul(t, t, u, FMPRB_RAD_PREC, FMPR_RND_UP);
r = fmpr_exp(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN);
fmpr_add_error_result(fmprb_radref(z), t, fmprb_midref(z), r,
FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_add(t, fmprb_midref(x), fmprb_radref(x), prec, FMPR_RND_CEIL);
fmpr_exp_ubound(t, t, FMPRB_RAD_PREC, maglim);
fmpr_mul(t, t, fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
fmprb_exp_fmpr(z, fmprb_midref(x), prec, maglim, m1);
fmpr_add(fmprb_radref(z), fmprb_radref(z), t, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_clear(t);
fmpr_clear(u);
}
fmprb_adjust(z);
}
void
fmprb_exp(fmprb_t z, const fmprb_t x, long prec)
{
_fmprb_exp(z, x, prec, 0);
}
void
fmprb_expm1(fmprb_t z, const fmprb_t x, long prec)
{
_fmprb_exp(z, x, prec, 1);
}

56
fmprb/test/t-exp.c
View file

@ -19,7 +19,7 @@
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
@ -78,6 +78,60 @@ int main()
mpfr_clear(t);
}
/* check large arguments */
for (iter = 0; iter < 100000; iter++)
{
fmprb_t a, b, c, d;
long prec1, prec2;
prec1 = 2 + n_randint(state, 1000);
prec2 = prec1 + 30;
fmprb_init(a);
fmprb_init(b);
fmprb_init(c);
fmprb_init(d);
fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 100);
fmprb_exp(b, a, prec1);
fmprb_exp(c, a, prec2);
if (!fmprb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
fmprb_randtest_precise(b, state, 1 + n_randint(state, 1000), 100);
/* check exp(a)*exp(b) = exp(a+b) */
fmprb_exp(c, a, prec1);
fmprb_exp(d, b, prec1);
fmprb_mul(c, c, d, prec1);
fmprb_add(d, a, b, prec1);
fmprb_exp(d, d, prec1);
if (!fmprb_overlaps(c, d))
{
printf("FAIL: functional equation\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
printf("d = "); fmprb_print(d); printf("\n\n");
abort();
}
fmprb_clear(a);
fmprb_clear(b);
fmprb_clear(c);
fmprb_clear(d);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");

137
fmprb/test/t-expm1.c Normal file
View file

@ -0,0 +1,137 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
#include "fmprb.h"
int main()
{
long iter;
flint_rand_t state;
printf("expm1....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmprb_t a, b;
fmpq_t q;
mpfr_t t;
long prec = 2 + n_randint(state, 200);
fmprb_init(a);
fmprb_init(b);
fmpq_init(q);
mpfr_init2(t, prec + 100);
fmprb_randtest(a, state, 1 + n_randint(state, 200), 3);
fmprb_randtest(b, state, 1 + n_randint(state, 200), 3);
fmprb_get_rand_fmpq(q, state, a, 1 + n_randint(state, 200));
fmpq_get_mpfr(t, q, MPFR_RNDN);
mpfr_expm1(t, t, MPFR_RNDN);
fmprb_expm1(b, a, prec);
if (!fmprb_contains_mpfr(b, t))
{
printf("FAIL: containment\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
abort();
}
fmprb_expm1(a, a, prec);
if (!fmprb_equal(a, b))
{
printf("FAIL: aliasing\n\n");
abort();
}
fmprb_clear(a);
fmprb_clear(b);
fmpq_clear(q);
mpfr_clear(t);
}
/* check large arguments */
for (iter = 0; iter < 100000; iter++)
{
fmprb_t a, b, c, d;
long prec1, prec2;
prec1 = 2 + n_randint(state, 1000);
prec2 = prec1 + 30;
fmprb_init(a);
fmprb_init(b);
fmprb_init(c);
fmprb_init(d);
fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 100);
fmprb_expm1(b, a, prec1);
fmprb_expm1(c, a, prec2);
if (!fmprb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
fmprb_randtest_precise(b, state, 1 + n_randint(state, 1000), 100);
/* compare with exp */
fmprb_expm1(c, a, prec1);
fmprb_exp(d, a, prec1);
fmprb_sub_ui(d, d, 1, prec1);
if (!fmprb_overlaps(c, d))
{
printf("FAIL: comparison with exp\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
printf("d = "); fmprb_print(d); printf("\n\n");
abort();
}
fmprb_clear(a);
fmprb_clear(b);
fmprb_clear(c);
fmprb_clear(d);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}