From 5b4e35b750a8fb41a7bf8e06450244f533620615 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 5 Feb 2013 15:10:27 +0100 Subject: [PATCH] tidy up float conversion functions --- doc/source/fmpr.rst | 4 ++ fmpr.h | 75 +++++++------------------------------- fmpr/get_d.c | 43 ++++++++++++++++++++++ fmpr/get_mpfr.c | 62 +++++++++++++++++++++++++++++++ fmpr/set_mpfr.c | 46 +++++++++++++++++++++++ fmprb/gamma.c | 18 ++------- hypgeom/bound.c | 4 +- hypgeom/estimate_terms_d.c | 4 +- 8 files changed, 174 insertions(+), 82 deletions(-) create mode 100644 fmpr/get_d.c create mode 100644 fmpr/get_mpfr.c create mode 100644 fmpr/set_mpfr.c diff --git a/doc/source/fmpr.rst b/doc/source/fmpr.rst index 71ea73b1..143f82cf 100644 --- a/doc/source/fmpr.rst +++ b/doc/source/fmpr.rst @@ -192,6 +192,10 @@ Assignment, rounding and conversions Sets *x* to the exact value of the MPFR variable *y*. +.. function:: double fmpr_get_d(const fmpr_t x, fmpr_rnd_t rnd) + + Returns *x* rounded to a *double* in the direction specified by *rnd*. + .. function:: void fmpr_set_ui(fmpr_t x, ulong c) .. function:: void fmpr_set_si(fmpr_t x, long c) diff --git a/fmpr.h b/fmpr.h index 1ac8512f..bba453ff 100644 --- a/fmpr.h +++ b/fmpr.h @@ -41,6 +41,15 @@ #define FMPR_RND_DOWN 3 #define FMPR_RND_NEAR 4 +static __inline__ mpfr_rnd_t rnd_to_mpfr(fmpr_rnd_t rnd) +{ + if (rnd == FMPR_RND_DOWN) return MPFR_RNDZ; + else if (rnd == FMPR_RND_UP) return MPFR_RNDA; + else if (rnd == FMPR_RND_FLOOR) return MPFR_RNDD; + else if (rnd == FMPR_RND_CEIL) return MPFR_RNDU; + else return MPFR_RNDN; +} + typedef struct { fmpz man; @@ -277,59 +286,11 @@ void fmpr_randtest_not_zero(fmpr_t x, flint_rand_t state, long bits, long exp_bi void fmpr_randtest_special(fmpr_t x, flint_rand_t state, long bits, long exp_bits); -static __inline__ int -fmpr_get_mpfr(mpfr_t x, const fmpr_t y, mpfr_rnd_t rnd) -{ - int r; +int fmpr_get_mpfr(mpfr_t x, const fmpr_t y, mpfr_rnd_t rnd); - if (fmpr_is_special(y)) - { - if (fmpr_is_zero(y)) mpfr_set_zero(x, 0); - else if (fmpr_is_pos_inf(y)) mpfr_set_inf(x, 1); - else if (fmpr_is_neg_inf(y)) mpfr_set_inf(x, -1); - else mpfr_set_nan(x); - r = 0; - } - else if (COEFF_IS_MPZ(*fmpr_expref(y))) - { - printf("exception: exponent too large to convert to mpfr"); - abort(); - } - else - { - if (!COEFF_IS_MPZ(*fmpr_manref(y))) - r = mpfr_set_si_2exp(x, *fmpr_manref(y), *fmpr_expref(y), rnd); - else - r = mpfr_set_z_2exp(x, COEFF_TO_PTR(*fmpr_manref(y)), *fmpr_expref(y), rnd); +void fmpr_set_mpfr(fmpr_t x, const mpfr_t y); - if (!mpfr_regular_p(x)) - { - printf("exception: exponent too large to convert to mpfr"); - abort(); - } - } - - return r; -} - -static __inline__ void -fmpr_set_mpfr(fmpr_t x, const mpfr_t y) -{ - if (!mpfr_regular_p(y)) - { - if (mpfr_zero_p(y)) fmpr_zero(x); - else if (mpfr_inf_p(y) && mpfr_sgn(y) > 0) fmpr_pos_inf(x); - else if (mpfr_inf_p(y) && mpfr_sgn(y) < 0) fmpr_neg_inf(x); - else fmpr_nan(x); - } - else - { - __mpz_struct * z = _fmpz_promote(fmpr_manref(x)); - fmpz_set_si(fmpr_expref(x), mpfr_get_z_2exp(z, y)); - _fmpz_demote_val(fmpr_manref(x)); - _fmpr_normalise(fmpr_manref(x), fmpr_expref(x), y->_mpfr_prec, FMPR_RND_DOWN); - } -} +double fmpr_get_d(const fmpr_t x, fmpr_rnd_t rnd); static __inline__ void fmpr_set_ui(fmpr_t x, ulong c) { @@ -671,11 +632,7 @@ static __inline__ int fmpr_cmpabs_2exp_si(const fmpr_t x, long e) do { \ mpfr_t __t, __u; \ mpfr_rnd_t __rnd; \ - if (rnd == FMPR_RND_DOWN) __rnd = MPFR_RNDZ; \ - else if (rnd == FMPR_RND_UP) __rnd = MPFR_RNDA; \ - else if (rnd == FMPR_RND_FLOOR) __rnd = MPFR_RNDD; \ - else if (rnd == FMPR_RND_CEIL) __rnd = MPFR_RNDU; \ - else __rnd = MPFR_RNDN; \ + __rnd = rnd_to_mpfr(rnd); \ mpfr_init2(__t, 2 + fmpz_bits(fmpr_manref(x))); \ mpfr_init2(__u, FLINT_MAX(2, prec)); \ fmpr_get_mpfr(__t, x, MPFR_RNDD); \ @@ -695,11 +652,7 @@ static __inline__ int fmpr_cmpabs_2exp_si(const fmpr_t x, long e) do { \ mpfr_t __t, __u, __v; \ mpfr_rnd_t __rnd; \ - if (rnd == FMPR_RND_DOWN) __rnd = MPFR_RNDZ; \ - else if (rnd == FMPR_RND_UP) __rnd = MPFR_RNDA; \ - else if (rnd == FMPR_RND_FLOOR) __rnd = MPFR_RNDD; \ - else if (rnd == FMPR_RND_CEIL) __rnd = MPFR_RNDU; \ - else __rnd = MPFR_RNDN; \ + __rnd = rnd_to_mpfr(rnd); \ mpfr_init2(__t, 2 + fmpz_bits(fmpr_manref(x))); \ mpfr_init2(__u, FLINT_MAX(2, prec)); \ mpfr_init2(__v, FLINT_MAX(2, prec)); \ diff --git a/fmpr/get_d.c b/fmpr/get_d.c new file mode 100644 index 00000000..481e435b --- /dev/null +++ b/fmpr/get_d.c @@ -0,0 +1,43 @@ +/*============================================================================= + + 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" + +double +fmpr_get_d(const fmpr_t x, fmpr_rnd_t rnd) +{ + double r; + mpfr_rnd_t mrnd; + mpfr_t t; + + mrnd = rnd_to_mpfr(rnd); + mpfr_init2(t, 53); + fmpr_get_mpfr(t, x, mrnd); + r = mpfr_get_d(t, mrnd); + + mpfr_clear(t); + return r; +} + diff --git a/fmpr/get_mpfr.c b/fmpr/get_mpfr.c new file mode 100644 index 00000000..7fc472a7 --- /dev/null +++ b/fmpr/get_mpfr.c @@ -0,0 +1,62 @@ +/*============================================================================= + + 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 +fmpr_get_mpfr(mpfr_t x, const fmpr_t y, mpfr_rnd_t rnd) +{ + int r; + + if (fmpr_is_special(y)) + { + if (fmpr_is_zero(y)) mpfr_set_zero(x, 0); + else if (fmpr_is_pos_inf(y)) mpfr_set_inf(x, 1); + else if (fmpr_is_neg_inf(y)) mpfr_set_inf(x, -1); + else mpfr_set_nan(x); + r = 0; + } + else if (COEFF_IS_MPZ(*fmpr_expref(y))) + { + printf("exception: exponent too large to convert to mpfr"); + abort(); + } + else + { + if (!COEFF_IS_MPZ(*fmpr_manref(y))) + r = mpfr_set_si_2exp(x, *fmpr_manref(y), *fmpr_expref(y), rnd); + else + r = mpfr_set_z_2exp(x, COEFF_TO_PTR(*fmpr_manref(y)), *fmpr_expref(y), rnd); + + if (!mpfr_regular_p(x)) + { + printf("exception: exponent too large to convert to mpfr"); + abort(); + } + } + + return r; +} + diff --git a/fmpr/set_mpfr.c b/fmpr/set_mpfr.c new file mode 100644 index 00000000..ae4e9428 --- /dev/null +++ b/fmpr/set_mpfr.c @@ -0,0 +1,46 @@ +/*============================================================================= + + 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" + +void +fmpr_set_mpfr(fmpr_t x, const mpfr_t y) +{ + if (!mpfr_regular_p(y)) + { + if (mpfr_zero_p(y)) fmpr_zero(x); + else if (mpfr_inf_p(y) && mpfr_sgn(y) > 0) fmpr_pos_inf(x); + else if (mpfr_inf_p(y) && mpfr_sgn(y) < 0) fmpr_neg_inf(x); + else fmpr_nan(x); + } + else + { + __mpz_struct * z = _fmpz_promote(fmpr_manref(x)); + fmpz_set_si(fmpr_expref(x), mpfr_get_z_2exp(z, y)); + _fmpz_demote_val(fmpr_manref(x)); + _fmpr_normalise(fmpr_manref(x), fmpr_expref(x), y->_mpfr_prec, FMPR_RND_DOWN); + } +} + diff --git a/fmprb/gamma.c b/fmprb/gamma.c index 96485905..0d75a2ec 100644 --- a/fmprb/gamma.c +++ b/fmprb/gamma.c @@ -39,18 +39,6 @@ fmprb_stirling_series_coeff(fmprb_t b, ulong k, long prec) fmpz_clear(d); } - -double fmpr_get_d(const fmpr_t x) -{ - double r; - mpfr_t t; - mpfr_init2(t, 53); - fmpr_get_mpfr(t, x, MPFR_RNDN); - r = mpfr_get_d(t, MPFR_RNDN); - mpfr_clear(t); - return r; -} - /* Heuristic: we use Stirling's series if abs(x) > beta * prec. @@ -63,7 +51,7 @@ argument reduction. long stirling_choose_r(const fmprb_t x, long wp) { - double t = fmpr_get_d(fmprb_midref(x)); + double t = fmpr_get_d(fmprb_midref(x), FMPR_RND_FLOOR); double want = FLINT_MAX(5, GAMMA_STIRLING_BETA * wp); return (long) FLINT_MAX(0, want - t + 1); @@ -77,7 +65,7 @@ stirling_choose_nterms(const fmprb_t x, long r, double bits) double t, logt; double mag; - t = fmpr_get_d(fmprb_midref(x)) + r; + t = fmpr_get_d(fmprb_midref(x), FMPR_RND_FLOOR) + r; logt = log(t); for (i = 1; ; i++) @@ -112,7 +100,7 @@ fmprb_gamma_log_stirling(fmprb_t s, const fmprb_t z, long nterms, long prec) { fmprb_ui_div(t, 1UL, z, prec); fmprb_mul(u, t, t, prec); - z_mag = fmpr_get_d(fmprb_midref(w)) * 1.44269504088896; + z_mag = fmpr_get_d(fmprb_midref(w), FMPR_RND_UP) * 1.44269504088896; for (k = nterms - 1; k >= 1; k--) { diff --git a/hypgeom/bound.c b/hypgeom/bound.c index 2801d3be..106b15a6 100644 --- a/hypgeom/bound.c +++ b/hypgeom/bound.c @@ -27,8 +27,6 @@ #include "double_extras.h" #include "hypgeom.h" -double fmpr_get_d(const fmpr_t x); - static __inline__ double d_root(double x, int r) { if (r == 1) @@ -124,7 +122,7 @@ long hypgeom_root_bound(const fmpr_t z, int r) else { double zd; - zd = fmpr_get_d(z); + zd = fmpr_get_d(z, FMPR_RND_UP); return d_root(zd, r) + 2; } } diff --git a/hypgeom/estimate_terms_d.c b/hypgeom/estimate_terms_d.c index a8c4dcda..4aa6e41b 100644 --- a/hypgeom/estimate_terms_d.c +++ b/hypgeom/estimate_terms_d.c @@ -39,14 +39,12 @@ static __inline__ double d_root(double x, int r) return pow(x, 1. / r); } -double fmpr_get_d(const fmpr_t x); - long hypgeom_estimate_terms(const fmpr_t z, int r, long prec) { double y, t; - t = fmpr_get_d(z); + t = fmpr_get_d(z, FMPR_RND_UP); t = fabs(t); if (t == 0)