diff options
| -rw-r--r-- | lngamma.c | 119 | ||||
| -rw-r--r-- | tests/Makefile.am | 2 | ||||
| -rw-r--r-- | tests/tlgamma.c | 275 | ||||
| -rw-r--r-- | tests/tlngamma.c | 12 |
4 files changed, 377 insertions, 31 deletions
@@ -146,12 +146,12 @@ unit_bit (mpfr_srcptr (x)) For z real, |K(z)| <= 1 thus R_n(z) is bounded by the first neglected term. */ #ifdef IS_GAMMA -static int #define GAMMA_FUNC mpfr_gamma_aux #else -int -#define GAMMA_FUNC mpfr_lngamma +#define GAMMA_FUNC mpfr_lngamma_aux #endif + +static int GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd) { mp_prec_t precy, w; /* working precision */ @@ -165,31 +165,6 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd) double d; MPFR_SAVE_EXPO_DECL (expo); -#ifndef IS_GAMMA - /* special cases */ - if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z0))) - { - if (MPFR_IS_NAN (z0) || MPFR_IS_NEG (z0)) - { - MPFR_SET_NAN (y); - MPFR_RET_NAN; - } - else /* lngamma(+Inf) = lngamma(+0) = +Inf */ - { - MPFR_SET_INF (y); - MPFR_SET_POS (y); - MPFR_RET (0); /* exact */ - } - } - - /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */ - if (MPFR_IS_NEG (z0) && (unit_bit (z0) == 0 || mpfr_integer_p (z0))) - { - MPFR_SET_NAN (y); - MPFR_RET_NAN; - } -#endif - precy = MPFR_PREC(y); compared = mpfr_cmp_ui (z0, 1); @@ -251,8 +226,7 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd) mpfr_sub_ui (v, z0, 1, GMP_RNDN); /* v = (x-1) * (1+u) */ mpfr_mul (v, v, t, GMP_RNDN); /* v = Pi*(x-1) * (1+u)^3 */ mpfr_div (v, v, s, GMP_RNDN); /* Pi*(x-1)/sin(Pi*(2-x)) */ - if (MPFR_IS_NEG (v)) /* can be explained by a too large error? */ - continue; /* so, start again with a larger precision */ + mpfr_abs (v, v, GMP_RNDN); /* (1+u)^(4+2^err_s+1) */ err_s = (err_s <= 2) ? 3 + (err_s / 2) : err_s + 1; mpfr_log (v, v, GMP_RNDN); @@ -486,3 +460,88 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd) MPFR_SAVE_EXPO_FREE (expo); return mpfr_check_range (y, inexact, rnd); } + +#ifndef IS_GAMMA + +int +mpfr_lngamma (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd) +{ + int inex; + + MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd), + ("lngamma[%#R]=%R inexact=%d", y, y, inex)); + + /* special cases */ + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) + { + if (MPFR_IS_NAN (x) || MPFR_IS_NEG (x)) + { + MPFR_SET_NAN (y); + MPFR_RET_NAN; + } + else /* lngamma(+Inf) = lngamma(+0) = +Inf */ + { + MPFR_SET_INF (y); + MPFR_SET_POS (y); + MPFR_RET (0); /* exact */ + } + } + + /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */ + if (MPFR_IS_NEG (x) && (unit_bit (x) == 0 || mpfr_integer_p (x))) + { + MPFR_SET_NAN (y); + MPFR_RET_NAN; + } + + inex = mpfr_lngamma_aux (y, x, rnd); + return inex; +} + +int +mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mp_rnd_t rnd) +{ + int inex; + + MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd), + ("lgamma[%#R]=%R inexact=%d", y, y, inex)); + + *signp = 1; /* most common case */ + + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) + { + if (MPFR_IS_NAN (x)) + { + MPFR_SET_NAN (y); + MPFR_RET_NAN; + } + else + { + *signp = MPFR_INT_SIGN (x); + MPFR_SET_INF (y); + MPFR_SET_POS (y); + MPFR_RET (0); + } + } + + if (MPFR_IS_NEG (x)) + { + if (mpfr_integer_p (x)) + { + MPFR_SET_INF (y); + MPFR_SET_POS (y); + MPFR_RET (0); + } + + if (unit_bit (x) == 0) + *signp = -1; + } + + /* TODO: add support for negative numbers with small exponent. + Re-enable the generic tests when this is done. */ + + inex = mpfr_lngamma_aux (y, x, rnd); + return inex; +} + +#endif diff --git a/tests/Makefile.am b/tests/Makefile.am index 381a42daf..dee71e676 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,6 +1,6 @@ AUTOMAKE_OPTIONS = 1.6 gnu $(top_builddir)/ansi2knr -check_PROGRAMS = tversion tinternals tinits tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tgmpop tsi_op tmul_2exp tfma tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal tlngamma tzeta_ui tget_ld_2exp tget_set_d64 tj0 tj1 tjn ty0 ty1 tyn +check_PROGRAMS = tversion tinternals tinits tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tgmpop tsi_op tmul_2exp tfma tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal tlngamma tlgamma tzeta_ui tget_ld_2exp tget_set_d64 tj0 tj1 tjn ty0 ty1 tyn EXTRA_DIST = tgeneric.c tgeneric_ui.c mpf_compat.h inp_str.data diff --git a/tests/tlgamma.c b/tests/tlgamma.c new file mode 100644 index 000000000..4a270b4ee --- /dev/null +++ b/tests/tlgamma.c @@ -0,0 +1,275 @@ +/* mpfr_tlgamma -- test file for lgamma function + +Copyright 2005, 2006, 2007 Free Software Foundation, Inc. +Contributed by the Arenaire and Cacao projects, INRIA. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library 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 Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, +MA 02110-1301, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +static int +mpfr_lgamma_nosign (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd) +{ + int inex, sign; + + inex = mpfr_lgamma (y, &sign, x, rnd); + if (!MPFR_IS_SINGULAR (y)) + { + MPFR_ASSERTN (sign == 1 || sign == -1); + if (sign == -1 && (rnd == GMP_RNDN || rnd == GMP_RNDZ)) + mpfr_neg (y, y, GMP_RNDN); + /* This is a way to check that the sign is consistent. */ + } + + return inex; +} + +#define TEST_FUNCTION mpfr_lgamma_nosign +#include "tgeneric.c" + +static void +special (void) +{ + mpfr_t x, y; + int inex; + int sign; + + mpfr_init (x); + mpfr_init (y); + + mpfr_set_nan (x); + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error for lgamma(NaN)\n"); + exit (1); + } + + mpfr_set_inf (x, -1); + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) + { + printf ("Error for lgamma(-Inf)\n"); + exit (1); + } + + mpfr_set_inf (x, 1); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1) + { + printf ("Error for lgamma(+Inf)\n"); + exit (1); + } + + mpfr_set_ui (x, 0, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1) + { + printf ("Error for lgamma(+0)\n"); + exit (1); + } + + mpfr_set_ui (x, 0, GMP_RNDN); + mpfr_neg (x, x, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != -1) + { + printf ("Error for lgamma(-0)\n"); + exit (1); + } + + mpfr_set_ui (x, 1, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1) + { + printf ("Error for lgamma(1)\n"); + exit (1); + } + + mpfr_set_si (x, -1, GMP_RNDN); + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) + { + printf ("Error for lgamma(-1)\n"); + exit (1); + } + + mpfr_set_ui (x, 2, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1) + { + printf ("Error for lgamma(2)\n"); + exit (1); + } + + mpfr_set_prec (x, 53); + mpfr_set_prec (y, 53); + +#define CHECK_X1 "1.0762904832837976166" +#define CHECK_Y1 "-0.039418362817587634939" + + mpfr_set_str (x, CHECK_X1, 10, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + mpfr_set_str (x, CHECK_Y1, 10, GMP_RNDN); + if (MPFR_IS_NAN (y) || mpfr_cmp (y, x) || sign != 1) + { + printf ("mpfr_lgamma("CHECK_X1") is wrong:\n" + "expected "); + mpfr_print_binary (x); putchar ('\n'); + printf ("got "); + mpfr_print_binary (y); putchar ('\n'); + exit (1); + } + +#define CHECK_X2 "9.23709516716202383435e-01" +#define CHECK_Y2 "0.049010669407893718563" + mpfr_set_str (x, CHECK_X2, 10, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + mpfr_set_str (x, CHECK_Y2, 10, GMP_RNDN); + if (MPFR_IS_NAN (y) || mpfr_cmp (y, x) || sign != 1) + { + printf ("mpfr_lgamma("CHECK_X2") is wrong:\n" + "expected "); + mpfr_print_binary (x); putchar ('\n'); + printf ("got "); + mpfr_print_binary (y); putchar ('\n'); + exit (1); + } + + mpfr_set_prec (x, 8); + mpfr_set_prec (y, 175); + mpfr_set_ui (x, 33, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDU); + mpfr_set_prec (x, 175); + mpfr_set_str_binary (x, "0.1010001100011101101011001101110010100001000001000001110011000001101100001111001001000101011011100100010101011110100111110101010100010011010010000101010111001100011000101111E7"); + if (MPFR_IS_NAN (y) || mpfr_cmp (x, y) || sign != 1) + { + printf ("Error in mpfr_lgamma (1)\n"); + exit (1); + } + + mpfr_set_prec (x, 21); + mpfr_set_prec (y, 8); + mpfr_set_ui (y, 120, GMP_RNDN); + sign = -17; + mpfr_lgamma (x, &sign, y, GMP_RNDZ); + mpfr_set_prec (y, 21); + mpfr_set_str_binary (y, "0.111000101000001100101E9"); + if (MPFR_IS_NAN (x) || mpfr_cmp (x, y) || sign != 1) + { + printf ("Error in mpfr_lgamma (120)\n"); + printf ("Expected "); mpfr_print_binary (y); puts (""); + printf ("Got "); mpfr_print_binary (x); puts (""); + exit (1); + } + + mpfr_set_prec (x, 3); + mpfr_set_prec (y, 206); + mpfr_set_str_binary (x, "0.110e10"); + sign = -17; + inex = mpfr_lgamma (y, &sign, x, GMP_RNDN); + mpfr_set_prec (x, 206); + mpfr_set_str_binary (x, "0.10000111011000000011100010101001100110001110000111100011000100100110110010001011011110101001111011110110000001010100111011010000000011100110110101100111000111010011110010000100010111101010001101000110101001E13"); + if (MPFR_IS_NAN (y) || mpfr_cmp (x, y) || sign != 1) + { + printf ("Error in mpfr_lgamma (768)\n"); + exit (1); + } + if (inex >= 0) + { + printf ("Wrong flag for mpfr_lgamma (768)\n"); + exit (1); + } + + mpfr_set_prec (x, 4); + mpfr_set_prec (y, 4); + mpfr_set_str_binary (x, "0.1100E-66"); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + mpfr_set_str_binary (x, "0.1100E6"); + if (MPFR_IS_NAN (y) || mpfr_cmp (x, y) || sign != 1) + { + printf ("Error for lgamma(0.1100E-66)\n"); + exit (1); + } + + mpfr_set_prec (x, 256); + mpfr_set_prec (y, 32); + mpfr_set_si_2exp (x, -1, 200, GMP_RNDN); + mpfr_add_ui (x, x, 1, GMP_RNDN); + mpfr_div_2ui (x, x, 1, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + mpfr_set_prec (x, 32); + mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207"); + if (MPFR_IS_NAN (y) || mpfr_cmp (x, y) || sign != 1) + { + printf ("Error for lgamma(-2^199+0.5)\n"); + printf ("Got "); + mpfr_dump (y); + printf ("instead of "); + mpfr_dump (x); + exit (1); + } + + mpfr_set_prec (x, 256); + mpfr_set_prec (y, 32); + mpfr_set_si_2exp (x, -1, 200, GMP_RNDN); + mpfr_sub_ui (x, x, 1, GMP_RNDN); + mpfr_div_2ui (x, x, 1, GMP_RNDN); + sign = -17; + mpfr_lgamma (y, &sign, x, GMP_RNDN); + mpfr_set_prec (x, 32); + mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207"); + if (MPFR_IS_NAN (y) || mpfr_cmp (x, y) || sign != -1) + { + printf ("Error for lgamma(-2^199-0.5)\n"); + printf ("Got "); + mpfr_dump (y); + printf ("with sign %d instead of ", sign); + mpfr_dump (x); + printf ("with sign -1.\n"); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); +} + +int +main (void) +{ + tests_start_mpfr (); + + special (); + /* test_generic (2, 100, 2); */ + + tests_end_mpfr (); + return 0; +} diff --git a/tests/tlngamma.c b/tests/tlngamma.c index 0c7f81351..30dc9d848 100644 --- a/tests/tlngamma.c +++ b/tests/tlngamma.c @@ -208,6 +208,18 @@ special (void) exit (1); } + mpfr_set_prec (x, 256); + mpfr_set_prec (y, 32); + mpfr_set_si_2exp (x, -1, 200, GMP_RNDN); + mpfr_sub_ui (x, x, 1, GMP_RNDN); + mpfr_div_2ui (x, x, 1, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error for lngamma(-2^199-0.5)\n"); + exit (1); + } + mpfr_clear (x); mpfr_clear (y); } |