diff options
Diffstat (limited to 'gamma.c')
-rw-r--r-- | gamma.c | 208 |
1 files changed, 101 insertions, 107 deletions
@@ -1,6 +1,6 @@ /* mpfr_gamma -- gamma function -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library, and was contributed by Mathieu Dutour. @@ -102,123 +102,117 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) while (!good) { - /* Precision stuff */ - if (compared < 0) - { - prec_nec = 2+realprec; /* We will use the reflexion formula! */ - } - else - { - prec_nec = realprec; - } - C = (double)(((double) prec_nec)*CST-0.5); - A = (int) C; - N = A - 1; + /* Precision stuff */ + prec_nec = compared < 0 ? + 2 + realprec /* We will use the reflexion formula! */ + : realprec; + C = (double)(((double) prec_nec)*CST-0.5); + A = (int) C; + N = A - 1; #ifdef DEBUG - printf("C=%u", (int)C); - printf(" A=%u", (int)A); - printf(" N=%u", (int)N); - printf("\n"); + printf("C=%u", (int)C); + printf(" A=%u", (int)A); + printf(" N=%u", (int)N); + printf("\n"); #endif - /* estimated_cancel is the amount of bit that will be flushed */ - estimated_cancel= (long) (ecCST * (double) A + 1.0); - Prec = prec_nec + estimated_cancel + 20; - - - mpfr_set_prec (xp, Prec); - if (compared < 0) - { - mpfr_ui_sub (xp, 1, x, GMP_RNDN); - } - else - { - mpfr_sub_ui (xp, x, 1, GMP_RNDN); - } - - /* Initialisation */ - mpfr_init2(tmp, Prec); - mpfr_init2(tmp2, Prec); - mpfr_init2(the_pi, Prec); - mpfr_init2(product, Prec); - mpfr_init2(GammaTrial, Prec); - - - mpfr_set_ui(GammaTrial, 0, GMP_RNDN); - sign=1; - for (k = 1; k<=N; k++) - { - mpfr_set_ui(tmp, A-k, GMP_RNDN); - mpfr_exp(product, tmp, GMP_RNDN); - mpfr_ui_pow_ui(tmp, A-k, k-1, GMP_RNDN); - mpfr_mul(product, product, tmp, GMP_RNDN); - mpfr_sqrt_ui(tmp, A-k, GMP_RNDN); - mpfr_mul(product, product, tmp, GMP_RNDN); - mpfr_fac_ui(tmp, k-1, GMP_RNDN); - mpfr_div(product, product, tmp, GMP_RNDN); - mpfr_add_ui(tmp, xp, k, GMP_RNDN); - mpfr_div(product, product, tmp, GMP_RNDN); - sign=-sign; - if (sign == 1) - { - mpfr_neg(product, product, GMP_RNDN); + /* estimated_cancel is the amount of bit that will be flushed */ + estimated_cancel= (long) (ecCST * (double) A + 1.0); + Prec = prec_nec + estimated_cancel + 20; + + mpfr_set_prec (xp, Prec); + if (compared < 0) + { + mpfr_ui_sub (xp, 1, x, GMP_RNDN); + } + else + { + mpfr_sub_ui (xp, x, 1, GMP_RNDN); + } + + /* Initialisation */ + mpfr_init2(tmp, Prec); + mpfr_init2(tmp2, Prec); + mpfr_init2(the_pi, Prec); + mpfr_init2(product, Prec); + mpfr_init2(GammaTrial, Prec); + + mpfr_set_ui(GammaTrial, 0, GMP_RNDN); + sign = 1; + for (k = 1; k <= N; k++) + { + mpfr_set_ui(tmp, A-k, GMP_RNDN); + mpfr_exp(product, tmp, GMP_RNDN); + mpfr_ui_pow_ui(tmp, A-k, k-1, GMP_RNDN); + mpfr_mul(product, product, tmp, GMP_RNDN); + mpfr_sqrt_ui(tmp, A-k, GMP_RNDN); + mpfr_mul(product, product, tmp, GMP_RNDN); + mpfr_fac_ui(tmp, k-1, GMP_RNDN); + mpfr_div(product, product, tmp, GMP_RNDN); + mpfr_add_ui(tmp, xp, k, GMP_RNDN); + mpfr_div(product, product, tmp, GMP_RNDN); + sign = -sign; + if (sign == 1) + { + mpfr_neg(product, product, GMP_RNDN); #ifdef DEBUG - /* printf(" k=%u", k); - printf("\n");*/ + /* printf(" k=%u", k); + printf("\n");*/ #endif - } - mpfr_add(GammaTrial, GammaTrial, product, GMP_RNDN); - } + } + mpfr_add(GammaTrial, GammaTrial, product, GMP_RNDN); + } #ifdef DEBUG - printf("GammaTrial ="); - mpfr_out_str (stdout, 10, 0, GammaTrial, GMP_RNDD); - printf ("\n"); + printf("GammaTrial ="); + mpfr_out_str (stdout, 10, 0, GammaTrial, GMP_RNDD); + printf ("\n"); #endif - mpfr_const_pi(the_pi, GMP_RNDN); - mpfr_const_pi(tmp, GMP_RNDN); - mpfr_mul_2ui(tmp, tmp, 1, GMP_RNDN); - mpfr_sqrt(tmp, tmp, GMP_RNDN); - mpfr_add(GammaTrial, GammaTrial, tmp, GMP_RNDN); - mpfr_add_ui(tmp2, xp, A, GMP_RNDN); - mpfr_set_ui(tmp, 1, GMP_RNDN); - mpfr_div_2ui(tmp, tmp, 1, GMP_RNDN); - mpfr_add(tmp, tmp, xp, GMP_RNDN); - mpfr_pow(tmp, tmp2, tmp, GMP_RNDN); - mpfr_mul(GammaTrial, GammaTrial, tmp, GMP_RNDN); - mpfr_neg(tmp, tmp2, GMP_RNDN); - mpfr_exp(tmp, tmp, GMP_RNDN); - mpfr_mul(GammaTrial, GammaTrial, tmp, GMP_RNDN); - if (compared < 0) - { - mpfr_sub_ui (tmp, x, 1, GMP_RNDN); - mpfr_mul (tmp, the_pi, tmp, GMP_RNDN); - mpfr_div (GammaTrial, tmp, GammaTrial, GMP_RNDN); - mpfr_sin (tmp, tmp, GMP_RNDN); - mpfr_div (GammaTrial, GammaTrial, tmp, GMP_RNDN); - } + mpfr_const_pi(the_pi, GMP_RNDN); + mpfr_const_pi(tmp, GMP_RNDN); + mpfr_mul_2ui(tmp, tmp, 1, GMP_RNDN); + mpfr_sqrt(tmp, tmp, GMP_RNDN); + mpfr_add(GammaTrial, GammaTrial, tmp, GMP_RNDN); + mpfr_add_ui(tmp2, xp, A, GMP_RNDN); + mpfr_set_ui(tmp, 1, GMP_RNDN); + mpfr_div_2ui(tmp, tmp, 1, GMP_RNDN); + mpfr_add(tmp, tmp, xp, GMP_RNDN); + mpfr_pow(tmp, tmp2, tmp, GMP_RNDN); + mpfr_mul(GammaTrial, GammaTrial, tmp, GMP_RNDN); + mpfr_neg(tmp, tmp2, GMP_RNDN); + mpfr_exp(tmp, tmp, GMP_RNDN); + mpfr_mul(GammaTrial, GammaTrial, tmp, GMP_RNDN); + if (compared < 0) + { + mpfr_sub_ui (tmp, x, 1, GMP_RNDN); + mpfr_mul (tmp, the_pi, tmp, GMP_RNDN); + mpfr_div (GammaTrial, tmp, GammaTrial, GMP_RNDN); + mpfr_sin (tmp, tmp, GMP_RNDN); + mpfr_div (GammaTrial, GammaTrial, tmp, GMP_RNDN); + } #ifdef DEBUG - printf("GammaTrial ="); - mpfr_out_str (stdout, 10, 0, GammaTrial, GMP_RNDD); - printf ("\n"); + printf("GammaTrial ="); + mpfr_out_str (stdout, 10, 0, GammaTrial, GMP_RNDD); + printf ("\n"); #endif - if (mpfr_can_round (GammaTrial, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(gamma))) - { - mpfr_set (gamma, GammaTrial, rnd_mode); - good = 1; - } - else - { - realprec += __gmpfr_ceil_log2 ((double) realprec); + if (mpfr_can_round (GammaTrial, realprec, GMP_RNDD, rnd_mode, + MPFR_PREC(gamma))) + { + mpfr_set (gamma, GammaTrial, rnd_mode); + good = 1; + } + else + { + realprec += __gmpfr_ceil_log2 ((double) realprec); #ifdef DEBUG - printf("RETRY\n"); + printf("RETRY\n"); #endif - } - mpfr_clear(tmp); - mpfr_clear(tmp2); - mpfr_clear(the_pi); - mpfr_clear(product); - mpfr_clear(GammaTrial); - } + } + mpfr_clear(tmp); + mpfr_clear(tmp2); + mpfr_clear(the_pi); + mpfr_clear(product); + mpfr_clear(GammaTrial); + } mpfr_clear (xp); |