362 lines
10 KiB
C++
362 lines
10 KiB
C++
/* Translated into C++ by SciPy developers in 2024.
|
|
* Original header with Copyright information appears below.
|
|
*/
|
|
|
|
/* hyperg.c
|
|
*
|
|
* Confluent hypergeometric function
|
|
*
|
|
*
|
|
*
|
|
* SYNOPSIS:
|
|
*
|
|
* double a, b, x, y, hyperg();
|
|
*
|
|
* y = hyperg( a, b, x );
|
|
*
|
|
*
|
|
*
|
|
* DESCRIPTION:
|
|
*
|
|
* Computes the confluent hypergeometric function
|
|
*
|
|
* 1 2
|
|
* a x a(a+1) x
|
|
* F ( a,b;x ) = 1 + ---- + --------- + ...
|
|
* 1 1 b 1! b(b+1) 2!
|
|
*
|
|
* Many higher transcendental functions are special cases of
|
|
* this power series.
|
|
*
|
|
* As is evident from the formula, b must not be a negative
|
|
* integer or zero unless a is an integer with 0 >= a > b.
|
|
*
|
|
* The routine attempts both a direct summation of the series
|
|
* and an asymptotic expansion. In each case error due to
|
|
* roundoff, cancellation, and nonconvergence is estimated.
|
|
* The result with smaller estimated error is returned.
|
|
*
|
|
*
|
|
*
|
|
* ACCURACY:
|
|
*
|
|
* Tested at random points (a, b, x), all three variables
|
|
* ranging from 0 to 30.
|
|
* Relative error:
|
|
* arithmetic domain # trials peak rms
|
|
* IEEE 0,30 30000 1.8e-14 1.1e-15
|
|
*
|
|
* Larger errors can be observed when b is near a negative
|
|
* integer or zero. Certain combinations of arguments yield
|
|
* serious cancellation error in the power series summation
|
|
* and also are not in the region of near convergence of the
|
|
* asymptotic series. An error message is printed if the
|
|
* self-estimated relative error is greater than 1.0e-12.
|
|
*
|
|
*/
|
|
|
|
/*
|
|
* Cephes Math Library Release 2.8: June, 2000
|
|
* Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
|
|
*/
|
|
|
|
#pragma once
|
|
|
|
#include "../config.h"
|
|
#include "../error.h"
|
|
|
|
#include "const.h"
|
|
#include "gamma.h"
|
|
#include "rgamma.h"
|
|
|
|
namespace xsf {
|
|
namespace cephes {
|
|
|
|
namespace detail {
|
|
|
|
/* the `type` parameter determines what converging factor to use */
|
|
XSF_HOST_DEVICE inline double hyp2f0(double a, double b, double x, int type, double *err) {
|
|
double a0, alast, t, tlast, maxt;
|
|
double n, an, bn, u, sum, temp;
|
|
|
|
an = a;
|
|
bn = b;
|
|
a0 = 1.0e0;
|
|
alast = 1.0e0;
|
|
sum = 0.0;
|
|
n = 1.0e0;
|
|
t = 1.0e0;
|
|
tlast = 1.0e9;
|
|
maxt = 0.0;
|
|
|
|
do {
|
|
if (an == 0)
|
|
goto pdone;
|
|
if (bn == 0)
|
|
goto pdone;
|
|
|
|
u = an * (bn * x / n);
|
|
|
|
/* check for blowup */
|
|
temp = std::abs(u);
|
|
if ((temp > 1.0) && (maxt > (std::numeric_limits<double>::max() / temp)))
|
|
goto error;
|
|
|
|
a0 *= u;
|
|
t = std::abs(a0);
|
|
|
|
/* terminating condition for asymptotic series:
|
|
* the series is divergent (if a or b is not a negative integer),
|
|
* but its leading part can be used as an asymptotic expansion
|
|
*/
|
|
if (t > tlast)
|
|
goto ndone;
|
|
|
|
tlast = t;
|
|
sum += alast; /* the sum is one term behind */
|
|
alast = a0;
|
|
|
|
if (n > 200)
|
|
goto ndone;
|
|
|
|
an += 1.0e0;
|
|
bn += 1.0e0;
|
|
n += 1.0e0;
|
|
if (t > maxt)
|
|
maxt = t;
|
|
} while (t > MACHEP);
|
|
|
|
pdone: /* series converged! */
|
|
|
|
/* estimate error due to roundoff and cancellation */
|
|
*err = std::abs(MACHEP * (n + maxt));
|
|
|
|
alast = a0;
|
|
goto done;
|
|
|
|
ndone: /* series did not converge */
|
|
|
|
/* The following "Converging factors" are supposed to improve accuracy,
|
|
* but do not actually seem to accomplish very much. */
|
|
|
|
n -= 1.0;
|
|
x = 1.0 / x;
|
|
|
|
switch (type) { /* "type" given as subroutine argument */
|
|
case 1:
|
|
alast *= (0.5 + (0.125 + 0.25 * b - 0.5 * a + 0.25 * x - 0.25 * n) / x);
|
|
break;
|
|
|
|
case 2:
|
|
alast *= 2.0 / 3.0 - b + 2.0 * a + x - n;
|
|
break;
|
|
|
|
default:;
|
|
}
|
|
|
|
/* estimate error due to roundoff, cancellation, and nonconvergence */
|
|
*err = MACHEP * (n + maxt) + std::abs(a0);
|
|
|
|
done:
|
|
sum += alast;
|
|
return (sum);
|
|
|
|
/* series blew up: */
|
|
error:
|
|
*err = std::numeric_limits<double>::infinity();
|
|
set_error("hyperg", SF_ERROR_NO_RESULT, NULL);
|
|
return (sum);
|
|
}
|
|
|
|
/* asymptotic formula for hypergeometric function:
|
|
*
|
|
* ( -a
|
|
* -- ( |z|
|
|
* | (b) ( -------- 2f0( a, 1+a-b, -1/x )
|
|
* ( --
|
|
* ( | (b-a)
|
|
*
|
|
*
|
|
* x a-b )
|
|
* e |x| )
|
|
* + -------- 2f0( b-a, 1-a, 1/x ) )
|
|
* -- )
|
|
* | (a) )
|
|
*/
|
|
|
|
XSF_HOST_DEVICE inline double hy1f1a(double a, double b, double x, double *err) {
|
|
double h1, h2, t, u, temp, acanc, asum, err1, err2;
|
|
|
|
if (x == 0) {
|
|
acanc = 1.0;
|
|
asum = std::numeric_limits<double>::infinity();
|
|
goto adone;
|
|
}
|
|
temp = std::log(std::abs(x));
|
|
t = x + temp * (a - b);
|
|
u = -temp * a;
|
|
|
|
if (b > 0) {
|
|
temp = xsf::cephes::lgam(b);
|
|
t += temp;
|
|
u += temp;
|
|
}
|
|
|
|
h1 = hyp2f0(a, a - b + 1, -1.0 / x, 1, &err1);
|
|
|
|
temp = std::exp(u) * xsf::cephes::rgamma(b - a);
|
|
h1 *= temp;
|
|
err1 *= temp;
|
|
|
|
h2 = hyp2f0(b - a, 1.0 - a, 1.0 / x, 2, &err2);
|
|
|
|
if (a < 0)
|
|
temp = std::exp(t) * xsf::cephes::rgamma(a);
|
|
else
|
|
temp = std::exp(t - xsf::cephes::lgam(a));
|
|
|
|
h2 *= temp;
|
|
err2 *= temp;
|
|
|
|
if (x < 0.0)
|
|
asum = h1;
|
|
else
|
|
asum = h2;
|
|
|
|
acanc = std::abs(err1) + std::abs(err2);
|
|
|
|
if (b < 0) {
|
|
temp = xsf::cephes::Gamma(b);
|
|
asum *= temp;
|
|
acanc *= std::abs(temp);
|
|
}
|
|
|
|
if (asum != 0.0)
|
|
acanc /= std::abs(asum);
|
|
|
|
if (acanc != acanc)
|
|
/* nan */
|
|
acanc = 1.0;
|
|
|
|
if (std::isinf(asum))
|
|
/* infinity */
|
|
acanc = 0;
|
|
|
|
acanc *= 30.0; /* fudge factor, since error of asymptotic formula
|
|
* often seems this much larger than advertised */
|
|
adone:
|
|
*err = acanc;
|
|
return (asum);
|
|
}
|
|
|
|
/* Power series summation for confluent hypergeometric function */
|
|
XSF_HOST_DEVICE inline double hy1f1p(double a, double b, double x, double *err) {
|
|
double n, a0, sum, t, u, temp, maxn;
|
|
double an, bn, maxt;
|
|
double y, c, sumc;
|
|
|
|
/* set up for power series summation */
|
|
an = a;
|
|
bn = b;
|
|
a0 = 1.0;
|
|
sum = 1.0;
|
|
c = 0.0;
|
|
n = 1.0;
|
|
t = 1.0;
|
|
maxt = 0.0;
|
|
*err = 1.0;
|
|
|
|
maxn = 200.0 + 2 * fabs(a) + 2 * fabs(b);
|
|
|
|
while (t > MACHEP) {
|
|
if (bn == 0) { /* check bn first since if both */
|
|
sf_error("hyperg", SF_ERROR_SINGULAR, NULL);
|
|
return (std::numeric_limits<double>::infinity()); /* an and bn are zero it is */
|
|
}
|
|
if (an == 0) /* a singularity */
|
|
return (sum);
|
|
if (n > maxn) {
|
|
/* too many terms; take the last one as error estimate */
|
|
c = std::abs(c) + std::abs(t) * 50.0;
|
|
goto pdone;
|
|
}
|
|
u = x * (an / (bn * n));
|
|
|
|
/* check for blowup */
|
|
temp = std::abs(u);
|
|
if ((temp > 1.0) && (maxt > (std::numeric_limits<double>::max() / temp))) {
|
|
*err = 1.0; /* blowup: estimate 100% error */
|
|
return sum;
|
|
}
|
|
|
|
a0 *= u;
|
|
|
|
y = a0 - c;
|
|
sumc = sum + y;
|
|
c = (sumc - sum) - y;
|
|
sum = sumc;
|
|
|
|
t = std::abs(a0);
|
|
|
|
an += 1.0;
|
|
bn += 1.0;
|
|
n += 1.0;
|
|
}
|
|
|
|
pdone:
|
|
|
|
/* estimate error due to roundoff and cancellation */
|
|
if (sum != 0.0) {
|
|
*err = std::abs(c / sum);
|
|
} else {
|
|
*err = std::abs(c);
|
|
}
|
|
|
|
if (*err != *err) {
|
|
/* nan */
|
|
*err = 1.0;
|
|
}
|
|
|
|
return (sum);
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
XSF_HOST_DEVICE inline double hyperg(double a, double b, double x) {
|
|
double asum, psum, acanc, pcanc, temp;
|
|
|
|
/* See if a Kummer transformation will help */
|
|
temp = b - a;
|
|
if (std::abs(temp) < 0.001 * std::abs(a))
|
|
return (exp(x) * hyperg(temp, b, -x));
|
|
|
|
/* Try power & asymptotic series, starting from the one that is likely OK */
|
|
if (std::abs(x) < 10 + std::abs(a) + std::abs(b)) {
|
|
psum = detail::hy1f1p(a, b, x, &pcanc);
|
|
if (pcanc < 1.0e-15)
|
|
goto done;
|
|
asum = detail::hy1f1a(a, b, x, &acanc);
|
|
} else {
|
|
psum = detail::hy1f1a(a, b, x, &pcanc);
|
|
if (pcanc < 1.0e-15)
|
|
goto done;
|
|
asum = detail::hy1f1p(a, b, x, &acanc);
|
|
}
|
|
|
|
/* Pick the result with less estimated error */
|
|
|
|
if (acanc < pcanc) {
|
|
pcanc = acanc;
|
|
psum = asum;
|
|
}
|
|
|
|
done:
|
|
if (pcanc > 1.0e-12)
|
|
set_error("hyperg", SF_ERROR_LOSS, NULL);
|
|
|
|
return (psum);
|
|
}
|
|
|
|
} // namespace cephes
|
|
} // namespace xsf
|