164 lines
4.5 KiB
C++
164 lines
4.5 KiB
C++
/* Translated into C++ by SciPy developers in 2024.
|
|
* Original header with Copyright information appears below.
|
|
*/
|
|
|
|
/* k1.c
|
|
*
|
|
* Modified Bessel function, third kind, order one
|
|
*
|
|
*
|
|
*
|
|
* SYNOPSIS:
|
|
*
|
|
* double x, y, k1();
|
|
*
|
|
* y = k1( x );
|
|
*
|
|
*
|
|
*
|
|
* DESCRIPTION:
|
|
*
|
|
* Computes the modified Bessel function of the third kind
|
|
* of order one of the argument.
|
|
*
|
|
* The range is partitioned into the two intervals [0,2] and
|
|
* (2, infinity). Chebyshev polynomial expansions are employed
|
|
* in each interval.
|
|
*
|
|
*
|
|
*
|
|
* ACCURACY:
|
|
*
|
|
* Relative error:
|
|
* arithmetic domain # trials peak rms
|
|
* IEEE 0, 30 30000 1.2e-15 1.6e-16
|
|
*
|
|
* ERROR MESSAGES:
|
|
*
|
|
* message condition value returned
|
|
* k1 domain x <= 0 INFINITY
|
|
*
|
|
*/
|
|
/* k1e.c
|
|
*
|
|
* Modified Bessel function, third kind, order one,
|
|
* exponentially scaled
|
|
*
|
|
*
|
|
*
|
|
* SYNOPSIS:
|
|
*
|
|
* double x, y, k1e();
|
|
*
|
|
* y = k1e( x );
|
|
*
|
|
*
|
|
*
|
|
* DESCRIPTION:
|
|
*
|
|
* Returns exponentially scaled modified Bessel function
|
|
* of the third kind of order one of the argument:
|
|
*
|
|
* k1e(x) = exp(x) * k1(x).
|
|
*
|
|
*
|
|
*
|
|
* ACCURACY:
|
|
*
|
|
* Relative error:
|
|
* arithmetic domain # trials peak rms
|
|
* IEEE 0, 30 30000 7.8e-16 1.2e-16
|
|
* See k1().
|
|
*
|
|
*/
|
|
|
|
/*
|
|
* Cephes Math Library Release 2.8: June, 2000
|
|
* Copyright 1984, 1987, 2000 by Stephen L. Moshier
|
|
*/
|
|
#pragma once
|
|
|
|
#include "../config.h"
|
|
#include "../error.h"
|
|
|
|
#include "chbevl.h"
|
|
#include "const.h"
|
|
|
|
namespace xsf {
|
|
namespace cephes {
|
|
|
|
namespace detail {
|
|
/* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
|
|
* in the interval [0,2].
|
|
*
|
|
* lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
|
|
*/
|
|
|
|
constexpr double k1_A[] = {
|
|
-7.02386347938628759343E-18, -2.42744985051936593393E-15, -6.66690169419932900609E-13,
|
|
-1.41148839263352776110E-10, -2.21338763073472585583E-8, -2.43340614156596823496E-6,
|
|
-1.73028895751305206302E-4, -6.97572385963986435018E-3, -1.22611180822657148235E-1,
|
|
-3.53155960776544875667E-1, 1.52530022733894777053E0};
|
|
|
|
/* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
|
|
* in the interval [2,infinity].
|
|
*
|
|
* lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
|
|
*/
|
|
constexpr double k1_B[] = {
|
|
-5.75674448366501715755E-18, 1.79405087314755922667E-17, -5.68946255844285935196E-17,
|
|
1.83809354436663880070E-16, -6.05704724837331885336E-16, 2.03870316562433424052E-15,
|
|
-7.01983709041831346144E-15, 2.47715442448130437068E-14, -8.97670518232499435011E-14,
|
|
3.34841966607842919884E-13, -1.28917396095102890680E-12, 5.13963967348173025100E-12,
|
|
-2.12996783842756842877E-11, 9.21831518760500529508E-11, -4.19035475934189648750E-10,
|
|
2.01504975519703286596E-9, -1.03457624656780970260E-8, 5.74108412545004946722E-8,
|
|
-3.50196060308781257119E-7, 2.40648494783721712015E-6, -1.93619797416608296024E-5,
|
|
1.95215518471351631108E-4, -2.85781685962277938680E-3, 1.03923736576817238437E-1,
|
|
2.72062619048444266945E0};
|
|
|
|
} // namespace detail
|
|
|
|
XSF_HOST_DEVICE inline double k1(double x) {
|
|
double y, z;
|
|
|
|
if (x == 0.0) {
|
|
set_error("k1", SF_ERROR_SINGULAR, NULL);
|
|
return std::numeric_limits<double>::infinity();
|
|
} else if (x < 0.0) {
|
|
set_error("k1", SF_ERROR_DOMAIN, NULL);
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
}
|
|
z = 0.5 * x;
|
|
|
|
if (x <= 2.0) {
|
|
y = x * x - 2.0;
|
|
y = std::log(z) * i1(x) + chbevl(y, detail::k1_A, 11) / x;
|
|
return (y);
|
|
}
|
|
|
|
return (std::exp(-x) * chbevl(8.0 / x - 2.0, detail::k1_B, 25) / std::sqrt(x));
|
|
}
|
|
|
|
XSF_HOST_DEVICE double k1e(double x) {
|
|
double y;
|
|
|
|
if (x == 0.0) {
|
|
set_error("k1e", SF_ERROR_SINGULAR, NULL);
|
|
return std::numeric_limits<double>::infinity();
|
|
} else if (x < 0.0) {
|
|
set_error("k1e", SF_ERROR_DOMAIN, NULL);
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
}
|
|
|
|
if (x <= 2.0) {
|
|
y = x * x - 2.0;
|
|
y = std::log(0.5 * x) * i1(x) + chbevl(y, detail::k1_A, 11) / x;
|
|
return (y * exp(x));
|
|
}
|
|
|
|
return (chbevl(8.0 / x - 2.0, detail::k1_B, 25) / std::sqrt(x));
|
|
}
|
|
|
|
} // namespace cephes
|
|
} // namespace xsf
|