From 147297748f050387bf6aa02e627b7f9ea246cf29 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sat, 22 Dec 2012 21:29:11 +0100 Subject: [PATCH] refine fmprb test functions; add gamma/rgamma/lgamma with test code --- doc/source/fmprb.rst | 27 ++++++++++++--- fmprb.h | 8 ++++- fmprb/gamma.c | 55 ++++++++++++++++++++++++++++- fmprb/randtest.c | 57 +++++++++++++++++++++--------- fmprb/test/t-gamma.c | 81 +++++++++++++++++++++++++++++++++++++++++++ fmprb/test/t-lgamma.c | 81 +++++++++++++++++++++++++++++++++++++++++++ fmprb/test/t-rgamma.c | 81 +++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 367 insertions(+), 23 deletions(-) create mode 100644 fmprb/test/t-gamma.c create mode 100644 fmprb/test/t-lgamma.c create mode 100644 fmprb/test/t-rgamma.c diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index 968eee57..a6679404 100644 --- a/doc/source/fmprb.rst +++ b/doc/source/fmprb.rst @@ -178,6 +178,20 @@ Random number generation Generates a random ball. The midpoint and radius will both be finite. +.. function:: void fmprb_randtest_exact(fmprb_t x, flint_rand_t state, long prec, long mag_bits) + + Generates a random number with zero radius. + +.. function:: void fmprb_randtest_precise(fmprb_t x, flint_rand_t state, long prec, long mag_bits) + + Generates a random number with radius at most `2^{-\mathrm{prec}}` + the magnitude of the midpoint. + +.. function:: void fmprb_randtest_wide(fmprb_t x, flint_rand_t state, long prec, long mag_bits) + + Generates a random number with midpoint and radius chosen independently, + possibly giving a very large interval. + .. function:: void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, long bits) Sets *q* to a random rational number from the interval represented by *x*. @@ -723,9 +737,14 @@ Gamma function be done here so that a working precision is selected which always is sufficient and also nearly optimal. -.. function:: void fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec) +.. function:: void fmprb_gamma(fmprb_t y, const fmprb_t x, long prec) - Sets `y = \log \Gamma(x)`, assuming that `x > 0`. +.. function:: void fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec) + +.. function:: void fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec) + + Sets, respectively, `y = \Gamma(x)`, `y = 1/\Gamma(x)`, + `y = \log \Gamma(x)`. For large `x`, uses Stirling's expansion @@ -744,8 +763,8 @@ Gamma function If `x` is too small for the asymptotic expansion to give sufficient accuracy directly, we translate to `x + r` - using the formula `\log \Gamma(x) = \log \Gamma(x+r) - - \log(x (x+1) (x+2) \cdots (x+r-1))`. + using the formula `\Gamma(x) = \Gamma(x+r) / + (x (x+1) (x+2) \cdots (x+r-1))`. To obtain a remainder smaller than `2^{-b}`, we must choose an `r` such that `x + r > \beta b`, where `\beta > \log(2) / (2 \pi) \approx 0.11`. diff --git a/fmprb.h b/fmprb.h index a165a537..27bb9d34 100644 --- a/fmprb.h +++ b/fmprb.h @@ -316,7 +316,10 @@ void fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec); void fmprb_zeta_ui_vec(fmprb_struct * x, ulong start, long num, long prec); void fmprb_gamma_fmpq_karatsuba(fmprb_struct * v, const fmpq_t a, long num, long prec); -void fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec); + +void fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec); +void fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec); +void fmprb_gamma(fmprb_t y, const fmprb_t x, long prec); static __inline__ void fmprb_print(const fmprb_t x) @@ -493,6 +496,9 @@ void fmprb_add_error_2exp_si(fmprb_t x, long 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); +void fmprb_randtest_exact(fmprb_t x, flint_rand_t state, long prec, long mag_bits); +void fmprb_randtest_wide(fmprb_t x, flint_rand_t state, long prec, long mag_bits); +void fmprb_randtest_precise(fmprb_t x, flint_rand_t state, long prec, long mag_bits); void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, long bits); diff --git a/fmprb/gamma.c b/fmprb/gamma.c index f5ee82d7..a38738b5 100644 --- a/fmprb/gamma.c +++ b/fmprb/gamma.c @@ -159,7 +159,7 @@ fmprb_gamma_log_stirling(fmprb_t s, const fmprb_t z, long nterms, long prec) } void -fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec) +fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec) { long r, n, wp; fmprb_t t, u; @@ -183,3 +183,56 @@ fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec) fmprb_clear(t); fmprb_clear(u); } + +void +fmprb_gamma(fmprb_t y, const fmprb_t x, long prec) +{ + long r, n, wp; + fmprb_t t, u; + + wp = prec + FLINT_BIT_COUNT(prec); + + r = stirling_choose_r(x, wp); + n = stirling_choose_nterms(x, r, wp); + + /* gamma(x) = gamma(x+r) / rf(x,r) */ + fmprb_init(t); + fmprb_init(u); + + fmprb_add_ui(t, x, r, wp); + fmprb_gamma_log_stirling(u, t, n, wp); + fmprb_exp(u, u, prec); + + fmprb_rfac_ui_bsplit(t, x, r, wp); + fmprb_div(y, u, t, prec); + + fmprb_clear(t); + fmprb_clear(u); +} + +void +fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec) +{ + long r, n, wp; + fmprb_t t, u; + + wp = prec + FLINT_BIT_COUNT(prec); + + r = stirling_choose_r(x, wp); + n = stirling_choose_nterms(x, r, wp); + + /* 1/gamma(x) = rf(x,r) / gamma(x+r) */ + fmprb_init(t); + fmprb_init(u); + + fmprb_add_ui(t, x, r, wp); + fmprb_gamma_log_stirling(u, t, n, wp); + fmprb_neg(u, u); + fmprb_exp(u, u, prec); + + fmprb_rfac_ui_bsplit(t, x, r, wp); + fmprb_mul(y, u, t, prec); + + fmprb_clear(t); + fmprb_clear(u); +} diff --git a/fmprb/randtest.c b/fmprb/randtest.c index e2e60cf2..e31382c3 100644 --- a/fmprb/randtest.c +++ b/fmprb/randtest.c @@ -26,26 +26,49 @@ #include "fmprb.h" void -fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits) +fmprb_randtest_exact(fmprb_t x, flint_rand_t state, long prec, long mag_bits) +{ + fmpr_randtest(fmprb_midref(x), state, prec, mag_bits); + fmpr_zero(fmprb_radref(x)); +} + +void +fmprb_randtest_wide(fmprb_t x, flint_rand_t state, long prec, long mag_bits) +{ + fmpr_randtest(fmprb_midref(x), state, prec, mag_bits); + fmpr_randtest(fmprb_radref(x), state, FMPRB_RAD_PREC, mag_bits); + fmpr_abs(fmprb_radref(x), fmprb_radref(x)); +} + +void +fmprb_randtest_precise(fmprb_t x, flint_rand_t state, long prec, long mag_bits) { fmpr_randtest(fmprb_midref(x), state, prec, mag_bits); - switch (n_randint(state, 8)) + if (fmpr_is_zero(fmprb_midref(x)) || (n_randint(state, 8) == 0)) { - /* exact */ - case 0: - fmpr_zero(fmprb_radref(x)); - break; - /* arbitrary radius */ - case 1: - fmpr_randtest(fmprb_radref(x), state, FMPRB_RAD_PREC, mag_bits); - fmpr_abs(fmprb_radref(x), fmprb_radref(x)); - break; - /* "typical" radius */ - default: - fmpr_randtest_not_zero(fmprb_radref(x), state, FMPRB_RAD_PREC, 4); - fmpr_abs(fmprb_radref(x), fmprb_radref(x)); - fmpz_add(fmpr_expref(fmprb_radref(x)), fmpr_expref(fmprb_radref(x)), - fmpr_expref(fmprb_midref(x))); + fmpr_zero(fmprb_radref(x)); + } + else + { + fmpz_randtest_not_zero(fmpr_manref(fmprb_radref(x)), state, FMPRB_RAD_PREC); + fmpz_abs(fmpr_manref(fmprb_radref(x)), fmpr_manref(fmprb_radref(x))); + fmpz_sub_ui(fmpr_expref(fmprb_radref(x)), fmpr_expref(fmprb_midref(x)), prec + FMPRB_RAD_PREC); + } +} + +void +fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits) +{ + switch (n_randint(state, 8)) + { + case 0: + fmprb_randtest_exact(x, state, prec, mag_bits); + break; + case 1: + fmprb_randtest_wide(x, state, prec, mag_bits); + break; + default: + fmprb_randtest_precise(x, state, prec, mag_bits); } } diff --git a/fmprb/test/t-gamma.c b/fmprb/test/t-gamma.c new file mode 100644 index 00000000..848e3996 --- /dev/null +++ b/fmprb/test/t-gamma.c @@ -0,0 +1,81 @@ +/*============================================================================= + + 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" + +int main() +{ + long iter; + flint_rand_t state; + + printf("gamma...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + fmprb_t a, b, c; + long prec1, prec2; + + prec1 = 2 + n_randint(state, 1000); + prec2 = prec1 + 30; + + fmprb_init(a); + fmprb_init(b); + fmprb_init(c); + + fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 3); + + fmprb_gamma(b, a, prec1); + fmprb_gamma(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_gamma(a, a, prec2); + + if (!fmprb_equal(a, c)) + { + printf("FAIL: aliasing\n\n"); + abort(); + } + + fmprb_clear(a); + fmprb_clear(b); + fmprb_clear(c); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmprb/test/t-lgamma.c b/fmprb/test/t-lgamma.c new file mode 100644 index 00000000..79da51ca --- /dev/null +++ b/fmprb/test/t-lgamma.c @@ -0,0 +1,81 @@ +/*============================================================================= + + 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" + +int main() +{ + long iter; + flint_rand_t state; + + printf("lgamma...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + fmprb_t a, b, c; + long prec1, prec2; + + prec1 = 2 + n_randint(state, 1000); + prec2 = prec1 + 30; + + fmprb_init(a); + fmprb_init(b); + fmprb_init(c); + + fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 3); + + fmprb_lgamma(b, a, prec1); + fmprb_lgamma(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_lgamma(a, a, prec2); + + if (!fmprb_equal(a, c)) + { + printf("FAIL: aliasing\n\n"); + abort(); + } + + fmprb_clear(a); + fmprb_clear(b); + fmprb_clear(c); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmprb/test/t-rgamma.c b/fmprb/test/t-rgamma.c new file mode 100644 index 00000000..2cf87a67 --- /dev/null +++ b/fmprb/test/t-rgamma.c @@ -0,0 +1,81 @@ +/*============================================================================= + + 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" + +int main() +{ + long iter; + flint_rand_t state; + + printf("rgamma...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + fmprb_t a, b, c; + long prec1, prec2; + + prec1 = 2 + n_randint(state, 1000); + prec2 = prec1 + 30; + + fmprb_init(a); + fmprb_init(b); + fmprb_init(c); + + fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 3); + + fmprb_rgamma(b, a, prec1); + fmprb_rgamma(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_rgamma(a, a, prec2); + + if (!fmprb_equal(a, c)) + { + printf("FAIL: aliasing\n\n"); + abort(); + } + + fmprb_clear(a); + fmprb_clear(b); + fmprb_clear(c); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}