The RcppBessel package exports an Rcpp interface for the
Bessel functions found in the Bessel package of Maechler by
wrapping the functions in the C code (translated and cleaned up from
Amos’ original FORTRAN code by Maechler) into C++. Not all functionality
is exposed, in order to keep things simple, but may be expanded in the
future to export more of the functions such as those for the asymptotic
expansion of BesselK and BesselI for large nu and x. The package also
exports the functions for use in R, as it was found to offer a speed up
of up to 10x.
The initial motivation for writing this package was so as to allow the modified Bessel function of the second kind to be used in C++ code which the author needed in another package. The Boost implementation does not allow for complex arguments whilst the special math functions in std library is not available for Mac OS (via the shipped libc++). Additionally, neither of those implementation has the option of exponential scaling.
The table below shows the functions currently available.
| RcppBessel Function | Description |
|---|---|
bessel_k |
Computes the modified Bessel function of the second kind. |
bessel_y |
Computes the Bessel function of the second kind (Neumann function). |
bessel_h |
Computes the Hankel function (Bessel function of the third kind). |
bessel_i |
Computes the modified Bessel function of the first kind. |
bessel_j |
Computes the Bessel function of the first kind. |
airy_a |
Computes the Airy function of the first kind. |
airy_b |
Computes the Airy function of the second kind. |
To distinguish between the CamelCase naming in the
Bessel package and the lowercamelCase naming in base R, the
package has adopted the snake_case syntax similar to the
std lib. Other than that, the arguments are the same as in
the Bessel package with the exception of nSeq, which
computes the result for a whole sequence of nu values, which we have
opted to skip in the implementation.
For real values less than zero, bessel_k and
bessel_y return complex values. In this case, if the vector
z has negative values, the whole output returned is complex
valued. The table below shows the general input/output types
expected.
| Function Name | Input for z |
Output Type |
|---|---|---|
bessel_k |
Real (all >= 0) | Real |
| Real (any < 0) | Complex | |
| Complex | Complex | |
bessel_y |
Real (all >= 0) | Real |
| Real (some < 0) | Complex | |
| Complex | Complex | |
bessel_h |
Real (any value) | Complex |
| Complex | Complex | |
bessel_i |
Real (any value) | Real |
| Complex | Complex | |
bessel_j |
Real (any value) | Real |
| Complex | Complex | |
airy_a |
Real (any value) | Real |
| Complex | Complex | |
airy_b |
Real (any value) | Real |
| Complex | Complex |
For values of zero, the output is similar to what the
Bessel package returns with the exception of
bessel_k which returns Inf similar to base R,
whereas Bessel::BesselK returns 0.
The next sections provide a individual summary of each of the
functions and a sample C++ code which can be pasted and
sourced to test the results against the Bessel function
implementation and also show how the package can be used inside other
Rcpp functions.
bessel_k)Computes the modified Bessel function of the second kind.
z: A numeric or complex vector of points at which to
evaluate the Bessel function.nu: A numeric value representing the order of the
Bessel function.expon_scaled: A logical value indicating whether to use
the exponentially scaled form.verbose: An integer specifying the verbosity level for
error messages.bessel_k function.// test_bessel_k.cpp
#include <Rcpp.h>
#include <RcppBessel.h>
// [[Rcpp::depends(Rcpp, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_bessel_k(SEXP x, double nu, bool expon_scaled, int verbose) {
SEXP z = RcppBessel::bessel_k(x, nu, expon_scaled, verbose);
return z;
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
library(Bessel)
# Source the C++ code
Rcpp::sourceCpp("test_bessel_k.cpp")
# Generate a sequence
x <- seq(1, 10, by = 1)
# Compare results with Bessel package
tmp <- Bessel::BesselK(x, 2, FALSE, 1, 0)
all.equal(tmp, test_bessel_k(x, 2, FALSE, 0))
*/bessel_i)Computes the modified Bessel function of the first kind.
z: A numeric or complex vector of points at which to
evaluate the Bessel function.nu: A numeric value representing the order of the
Bessel function.expon_scaled: A logical value indicating whether to use
the exponentially scaled form.verbose: An integer specifying the verbosity level for
error messages.bessel_i function.// test_bessel_i.cpp
#include <Rcpp.h>
#include <RcppBessel.h>
// [[Rcpp::depends(Rcpp, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_bessel_i(SEXP x, double nu, bool expon_scaled, int verbose) {
SEXP z = RcppBessel::bessel_i(x, nu, expon_scaled, verbose);
return z;
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
library(Bessel)
# Source the C++ code
Rcpp::sourceCpp("test_bessel_i.cpp")
# Generate a sequence
x <- seq(-10, 10, by = 1)
# Compare results with Bessel package
tmp <- Bessel::BesselI(x, 2, FALSE, 1, 0)
all.equal(tmp, test_bessel_i(x, 2, FALSE, 0))
*/bessel_j)Computes the Bessel function of the first kind.
z: A numeric or complex vector of points at which to
evaluate the Bessel function.nu: A numeric value representing the order of the
Bessel function.expon_scaled: A logical value indicating whether to use
the exponentially scaled form.verbose: An integer specifying the verbosity level for
error messages.bessel_j function.// test_bessel_j.cpp
#include <Rcpp.h>
#include <RcppBessel.h>
// [[Rcpp::depends(Rcpp, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_bessel_j(SEXP x, double nu, bool expon_scaled, int verbose) {
SEXP z = RcppBessel::bessel_j(x, nu, expon_scaled, verbose);
return z;
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
library(Bessel)
# Source the C++ code
Rcpp::sourceCpp("test_bessel_j.cpp")
# Generate a sequence
x <- seq(-10, 10, by = 1)
# Compare results with Bessel package
tmp <- Bessel::BesselJ(x, 2, FALSE, 1, 0)
all.equal(tmp, test_bessel_j(x, 2, FALSE, 0))
*/bessel_y)Computes the Bessel function of the second kind (Neumann function).
z: A (strictly positive) numeric or complex vector of
points at which to evaluate the Bessel function.nu: A numeric value representing the order of the
Bessel function.expon_scaled: A logical value indicating whether to use
the exponentially scaled form.verbose: An integer specifying the verbosity level for
error messages.bessel_y function.// test_bessel_y.cpp
#include <Rcpp.h>
#include <RcppBessel.h>
// [[Rcpp::depends(Rcpp, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_bessel_y(SEXP x, double nu, bool expon_scaled, int verbose) {
SEXP z = RcppBessel::bessel_y(x, nu, expon_scaled, verbose);
return z;
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
library(Bessel)
# Source the C++ code
Rcpp::sourceCpp("test_bessel_y.cpp")
# Generate a sequence
x <- seq(-10, 10, by = 0.35)
# skipping zero since bessel_y returns -Inf+0i whereas BesselY return -Inf+NaNi
# Compare results with Bessel package
tmp <- Bessel::BesselY(x, 2, FALSE, 1, 0)
all.equal(tmp, test_bessel_y(x, 2, FALSE, 0))
*/bessel_h)Computes the Hankel function (Bessel function of the third kind) of the first or second kind.
m: An integer specifying the type of Hankel function.
Must be either 1 (for the first kind) or 2
(for the second kind).z: A numeric or complex vector of points at which to
evaluate the Hankel function.nu: A numeric value representing the order of the
Hankel function.expon_scaled: A logical value indicating whether to use
the exponentially scaled form.verbose: An integer specifying the verbosity level for
error messages.bessel_h
function evaluated at the points in z.// test_bessel_h.cpp
#include <Rcpp.h>
#include <RcppBessel.h>
// [[Rcpp::depends(Rcpp, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_bessel_h(int m, SEXP x, double nu, bool expon_scaled, int verbose) {
SEXP z = RcppBessel::bessel_h(m, x, nu, expon_scaled, verbose);
return z;
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
library(Bessel)
# Source the C++ code
Rcpp::sourceCpp("test_bessel_h.cpp")
# Generate a sequence
x <- seq(-10, 10, by = 0.35)
# skipping again zero
# Compare results with Bessel package
tmp <- Bessel::BesselH(1, x, 2, FALSE, 1, 0)
all.equal(tmp, test_bessel_h(1, x, 2, FALSE, 0))
*/airy_a)Computes the Airy function Ai for real or complex inputs.
z: A numeric or complex vector of points at which to
evaluate the Airy function.deriv: An integer indicating whether to compute the
function itself (0) or its first derivative
(1). Defaults to 0.expon_scaled: A logical value indicating whether to use
the exponentially scaled form of the Airy function. Defaults to
FALSE.verbose: An integer specifying the verbosity level for
error messages.airy_a function evaluated at the points in
z.// test_airy_a.cpp
#include <Rcpp.h>
#include <RcppBessel.h>
// [[Rcpp::depends(Rcpp, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_airy_a(SEXP x, int deriv = 0, bool expon_scaled = false, int verbose = 0) {
SEXP z = RcppBessel::airy_a(x, deriv, expon_scaled, verbose);
return z;
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
# Source the C++ code
Rcpp::sourceCpp("test_airy_a.cpp")
# Generate a sequence
x <- seq(-10, 10, by = 1)
tmp <- Bessel::AiryA(x, 1)
all.equal(tmp, test_airy_a(x, 1, FALSE, 0))
*/airy_b)Computes the Airy function Bi for real or complex inputs.
z: A numeric or complex vector of points at which to
evaluate the Airy function.deriv: An integer indicating whether to compute the
function itself (0) or its first derivative
(1). Defaults to 0.expon_scaled: A logical value indicating whether to use
the exponentially scaled form of the Airy function. Defaults to
FALSE.verbose: An integer specifying the verbosity level for
error messages.airy_b function evaluated at the points in
z.// test_airy_b.cpp
#include <Rcpp.h>
#include <RcppBessel.h>
// [[Rcpp::depends(Rcpp, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_airy_b(SEXP x, int deriv = 0, bool expon_scaled = false, int verbose = 0) {
SEXP z = RcppBessel::airy_b(x, deriv, expon_scaled, verbose);
return z;
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
# Source the C++ code
Rcpp::sourceCpp("test_airy_b.cpp")
# Generate a sequence
x <- seq(-10, 10, by = 1)
tmp <- Bessel::AiryB(x, 1)
all.equal(tmp, test_airy_b(x, 1, FALSE, 0))
*/#include <RcppArmadillo.h>
#include <RcppBessel.h>
// [[Rcpp::depends(RcppArmadillo, RcppBessel)]]
// [[Rcpp::export]]
SEXP test_bessel_k(const arma::vec& x, double nu, bool expon_scaled = false, int verbose = 0) {
// Convert arma::vec to Rcpp::NumericVector
Rcpp::NumericVector rx = Rcpp::wrap(x);
// Call the BesselK function
SEXP z = RcppBessel::bessel_k(rx, nu, expon_scaled, verbose);
// Check if the result is complex
if (Rcpp::is<Rcpp::ComplexVector>(z)) {
Rcpp::ComplexVector cz(z);
arma::cx_vec result = Rcpp::as<arma::cx_vec>(cz);
return Rcpp::wrap(result);
} else {
Rcpp::NumericVector rz(z);
arma::vec result = Rcpp::as<arma::vec>(rz);
return Rcpp::wrap(result);
}
}
/*** R
# Example usage in R
library(Rcpp)
library(RcppBessel)
library(Bessel)
library(RcppArmadillo)
# Source the C++ code
#Rcpp::sourceCpp("test_bessel_k.cpp")
# Generate a sequence
x <- seq(1, 10, by = 1)
# Compare results with Bessel package
tmp <- Bessel::BesselK(x, 2, FALSE, 1, 0)
all.equal(tmp, as.vector(test_bessel_k(x, 2, FALSE, 0)))
*/The following code and output provides a quick timing test between
the Bessel and RcppBessel packages for the
bessel_k function.
library(microbenchmark)
library(Bessel)
library(RcppBessel)
z <- seq(-20, 20, by = 0.001)
b <- microbenchmark(BesselK(z, 2, TRUE), bessel_k(z, 2, TRUE), times = 5L)
#> Warning in microbenchmark(BesselK(z, 2, TRUE), bessel_k(z, 2, TRUE), times =
#> 5L): less accurate nanosecond times to avoid potential integer overflows
print(b)
#> Unit: milliseconds
#> expr min lq mean median uq max
#> BesselK(z, 2, TRUE) 94.64010 95.94836 96.72466 96.09502 98.44621 98.49360
#> bessel_k(z, 2, TRUE) 15.50489 15.53982 15.55122 15.54503 15.55532 15.61104
#> neval
#> 5
#> 5