diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index fa8161e7..3cbf4236 100644 --- a/doc/source/fmprb.rst +++ b/doc/source/fmprb.rst @@ -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 ------------------------------------------------------------------------------- diff --git a/fmprb.h b/fmprb.h index ff1c5e72..cc4ba242 100644 --- a/fmprb.h +++ b/fmprb.h @@ -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); diff --git a/fmprb/add_error.c b/fmprb/add_error.c index 1640c645..c7d846b8 100644 --- a/fmprb/add_error.c +++ b/fmprb/add_error.c @@ -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) { diff --git a/fmprb/exp.c b/fmprb/exp.c index 65545172..7fff9d69 100644 --- a/fmprb/exp.c +++ b/fmprb/exp.c @@ -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); +} + diff --git a/fmprb/test/t-exp.c b/fmprb/test/t-exp.c index 4e7c8157..1a5c2a84 100644 --- a/fmprb/test/t-exp.c +++ b/fmprb/test/t-exp.c @@ -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"); diff --git a/fmprb/test/t-expm1.c b/fmprb/test/t-expm1.c new file mode 100644 index 00000000..174064d4 --- /dev/null +++ b/fmprb/test/t-expm1.c @@ -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; +}