From 0847d21f66aba8373ffee194bcdd879aa4a9e8c8 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Fri, 7 Sep 2012 12:19:27 +0200 Subject: [PATCH] exp and log --- fmpr.h | 24 +++++++++++ fmpr/exp.c | 74 +++++++++++++++++++++++++++++++ fmpr/log.c | 57 +++++++++++------------- fmpr/test/t-exp.c | 103 ++++++++++++++++++++++++++++++++++++++++++++ fmpr/test/t-expm1.c | 103 ++++++++++++++++++++++++++++++++++++++++++++ fmpr/test/t-log1p.c | 103 ++++++++++++++++++++++++++++++++++++++++++++ fmprb.h | 2 + fmprb/exp.c | 61 ++++++++++++++++++++++++++ fmprb/log.c | 31 +++++++++++-- 9 files changed, 524 insertions(+), 34 deletions(-) create mode 100644 fmpr/exp.c create mode 100644 fmpr/test/t-exp.c create mode 100644 fmpr/test/t-expm1.c create mode 100644 fmpr/test/t-log1p.c create mode 100644 fmprb/exp.c diff --git a/fmpr.h b/fmpr.h index 427d8f59..124b7fe7 100644 --- a/fmpr.h +++ b/fmpr.h @@ -410,6 +410,10 @@ long fmpr_sqrt_ui(fmpr_t z, ulong x, long prec, fmpr_rnd_t rnd); long fmpr_sqrt_fmpz(fmpr_t z, const fmpz_t x, long prec, fmpr_rnd_t rnd); long fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd); +long fmpr_log1p(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd); + +long fmpr_exp(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd); +long fmpr_expm1(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd); static __inline__ void fmpr_neg(fmpr_t y, const fmpr_t x) @@ -529,5 +533,25 @@ void fmpr_get_fmpq(fmpq_t y, const fmpr_t x); long fmpr_set_fmpq(fmpr_t x, const fmpq_t y, long prec, fmpr_rnd_t rnd); +#define CALL_MPFR_FUNC(r, func, y, x, prec, rnd) \ + do { \ + mpfr_t t, u; \ + long r; \ + mpfr_rnd_t rrnd; \ + if (rnd == FMPR_RND_DOWN) rrnd = MPFR_RNDZ; \ + else if (rnd == FMPR_RND_UP) rrnd = MPFR_RNDA; \ + else if (rnd == FMPR_RND_FLOOR) rrnd = MPFR_RNDD; \ + else if (rnd == FMPR_RND_CEIL) rrnd = MPFR_RNDU; \ + else rrnd = MPFR_RNDN; \ + mpfr_init2(t, 2 + fmpz_bits(fmpr_manref(x))); \ + mpfr_init2(u, FLINT_MAX(2, prec)); \ + fmpr_get_mpfr(t, x, MPFR_RNDD); \ + func(u, t, rrnd); \ + fmpr_set_mpfr(y, u); \ + r = prec - fmpz_bits(fmpr_manref(y)); \ + mpfr_clear(t); \ + mpfr_clear(u); \ + } while (0); + #endif diff --git a/fmpr/exp.c b/fmpr/exp.c new file mode 100644 index 00000000..07ebf2fe --- /dev/null +++ b/fmpr/exp.c @@ -0,0 +1,74 @@ +/*============================================================================= + + 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 Fredrik Johansson + +******************************************************************************/ + +#include "fmpr.h" + +long +fmpr_exp(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd) +{ + if (fmpr_is_special(x)) + { + if (fmpr_is_zero(x)) + fmpr_one(y); + else if (fmpr_is_pos_inf(x)) + fmpr_pos_inf(y); + else if (fmpr_is_neg_inf(x)) + fmpr_zero(y); + else + fmpr_nan(y); + + return FMPR_RESULT_EXACT; + } + else + { + long r; + CALL_MPFR_FUNC(r, mpfr_exp, y, x, prec, rnd); + return r; + } +} + +long +fmpr_expm1(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd) +{ + if (fmpr_is_special(x)) + { + if (fmpr_is_zero(x)) + fmpr_zero(y); + else if (fmpr_is_pos_inf(x)) + fmpr_pos_inf(y); + else if (fmpr_is_neg_inf(x)) + fmpr_set_si(y, -1L); + else + fmpr_nan(y); + + return FMPR_RESULT_EXACT; + } + else + { + long r; + CALL_MPFR_FUNC(r, mpfr_expm1, y, x, prec, rnd); + return r; + } +} diff --git a/fmpr/log.c b/fmpr/log.c index 637b9055..7b8c24f9 100644 --- a/fmpr/log.c +++ b/fmpr/log.c @@ -25,16 +25,6 @@ #include "fmpr.h" -static mpfr_rnd_t fmpr_rnd_to_mpfr(fmpr_rnd_t rnd) -{ - if (rnd == FMPR_RND_NEAR) return MPFR_RNDN; - if (rnd == FMPR_RND_FLOOR) return MPFR_RNDD; - if (rnd == FMPR_RND_CEIL) return MPFR_RNDU; - if (rnd == FMPR_RND_DOWN) return MPFR_RNDZ; - if (rnd == FMPR_RND_UP) return MPFR_RNDA; - abort(); -} - long fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd) { @@ -49,37 +39,42 @@ fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd) return FMPR_RESULT_EXACT; } - - if (fmpr_sgn(x) < 0) + else if (fmpr_sgn(x) < 0) { fmpr_nan(y); return FMPR_RESULT_EXACT; } - - if (fmpr_is_one(x)) + else if (fmpr_is_one(x)) { fmpr_zero(y); return FMPR_RESULT_EXACT; } - + else { - mpfr_t t, u; long r; - - mpfr_init2(t, 1 + fmpz_bits(fmpr_manref(x))); - mpfr_init2(u, FLINT_MAX(2, prec)); - - fmpr_get_mpfr(t, x, MPFR_RNDD); - - mpfr_log(u, t, fmpr_rnd_to_mpfr(rnd)); - - fmpr_set_mpfr(y, u); - - r = prec - fmpz_bits(fmpr_manref(y)); - - mpfr_clear(t); - mpfr_clear(u); - + CALL_MPFR_FUNC(r, mpfr_log, y, x, prec, rnd); + return r; + } +} + +long +fmpr_log1p(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd) +{ + if (fmpr_is_special(x)) + { + if (fmpr_is_zero(x)) + fmpr_zero(y); + else if (fmpr_is_pos_inf(x)) + fmpr_pos_inf(y); + else + fmpr_nan(y); + + return FMPR_RESULT_EXACT; + } + else + { + long r; + CALL_MPFR_FUNC(r, mpfr_log1p, y, x, prec, rnd); return r; } } diff --git a/fmpr/test/t-exp.c b/fmpr/test/t-exp.c new file mode 100644 index 00000000..1ed19d80 --- /dev/null +++ b/fmpr/test/t-exp.c @@ -0,0 +1,103 @@ +/*============================================================================= + + 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 Fredrik Johansson + +******************************************************************************/ + +#include "fmpr.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("exp...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + long bits; + fmpr_t x, z, w; + mpfr_t X, Z; + + bits = 2 + n_randint(state, 200); + + fmpr_init(x); + fmpr_init(z); + fmpr_init(w); + + mpfr_init2(X, bits + 100); + mpfr_init2(Z, bits); + + fmpr_randtest_special(x, state, bits + n_randint(state, 100), 10); + fmpr_randtest_special(z, state, bits + n_randint(state, 100), 10); + + fmpr_get_mpfr(X, x, MPFR_RNDN); + + switch (n_randint(state, 4)) + { + case 0: + mpfr_exp(Z, X, MPFR_RNDZ); + fmpr_exp(z, x, bits, FMPR_RND_DOWN); + break; + case 1: + mpfr_exp(Z, X, MPFR_RNDA); + fmpr_exp(z, x, bits, FMPR_RND_UP); + break; + case 2: + mpfr_exp(Z, X, MPFR_RNDD); + fmpr_exp(z, x, bits, FMPR_RND_FLOOR); + break; + case 3: + mpfr_exp(Z, X, MPFR_RNDU); + fmpr_exp(z, x, bits, FMPR_RND_CEIL); + break; + } + + fmpr_set_mpfr(w, Z); + + if (!fmpr_equal(z, w)) + { + printf("FAIL\n\n"); + printf("bits = %ld\n", bits); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("z = "); fmpr_print(z); printf("\n\n"); + printf("w = "); fmpr_print(w); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(z); + fmpr_clear(w); + + mpfr_clear(X); + mpfr_clear(Z); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpr/test/t-expm1.c b/fmpr/test/t-expm1.c new file mode 100644 index 00000000..22615809 --- /dev/null +++ b/fmpr/test/t-expm1.c @@ -0,0 +1,103 @@ +/*============================================================================= + + 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 Fredrik Johansson + +******************************************************************************/ + +#include "fmpr.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("expm1...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + long bits; + fmpr_t x, z, w; + mpfr_t X, Z; + + bits = 2 + n_randint(state, 200); + + fmpr_init(x); + fmpr_init(z); + fmpr_init(w); + + mpfr_init2(X, bits + 100); + mpfr_init2(Z, bits); + + fmpr_randtest_special(x, state, bits + n_randint(state, 100), 10); + fmpr_randtest_special(z, state, bits + n_randint(state, 100), 10); + + fmpr_get_mpfr(X, x, MPFR_RNDN); + + switch (n_randint(state, 4)) + { + case 0: + mpfr_expm1(Z, X, MPFR_RNDZ); + fmpr_expm1(z, x, bits, FMPR_RND_DOWN); + break; + case 1: + mpfr_expm1(Z, X, MPFR_RNDA); + fmpr_expm1(z, x, bits, FMPR_RND_UP); + break; + case 2: + mpfr_expm1(Z, X, MPFR_RNDD); + fmpr_expm1(z, x, bits, FMPR_RND_FLOOR); + break; + case 3: + mpfr_expm1(Z, X, MPFR_RNDU); + fmpr_expm1(z, x, bits, FMPR_RND_CEIL); + break; + } + + fmpr_set_mpfr(w, Z); + + if (!fmpr_equal(z, w)) + { + printf("FAIL\n\n"); + printf("bits = %ld\n", bits); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("z = "); fmpr_print(z); printf("\n\n"); + printf("w = "); fmpr_print(w); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(z); + fmpr_clear(w); + + mpfr_clear(X); + mpfr_clear(Z); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpr/test/t-log1p.c b/fmpr/test/t-log1p.c new file mode 100644 index 00000000..c69f6af9 --- /dev/null +++ b/fmpr/test/t-log1p.c @@ -0,0 +1,103 @@ +/*============================================================================= + + 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 Fredrik Johansson + +******************************************************************************/ + +#include "fmpr.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("log1p...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + long bits; + fmpr_t x, z, w; + mpfr_t X, Z; + + bits = 2 + n_randint(state, 200); + + fmpr_init(x); + fmpr_init(z); + fmpr_init(w); + + mpfr_init2(X, bits + 100); + mpfr_init2(Z, bits); + + fmpr_randtest_special(x, state, bits + n_randint(state, 100), 10); + fmpr_randtest_special(z, state, bits + n_randint(state, 100), 10); + + fmpr_get_mpfr(X, x, MPFR_RNDN); + + switch (n_randint(state, 4)) + { + case 0: + mpfr_log1p(Z, X, MPFR_RNDZ); + fmpr_log1p(z, x, bits, FMPR_RND_DOWN); + break; + case 1: + mpfr_log1p(Z, X, MPFR_RNDA); + fmpr_log1p(z, x, bits, FMPR_RND_UP); + break; + case 2: + mpfr_log1p(Z, X, MPFR_RNDD); + fmpr_log1p(z, x, bits, FMPR_RND_FLOOR); + break; + case 3: + mpfr_log1p(Z, X, MPFR_RNDU); + fmpr_log1p(z, x, bits, FMPR_RND_CEIL); + break; + } + + fmpr_set_mpfr(w, Z); + + if (!fmpr_equal(z, w)) + { + printf("FAIL\n\n"); + printf("bits = %ld\n", bits); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("z = "); fmpr_print(z); printf("\n\n"); + printf("w = "); fmpr_print(w); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(z); + fmpr_clear(w); + + mpfr_clear(X); + mpfr_clear(Z); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmprb.h b/fmprb.h index 5d991407..85213a23 100644 --- a/fmprb.h +++ b/fmprb.h @@ -173,6 +173,8 @@ void fmprb_log(fmprb_t z, const fmprb_t x, long prec); 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_fac_ui(fmprb_t x, ulong n, long prec); void fmprb_bin_ui(fmprb_t x, const fmprb_t n, ulong k, long prec); diff --git a/fmprb/exp.c b/fmprb/exp.c new file mode 100644 index 00000000..65545172 --- /dev/null +++ b/fmprb/exp.c @@ -0,0 +1,61 @@ +/*============================================================================= + + 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 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb.h" + +void +fmprb_exp(fmprb_t z, const fmprb_t x, long prec) +{ + long r; + + 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); + } + else + { + /* + exp(a+b) - exp(a) = exp(a) * (exp(b)-1) + */ + fmpr_t t, u; + + 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_clear(t); + fmpr_clear(u); + } + + fmprb_adjust(z); +} diff --git a/fmprb/log.c b/fmprb/log.c index 57a13cbe..2f51b35f 100644 --- a/fmprb/log.c +++ b/fmprb/log.c @@ -30,12 +30,37 @@ fmprb_log(fmprb_t z, const fmprb_t x, long prec) { long r; - if (!fmprb_is_exact(x)) + if (fmprb_contains_zero(x) || fmpr_sgn(fmprb_midref(x)) < 0) + { + printf("log: interval contains zero or negative numbers\n"); abort(); + } - r = fmpr_log(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN); + if (fmprb_is_exact(x)) + { + r = fmpr_log(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN); + fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r); + } + else + { + /* + Let the input be [a-b, a+b]. We require a > b >= 0 (otherwise the + interval contains zero or a negative number and the logarithm is not + defined). The error is largest at a-b, and we have - fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r); + log(a) - log(a-b) = log(1 + b/(a-b)). + */ + fmpr_t err; + fmpr_init(err); + fmpr_sub(err, fmprb_midref(x), fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_DOWN); + fmpr_div(err, fmprb_radref(x), err, FMPRB_RAD_PREC, FMPR_RND_UP); + fmpr_log1p(err, err, FMPRB_RAD_PREC, FMPR_RND_UP); + + r = fmpr_log(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN); + fmpr_add_error_result(fmprb_radref(z), err, fmprb_midref(z), r, + FMPRB_RAD_PREC, FMPR_RND_UP); + fmpr_clear(err); + } fmprb_adjust(z); }