|
Packit |
67cb25 |
/* specfunc/exp.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_gamma.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_exp.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate the continued fraction for exprel.
|
|
Packit |
67cb25 |
* [Abramowitz+Stegun, 4.2.41]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
exprel_n_CF(const double N, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double RECUR_BIG = GSL_SQRT_DBL_MAX;
|
|
Packit |
67cb25 |
const int maxiter = 5000;
|
|
Packit |
67cb25 |
int n = 1;
|
|
Packit |
67cb25 |
double Anm2 = 1.0;
|
|
Packit |
67cb25 |
double Bnm2 = 0.0;
|
|
Packit |
67cb25 |
double Anm1 = 0.0;
|
|
Packit |
67cb25 |
double Bnm1 = 1.0;
|
|
Packit |
67cb25 |
double a1 = 1.0;
|
|
Packit |
67cb25 |
double b1 = 1.0;
|
|
Packit |
67cb25 |
double a2 = -x;
|
|
Packit |
67cb25 |
double b2 = N+1;
|
|
Packit |
67cb25 |
double an, bn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double fn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double An = b1*Anm1 + a1*Anm2; /* A1 */
|
|
Packit |
67cb25 |
double Bn = b1*Bnm1 + a1*Bnm2; /* B1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* One explicit step, before we get to the main pattern. */
|
|
Packit |
67cb25 |
n++;
|
|
Packit |
67cb25 |
Anm2 = Anm1;
|
|
Packit |
67cb25 |
Bnm2 = Bnm1;
|
|
Packit |
67cb25 |
Anm1 = An;
|
|
Packit |
67cb25 |
Bnm1 = Bn;
|
|
Packit |
67cb25 |
An = b2*Anm1 + a2*Anm2; /* A2 */
|
|
Packit |
67cb25 |
Bn = b2*Bnm1 + a2*Bnm2; /* B2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn = An/Bn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while(n < maxiter) {
|
|
Packit |
67cb25 |
double old_fn;
|
|
Packit |
67cb25 |
double del;
|
|
Packit |
67cb25 |
n++;
|
|
Packit |
67cb25 |
Anm2 = Anm1;
|
|
Packit |
67cb25 |
Bnm2 = Bnm1;
|
|
Packit |
67cb25 |
Anm1 = An;
|
|
Packit |
67cb25 |
Bnm1 = Bn;
|
|
Packit |
67cb25 |
an = ( GSL_IS_ODD(n) ? ((n-1)/2)*x : -(N+(n/2)-1)*x );
|
|
Packit |
67cb25 |
bn = N + n - 1;
|
|
Packit |
67cb25 |
An = bn*Anm1 + an*Anm2;
|
|
Packit |
67cb25 |
Bn = bn*Bnm1 + an*Bnm2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
|
|
Packit |
67cb25 |
An /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bn /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
old_fn = fn;
|
|
Packit |
67cb25 |
fn = An/Bn;
|
|
Packit |
67cb25 |
del = old_fn/fn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = fn;
|
|
Packit |
67cb25 |
result->err = 4.0*(n+1.0)*GSL_DBL_EPSILON*fabs(fn);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(n == maxiter)
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exp_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x > GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < GSL_LOG_DBL_MIN) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = exp(x);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x > INT_MAX-1) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < INT_MIN+1) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const int N = (x > GSL_LOG_DBL_MAX || x < GSL_LOG_DBL_MIN) ? (int) floor(x/M_LN10) : 0;
|
|
Packit |
67cb25 |
result->val = exp(x-N*M_LN10);
|
|
Packit |
67cb25 |
result->err = 2.0 * (fabs(x)+1.0) * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
result->e10 = N;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ay = fabs(y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(y == 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
|
|
Packit |
67cb25 |
&& (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
const double ex = exp(x);
|
|
Packit |
67cb25 |
result->val = y * ex;
|
|
Packit |
67cb25 |
result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double ly = log(ay);
|
|
Packit |
67cb25 |
const double lnr = x + ly;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(lnr > GSL_LOG_DBL_MAX - 0.01) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double sy = GSL_SIGN(y);
|
|
Packit |
67cb25 |
const double M = floor(x);
|
|
Packit |
67cb25 |
const double N = floor(ly);
|
|
Packit |
67cb25 |
const double a = x - M;
|
|
Packit |
67cb25 |
const double b = ly - N;
|
|
Packit |
67cb25 |
const double berr = 2.0 * GSL_DBL_EPSILON * (fabs(ly) + fabs(N));
|
|
Packit |
67cb25 |
result->val = sy * exp(M+N) * exp(a+b);
|
|
Packit |
67cb25 |
result->err = berr * fabs(result->val);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * (M + N + 1.0) * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ay = fabs(y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(y == 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
result->e10 = 0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
|
|
Packit |
67cb25 |
&& (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
const double ex = exp(x);
|
|
Packit |
67cb25 |
result->val = y * ex;
|
|
Packit |
67cb25 |
result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
result->e10 = 0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double ly = log(ay);
|
|
Packit |
67cb25 |
const double l10_val = (x + ly)/M_LN10;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(l10_val > INT_MAX-1) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l10_val < INT_MIN+1) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double sy = GSL_SIGN(y);
|
|
Packit |
67cb25 |
const int N = (int) floor(l10_val);
|
|
Packit |
67cb25 |
const double arg_val = (l10_val - N) * M_LN10;
|
|
Packit |
67cb25 |
const double arg_err = 2.0 * GSL_DBL_EPSILON * (fabs(x) + fabs(ly) + M_LN10*fabs(N));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = sy * exp(arg_val);
|
|
Packit |
67cb25 |
result->err = arg_err * fabs(result->val);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
result->e10 = N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exp_mult_err_e(const double x, const double dx,
|
|
Packit |
67cb25 |
const double y, const double dy,
|
|
Packit |
67cb25 |
gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ay = fabs(y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(y == 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = fabs(dy * exp(x));
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
|
|
Packit |
67cb25 |
&& (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
double ex = exp(x);
|
|
Packit |
67cb25 |
result->val = y * ex;
|
|
Packit |
67cb25 |
result->err = ex * (fabs(dy) + fabs(y*dx));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double ly = log(ay);
|
|
Packit |
67cb25 |
const double lnr = x + ly;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(lnr > GSL_LOG_DBL_MAX - 0.01) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double sy = GSL_SIGN(y);
|
|
Packit |
67cb25 |
const double M = floor(x);
|
|
Packit |
67cb25 |
const double N = floor(ly);
|
|
Packit |
67cb25 |
const double a = x - M;
|
|
Packit |
67cb25 |
const double b = ly - N;
|
|
Packit |
67cb25 |
const double eMN = exp(M+N);
|
|
Packit |
67cb25 |
const double eab = exp(a+b);
|
|
Packit |
67cb25 |
result->val = sy * eMN * eab;
|
|
Packit |
67cb25 |
result->err = eMN * eab * 2.0*GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
result->err += eMN * eab * fabs(dy/y);
|
|
Packit |
67cb25 |
result->err += eMN * eab * fabs(dx);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exp_mult_err_e10_e(const double x, const double dx,
|
|
Packit |
67cb25 |
const double y, const double dy,
|
|
Packit |
67cb25 |
gsl_sf_result_e10 * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ay = fabs(y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(y == 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = fabs(dy * exp(x));
|
|
Packit |
67cb25 |
result->e10 = 0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
|
|
Packit |
67cb25 |
&& (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
const double ex = exp(x);
|
|
Packit |
67cb25 |
result->val = y * ex;
|
|
Packit |
67cb25 |
result->err = ex * (fabs(dy) + fabs(y*dx));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
result->e10 = 0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double ly = log(ay);
|
|
Packit |
67cb25 |
const double l10_val = (x + ly)/M_LN10;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(l10_val > INT_MAX-1) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l10_val < INT_MIN+1) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double sy = GSL_SIGN(y);
|
|
Packit |
67cb25 |
const int N = (int) floor(l10_val);
|
|
Packit |
67cb25 |
const double arg_val = (l10_val - N) * M_LN10;
|
|
Packit |
67cb25 |
const double arg_err = dy/fabs(y) + dx + 2.0*GSL_DBL_EPSILON*fabs(arg_val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = sy * exp(arg_val);
|
|
Packit |
67cb25 |
result->err = arg_err * fabs(result->val);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
result->e10 = N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_expm1_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double cut = 0.002;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < GSL_LOG_DBL_MIN) {
|
|
Packit |
67cb25 |
result->val = -1.0;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -cut) {
|
|
Packit |
67cb25 |
result->val = exp(x) - 1.0;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < cut) {
|
|
Packit |
67cb25 |
result->val = x * (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
result->val = exp(x) - 1.0;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exprel_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double cut = 0.002;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < GSL_LOG_DBL_MIN) {
|
|
Packit |
67cb25 |
result->val = -1.0/x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -cut) {
|
|
Packit |
67cb25 |
result->val = (exp(x) - 1.0)/x;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < cut) {
|
|
Packit |
67cb25 |
result->val = (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
result->val = (exp(x) - 1.0)/x;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_exprel_2_e(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double cut = 0.002;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < GSL_LOG_DBL_MIN) {
|
|
Packit |
67cb25 |
result->val = -2.0/x*(1.0 + 1.0/x);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -cut) {
|
|
Packit |
67cb25 |
result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < cut) {
|
|
Packit |
67cb25 |
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))));
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_exprel_n_CF_e(const double N, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return exprel_n_CF(N, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_exprel_n_e(const int N, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(N < 0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == 0.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(x) < GSL_ROOT3_DBL_EPSILON * N) {
|
|
Packit |
67cb25 |
result->val = 1.0 + x/(N+1) * (1.0 + x/(N+2));
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(N == 0) {
|
|
Packit |
67cb25 |
return gsl_sf_exp_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(N == 1) {
|
|
Packit |
67cb25 |
return gsl_sf_exprel_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(N == 2) {
|
|
Packit |
67cb25 |
return gsl_sf_exprel_2_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
if(x > N && (-x + N*(1.0 + log(x/N)) < GSL_LOG_DBL_EPSILON)) {
|
|
Packit |
67cb25 |
/* x is much larger than n.
|
|
Packit |
67cb25 |
* Ignore polynomial part, so
|
|
Packit |
67cb25 |
* exprel_N(x) ~= e^x N!/x^N
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result lnf_N;
|
|
Packit |
67cb25 |
double lnr_val;
|
|
Packit |
67cb25 |
double lnr_err;
|
|
Packit |
67cb25 |
double lnterm;
|
|
Packit |
67cb25 |
gsl_sf_lnfact_e(N, &lnf_N);
|
|
Packit |
67cb25 |
lnterm = N*log(x);
|
|
Packit |
67cb25 |
lnr_val = x + lnf_N.val - lnterm;
|
|
Packit |
67cb25 |
lnr_err = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(lnterm));
|
|
Packit |
67cb25 |
lnr_err += lnf_N.err;
|
|
Packit |
67cb25 |
return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > N) {
|
|
Packit |
67cb25 |
/* Write the identity
|
|
Packit |
67cb25 |
* exprel_n(x) = e^x n! / x^n (1 - Gamma[n,x]/Gamma[n])
|
|
Packit |
67cb25 |
* then use the asymptotic expansion
|
|
Packit |
67cb25 |
* Gamma[n,x] ~ x^(n-1) e^(-x) (1 + (n-1)/x + (n-1)(n-2)/x^2 + ...)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double ln_x = log(x);
|
|
Packit |
67cb25 |
gsl_sf_result lnf_N;
|
|
Packit |
67cb25 |
double lg_N;
|
|
Packit |
67cb25 |
double lnpre_val;
|
|
Packit |
67cb25 |
double lnpre_err;
|
|
Packit |
67cb25 |
gsl_sf_lnfact_e(N, &lnf_N); /* log(N!) */
|
|
Packit |
67cb25 |
lg_N = lnf_N.val - log(N); /* log(Gamma(N)) */
|
|
Packit |
67cb25 |
lnpre_val = x + lnf_N.val - N*ln_x;
|
|
Packit |
67cb25 |
lnpre_err = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(N*ln_x));
|
|
Packit |
67cb25 |
lnpre_err += lnf_N.err;
|
|
Packit |
67cb25 |
if(lnpre_val < GSL_LOG_DBL_MAX - 5.0) {
|
|
Packit |
67cb25 |
int stat_eG;
|
|
Packit |
67cb25 |
gsl_sf_result bigG_ratio;
|
|
Packit |
67cb25 |
gsl_sf_result pre;
|
|
Packit |
67cb25 |
int stat_ex = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &pre);
|
|
Packit |
67cb25 |
double ln_bigG_ratio_pre = -x + (N-1)*ln_x - lg_N;
|
|
Packit |
67cb25 |
double bigGsum = 1.0;
|
|
Packit |
67cb25 |
double term = 1.0;
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=1; k
|
|
Packit |
67cb25 |
term *= (N-k)/x;
|
|
Packit |
67cb25 |
bigGsum += term;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
stat_eG = gsl_sf_exp_mult_e(ln_bigG_ratio_pre, bigGsum, &bigG_ratio);
|
|
Packit |
67cb25 |
if(stat_eG == GSL_SUCCESS) {
|
|
Packit |
67cb25 |
result->val = pre.val * (1.0 - bigG_ratio.val);
|
|
Packit |
67cb25 |
result->err = pre.val * (2.0*GSL_DBL_EPSILON + bigG_ratio.err);
|
|
Packit |
67cb25 |
result->err += pre.err * fabs(1.0 - bigG_ratio.val);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_ex;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return stat_eG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > -10.0*N) {
|
|
Packit |
67cb25 |
return exprel_n_CF(N, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* x -> -Inf asymptotic:
|
|
Packit |
67cb25 |
* exprel_n(x) ~ e^x n!/x^n - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
|
|
Packit |
67cb25 |
* ~ - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double sum = 1.0;
|
|
Packit |
67cb25 |
double term = 1.0;
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=1; k
|
|
Packit |
67cb25 |
term *= (N-k)/x;
|
|
Packit |
67cb25 |
sum += term;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = -N/x * sum;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double adx = fabs(dx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x + adx > GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x - adx < GSL_LOG_DBL_MIN) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double ex = exp(x);
|
|
Packit |
67cb25 |
const double edx = exp(adx);
|
|
Packit |
67cb25 |
result->val = ex;
|
|
Packit |
67cb25 |
result->err = ex * GSL_MAX_DBL(GSL_DBL_EPSILON, edx - 1.0/edx);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double adx = fabs(dx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x + adx > INT_MAX - 1) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x - adx < INT_MIN + 1) {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR_E10(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const int N = (int)floor(x/M_LN10);
|
|
Packit |
67cb25 |
const double ex = exp(x-N*M_LN10);
|
|
Packit |
67cb25 |
result->val = ex;
|
|
Packit |
67cb25 |
result->err = ex * (2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) + adx);
|
|
Packit |
67cb25 |
result->e10 = N;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_exp(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_exp_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_exp_mult(const double x, const double y)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_exp_mult_e(x, y, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_expm1(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_expm1_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_exprel(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_exprel_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_exprel_2(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_exprel_2_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_exprel_n(const int n, const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_exprel_n_e(n, x, &result));
|
|
Packit |
67cb25 |
}
|