diff options
Diffstat (limited to 'gsl-1.9/specfunc/exp.c')
-rw-r--r-- | gsl-1.9/specfunc/exp.c | 608 |
1 files changed, 608 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/exp.c b/gsl-1.9/specfunc/exp.c new file mode 100644 index 0000000..02030d0 --- /dev/null +++ b/gsl-1.9/specfunc/exp.c @@ -0,0 +1,608 @@ +/* specfunc/exp.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * This program 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. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sf_gamma.h> +#include <gsl/gsl_sf_exp.h> + +#include "error.h" + +/* Evaluate the continued fraction for exprel. + * [Abramowitz+Stegun, 4.2.41] + */ +static +int +exprel_n_CF(const int N, const double x, gsl_sf_result * result) +{ + const double RECUR_BIG = GSL_SQRT_DBL_MAX; + const int maxiter = 5000; + int n = 1; + double Anm2 = 1.0; + double Bnm2 = 0.0; + double Anm1 = 0.0; + double Bnm1 = 1.0; + double a1 = 1.0; + double b1 = 1.0; + double a2 = -x; + double b2 = N+1; + double an, bn; + + double fn; + + double An = b1*Anm1 + a1*Anm2; /* A1 */ + double Bn = b1*Bnm1 + a1*Bnm2; /* B1 */ + + /* One explicit step, before we get to the main pattern. */ + n++; + Anm2 = Anm1; + Bnm2 = Bnm1; + Anm1 = An; + Bnm1 = Bn; + An = b2*Anm1 + a2*Anm2; /* A2 */ + Bn = b2*Bnm1 + a2*Bnm2; /* B2 */ + + fn = An/Bn; + + while(n < maxiter) { + double old_fn; + double del; + n++; + Anm2 = Anm1; + Bnm2 = Bnm1; + Anm1 = An; + Bnm1 = Bn; + an = ( GSL_IS_ODD(n) ? ((n-1)/2)*x : -(N+(n/2)-1)*x ); + bn = N + n - 1; + An = bn*Anm1 + an*Anm2; + Bn = bn*Bnm1 + an*Bnm2; + + if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) { + An /= RECUR_BIG; + Bn /= RECUR_BIG; + Anm1 /= RECUR_BIG; + Bnm1 /= RECUR_BIG; + Anm2 /= RECUR_BIG; + Bnm2 /= RECUR_BIG; + } + + old_fn = fn; + fn = An/Bn; + del = old_fn/fn; + + if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break; + } + + result->val = fn; + result->err = 2.0*(n+1.0)*GSL_DBL_EPSILON*fabs(fn); + + if(n == maxiter) + GSL_ERROR ("error", GSL_EMAXITER); + else + return GSL_SUCCESS; +} + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int gsl_sf_exp_e(const double x, gsl_sf_result * result) +{ + if(x > GSL_LOG_DBL_MAX) { + OVERFLOW_ERROR(result); + } + else if(x < GSL_LOG_DBL_MIN) { + UNDERFLOW_ERROR(result); + } + else { + result->val = exp(x); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } +} + +int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result) +{ + if(x > INT_MAX-1) { + OVERFLOW_ERROR_E10(result); + } + else if(x < INT_MIN+1) { + UNDERFLOW_ERROR_E10(result); + } + else { + const int N = (int) floor(x/M_LN10); + result->val = exp(x-N*M_LN10); + result->err = 2.0 * (fabs(x)+1.0) * GSL_DBL_EPSILON * fabs(result->val); + result->e10 = N; + return GSL_SUCCESS; + } +} + + +int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result) +{ + const double ay = fabs(y); + + if(y == 0.0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN) + && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN) + ) { + const double ex = exp(x); + result->val = y * ex; + result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + const double ly = log(ay); + const double lnr = x + ly; + + if(lnr > GSL_LOG_DBL_MAX - 0.01) { + OVERFLOW_ERROR(result); + } + else if(lnr < GSL_LOG_DBL_MIN + 0.01) { + UNDERFLOW_ERROR(result); + } + else { + const double sy = GSL_SIGN(y); + const double M = floor(x); + const double N = floor(ly); + const double a = x - M; + const double b = ly - N; + const double berr = 2.0 * GSL_DBL_EPSILON * (fabs(ly) + fabs(N)); + result->val = sy * exp(M+N) * exp(a+b); + result->err = berr * fabs(result->val); + result->err += 2.0 * GSL_DBL_EPSILON * (M + N + 1.0) * fabs(result->val); + return GSL_SUCCESS; + } + } +} + + +int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result) +{ + const double ay = fabs(y); + + if(y == 0.0) { + result->val = 0.0; + result->err = 0.0; + result->e10 = 0; + return GSL_SUCCESS; + } + else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN) + && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN) + ) { + const double ex = exp(x); + result->val = y * ex; + result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val); + result->e10 = 0; + return GSL_SUCCESS; + } + else { + const double ly = log(ay); + const double l10_val = (x + ly)/M_LN10; + + if(l10_val > INT_MAX-1) { + OVERFLOW_ERROR_E10(result); + } + else if(l10_val < INT_MIN+1) { + UNDERFLOW_ERROR_E10(result); + } + else { + const double sy = GSL_SIGN(y); + const int N = (int) floor(l10_val); + const double arg_val = (l10_val - N) * M_LN10; + const double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(ly); + + result->val = sy * exp(arg_val); + result->err = arg_err * fabs(result->val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + result->e10 = N; + + return GSL_SUCCESS; + } + } +} + + +int gsl_sf_exp_mult_err_e(const double x, const double dx, + const double y, const double dy, + gsl_sf_result * result) +{ + const double ay = fabs(y); + + if(y == 0.0) { + result->val = 0.0; + result->err = fabs(dy * exp(x)); + return GSL_SUCCESS; + } + else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN) + && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN) + ) { + double ex = exp(x); + result->val = y * ex; + result->err = ex * (fabs(dy) + fabs(y*dx)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + const double ly = log(ay); + const double lnr = x + ly; + + if(lnr > GSL_LOG_DBL_MAX - 0.01) { + OVERFLOW_ERROR(result); + } + else if(lnr < GSL_LOG_DBL_MIN + 0.01) { + UNDERFLOW_ERROR(result); + } + else { + const double sy = GSL_SIGN(y); + const double M = floor(x); + const double N = floor(ly); + const double a = x - M; + const double b = ly - N; + const double eMN = exp(M+N); + const double eab = exp(a+b); + result->val = sy * eMN * eab; + result->err = eMN * eab * 2.0*GSL_DBL_EPSILON; + result->err += eMN * eab * fabs(dy/y); + result->err += eMN * eab * fabs(dx); + return GSL_SUCCESS; + } + } +} + + +int gsl_sf_exp_mult_err_e10_e(const double x, const double dx, + const double y, const double dy, + gsl_sf_result_e10 * result) +{ + const double ay = fabs(y); + + if(y == 0.0) { + result->val = 0.0; + result->err = fabs(dy * exp(x)); + result->e10 = 0; + return GSL_SUCCESS; + } + else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN) + && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN) + ) { + const double ex = exp(x); + result->val = y * ex; + result->err = ex * (fabs(dy) + fabs(y*dx)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + result->e10 = 0; + return GSL_SUCCESS; + } + else { + const double ly = log(ay); + const double l10_val = (x + ly)/M_LN10; + + if(l10_val > INT_MAX-1) { + OVERFLOW_ERROR_E10(result); + } + else if(l10_val < INT_MIN+1) { + UNDERFLOW_ERROR_E10(result); + } + else { + const double sy = GSL_SIGN(y); + const int N = (int) floor(l10_val); + const double arg_val = (l10_val - N) * M_LN10; + const double arg_err = dy/fabs(y) + dx + 2.0*GSL_DBL_EPSILON*fabs(arg_val); + + result->val = sy * exp(arg_val); + result->err = arg_err * fabs(result->val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + result->e10 = N; + + return GSL_SUCCESS; + } + } +} + + +int gsl_sf_expm1_e(const double x, gsl_sf_result * result) +{ + const double cut = 0.002; + + if(x < GSL_LOG_DBL_MIN) { + result->val = -1.0; + result->err = GSL_DBL_EPSILON; + return GSL_SUCCESS; + } + else if(x < -cut) { + result->val = exp(x) - 1.0; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < cut) { + result->val = x * (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x)))); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < GSL_LOG_DBL_MAX) { + result->val = exp(x) - 1.0; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + OVERFLOW_ERROR(result); + } +} + + +int gsl_sf_exprel_e(const double x, gsl_sf_result * result) +{ + const double cut = 0.002; + + if(x < GSL_LOG_DBL_MIN) { + result->val = -1.0/x; + result->err = GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < -cut) { + result->val = (exp(x) - 1.0)/x; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < cut) { + result->val = (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x)))); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < GSL_LOG_DBL_MAX) { + result->val = (exp(x) - 1.0)/x; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + OVERFLOW_ERROR(result); + } +} + + +int gsl_sf_exprel_2_e(double x, gsl_sf_result * result) +{ + const double cut = 0.002; + + if(x < GSL_LOG_DBL_MIN) { + result->val = -2.0/x*(1.0 + 1.0/x); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < -cut) { + result->val = 2.0*(exp(x) - 1.0 - x)/(x*x); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < cut) { + result->val = (1.0 + 1.0/3.0*x*(1.0 + 0.25*x*(1.0 + 0.2*x*(1.0 + 1.0/6.0*x)))); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < GSL_LOG_DBL_MAX) { + result->val = 2.0*(exp(x) - 1.0 - x)/(x*x); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + OVERFLOW_ERROR(result); + } +} + + +int +gsl_sf_exprel_n_e(const int N, const double x, gsl_sf_result * result) +{ + if(N < 0) { + DOMAIN_ERROR(result); + } + else if(x == 0.0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(fabs(x) < GSL_ROOT3_DBL_EPSILON * N) { + result->val = 1.0 + x/(N+1) * (1.0 + x/(N+2)); + result->err = 2.0 * GSL_DBL_EPSILON; + return GSL_SUCCESS; + } + else if(N == 0) { + return gsl_sf_exp_e(x, result); + } + else if(N == 1) { + return gsl_sf_exprel_e(x, result); + } + else if(N == 2) { + return gsl_sf_exprel_2_e(x, result); + } + else { + if(x > N && (-x + N*(1.0 + log(x/N)) < GSL_LOG_DBL_EPSILON)) { + /* x is much larger than n. + * Ignore polynomial part, so + * exprel_N(x) ~= e^x N!/x^N + */ + gsl_sf_result lnf_N; + double lnr_val; + double lnr_err; + double lnterm; + gsl_sf_lnfact_e(N, &lnf_N); + lnterm = N*log(x); + lnr_val = x + lnf_N.val - lnterm; + lnr_err = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(lnterm)); + lnr_err += lnf_N.err; + return gsl_sf_exp_err_e(lnr_val, lnr_err, result); + } + else if(x > N) { + /* Write the identity + * exprel_n(x) = e^x n! / x^n (1 - Gamma[n,x]/Gamma[n]) + * then use the asymptotic expansion + * Gamma[n,x] ~ x^(n-1) e^(-x) (1 + (n-1)/x + (n-1)(n-2)/x^2 + ...) + */ + double ln_x = log(x); + gsl_sf_result lnf_N; + double lg_N; + double lnpre_val; + double lnpre_err; + gsl_sf_lnfact_e(N, &lnf_N); /* log(N!) */ + lg_N = lnf_N.val - log(N); /* log(Gamma(N)) */ + lnpre_val = x + lnf_N.val - N*ln_x; + lnpre_err = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(N*ln_x)); + lnpre_err += lnf_N.err; + if(lnpre_val < GSL_LOG_DBL_MAX - 5.0) { + int stat_eG; + gsl_sf_result bigG_ratio; + gsl_sf_result pre; + int stat_ex = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &pre); + double ln_bigG_ratio_pre = -x + (N-1)*ln_x - lg_N; + double bigGsum = 1.0; + double term = 1.0; + int k; + for(k=1; k<N; k++) { + term *= (N-k)/x; + bigGsum += term; + } + stat_eG = gsl_sf_exp_mult_e(ln_bigG_ratio_pre, bigGsum, &bigG_ratio); + if(stat_eG == GSL_SUCCESS) { + result->val = pre.val * (1.0 - bigG_ratio.val); + result->err = pre.val * (2.0*GSL_DBL_EPSILON + bigG_ratio.err); + result->err += pre.err * fabs(1.0 - bigG_ratio.val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_ex; + } + else { + result->val = 0.0; + result->err = 0.0; + return stat_eG; + } + } + else { + OVERFLOW_ERROR(result); + } + } + else if(x > -10.0*N) { + return exprel_n_CF(N, x, result); + } + else { + /* x -> -Inf asymptotic: + * exprel_n(x) ~ e^x n!/x^n - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...) + * ~ - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...) + */ + double sum = 1.0; + double term = 1.0; + int k; + for(k=1; k<N; k++) { + term *= (N-k)/x; + sum += term; + } + result->val = -N/x * sum; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + } +} + + +int +gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result) +{ + const double adx = fabs(dx); + + /* CHECK_POINTER(result) */ + + if(x + adx > GSL_LOG_DBL_MAX) { + OVERFLOW_ERROR(result); + } + else if(x - adx < GSL_LOG_DBL_MIN) { + UNDERFLOW_ERROR(result); + } + else { + const double ex = exp(x); + const double edx = exp(adx); + result->val = ex; + result->err = ex * GSL_MAX_DBL(GSL_DBL_EPSILON, edx - 1.0/edx); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } +} + + +int +gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result) +{ + const double adx = fabs(dx); + + /* CHECK_POINTER(result) */ + + if(x + adx > INT_MAX - 1) { + OVERFLOW_ERROR_E10(result); + } + else if(x - adx < INT_MIN + 1) { + UNDERFLOW_ERROR_E10(result); + } + else { + const int N = (int)floor(x/M_LN10); + const double ex = exp(x-N*M_LN10); + result->val = ex; + result->err = ex * (2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) + adx); + result->e10 = N; + return GSL_SUCCESS; + } +} + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_exp(const double x) +{ + EVAL_RESULT(gsl_sf_exp_e(x, &result)); +} + +double gsl_sf_exp_mult(const double x, const double y) +{ + EVAL_RESULT(gsl_sf_exp_mult_e(x, y, &result)); +} + +double gsl_sf_expm1(const double x) +{ + EVAL_RESULT(gsl_sf_expm1_e(x, &result)); +} + +double gsl_sf_exprel(const double x) +{ + EVAL_RESULT(gsl_sf_exprel_e(x, &result)); +} + +double gsl_sf_exprel_2(const double x) +{ + EVAL_RESULT(gsl_sf_exprel_2_e(x, &result)); +} + +double gsl_sf_exprel_n(const int n, const double x) +{ + EVAL_RESULT(gsl_sf_exprel_n_e(n, x, &result)); +} |