Skip to content

Commit

Permalink
Revisit bessel functions
Browse files Browse the repository at this point in the history
  • Loading branch information
jtlap authored Dec 16, 2024
1 parent ac528c9 commit 0189a47
Show file tree
Hide file tree
Showing 183 changed files with 2,527 additions and 9,091 deletions.
60 changes: 6 additions & 54 deletions include/kyosu/bessel.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,60 +13,12 @@
#include <kyosu/math.hpp>
#include <kyosu/details/bessel.hpp>

#include <kyosu/functions/bessel_j.hpp>
#include <kyosu/functions/bessel_y.hpp>
#include <kyosu/functions/bessel_h.hpp>
#include <kyosu/functions/bessel_i.hpp>
#include <kyosu/functions/bessel_k.hpp>

#include <kyosu/functions/airy.hpp>
#include <kyosu/functions/airy_ai.hpp>
#include <kyosu/functions/airy_bi.hpp>

#include <kyosu/functions/cyl_bessel_j0.hpp>
#include <kyosu/functions/cyl_bessel_j1.hpp>
#include <kyosu/functions/cyl_bessel_jn.hpp>
#include <kyosu/functions/cyl_bessel_y0.hpp>
#include <kyosu/functions/cyl_bessel_y1.hpp>
#include <kyosu/functions/cyl_bessel_yn.hpp>
#include <kyosu/functions/cyl_bessel_h1_0.hpp>
#include <kyosu/functions/cyl_bessel_h1_1.hpp>
#include <kyosu/functions/cyl_bessel_h1n.hpp>
#include <kyosu/functions/cyl_bessel_h2_0.hpp>
#include <kyosu/functions/cyl_bessel_h2_1.hpp>
#include <kyosu/functions/cyl_bessel_h2n.hpp>
#include <kyosu/functions/cyl_bessel_h12n.hpp>
#include <kyosu/functions/cyl_bessel_i0.hpp>
#include <kyosu/functions/cyl_bessel_i1.hpp>
#include <kyosu/functions/cyl_bessel_in.hpp>
#include <kyosu/functions/cyl_bessel_k0.hpp>
#include <kyosu/functions/cyl_bessel_k1.hpp>
#include <kyosu/functions/cyl_bessel_kn.hpp>
#include <kyosu/functions/cyl_bessel_jyn.hpp>
#include <kyosu/functions/cyl_bessel_ikn.hpp>

#include <kyosu/functions/cyl_bessel_j.hpp>
#include <kyosu/functions/cyl_bessel_y.hpp>
#include <kyosu/functions/cyl_bessel_jy.hpp>
#include <kyosu/functions/cyl_bessel_h1.hpp>
#include <kyosu/functions/cyl_bessel_h2.hpp>
#include <kyosu/functions/cyl_bessel_h12.hpp>
#include <kyosu/functions/cyl_bessel_i.hpp>
#include <kyosu/functions/cyl_bessel_k.hpp>
#include <kyosu/functions/cyl_bessel_ik.hpp>

#include <kyosu/functions/sph_bessel_j0.hpp>
#include <kyosu/functions/sph_bessel_j1.hpp>
#include <kyosu/functions/sph_bessel_jn.hpp>
#include <kyosu/functions/sph_bessel_y0.hpp>
#include <kyosu/functions/sph_bessel_y1.hpp>
#include <kyosu/functions/sph_bessel_yn.hpp>
#include <kyosu/functions/sph_bessel_h1_0.hpp>
#include <kyosu/functions/sph_bessel_h1_1.hpp>
#include <kyosu/functions/sph_bessel_h1n.hpp>
#include <kyosu/functions/sph_bessel_h2_0.hpp>
#include <kyosu/functions/sph_bessel_h2_1.hpp>
#include <kyosu/functions/sph_bessel_h2n.hpp>
#include <kyosu/functions/sph_bessel_i1_0.hpp>
#include <kyosu/functions/sph_bessel_i1_1.hpp>
#include <kyosu/functions/sph_bessel_i1n.hpp>
#include <kyosu/functions/sph_bessel_i2_0.hpp>
#include <kyosu/functions/sph_bessel_i2_1.hpp>
#include <kyosu/functions/sph_bessel_i2n.hpp>
#include <kyosu/functions/sph_bessel_k0.hpp>
#include <kyosu/functions/sph_bessel_k1.hpp>
#include <kyosu/functions/sph_bessel_kn.hpp>
41 changes: 20 additions & 21 deletions include/kyosu/details/bessel.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@
//!
//! spherical functions
//!
//! As in DMLF we restrict the implementation to non-negative integral orders
//! Contrarily to DMLF we do not restrict the implementation to non-negative integral orders
//!
//! The computation scheme of \f$j_n\f$, \f$y_n\f$, \f$h_n^{(1)}\f$ and \f$h_n^{(2)}\f$
//! follows the article of Liang-Wu Kai (On the computation of spherical Bessel functions of complex arguments,
Expand Down Expand Up @@ -62,45 +62,41 @@
//!
//! spherical functions
//!
//! Following DLMF 10.47 spherical functions are only implemented for integral non negative orders.
//! Contrarily to DMLF we do not restrict the implementation to non-negative integral orders
//! Non-integral and negative orders are computed via their relation to cylindrical
//!
//!======================================================================================================================
//!
//! Rationale of the interface
//!
//! the functions are all named [cyl, sph]_bessel_xxx xxx designing the chosen bessel function(s).
//! the functions are all named bessel_x x designing the chosen bessel function(s) i.e. j, i, h, k.
//! * the options `cylindrical` and `spherical` are always available.
//! * the options `kind_1` and `kind_2` are also available for some functions : i.e.
//! no option means `cylindrical` and `kind_1`
//!
//! The choice not to make order an simd parameter is driven by two causes
//!
//! * The complexity involved in the merging of test on order and on evaluation arguments preclude the
//! SIMD expected performances
//! * many bessel functions of arithmetic range order (with step one) can be computed with almost no additionnal cost
//! provided that the returned value is complex (scalar or simd). So for each possible call another is provided
//! allowing additionnal output parameters, the rule being that if aaa_bessel_xxx(n, z, output) is called:
//! allowing additionnal output parameters, the rule being that if bessel_x(n, z, output) is called:
//!
//! * the caller must provide output as a random access range of kyosu::complex allocated to receive
//! N = int(trunc(abs(n)+1)) values.
//! * the caller must provide output as a std::span of kyosu::complex allocated to receive N values
//! values.
//! * These values will be the same as the result of calling :
//! `for (int i = 0; i <= N; ++i) output[i] = aaa_bessel_xxx(n0+sgn*i, z)`
//! `for (int i = 0; i < N; ++i) output[i] = bessel_x(n0+sgn*i, z)`
//! where n0 is the fractionnal part of n and sgn is the sign of n.
//!
//! * Moreover jy and ik functions can be computed together and two ouput parameters can be used
//!
//! values of xxx
//!
//! * one parameter
//! j0 j1 i0 i1 y0 y1 k0 k1 h10 h11 h20 h21 only admit the z parameter which can be any cayley-dickson type
//! * two parameters or three parameters
//! * jn in yn kn h1n h2n the first one is the integer order, the second z which can be any cayley-dickson type
//! * j i y k h1 h2 the first one is the floating order, the second z which can be any cayley-dickson type
//! * if present the third parameter will be the output parameter as described above and z must be floating or complex
//! * two parameters or four parameters
//! * jyn ikn h12n the first one is the integer order, the second z which can be any cayley-dickson type
//! * jy ik h12 the first one is the floating order, the second z which can be any cayley-dickson type
//! * if present the third and fourth parameter will be the output parameters as described above
//! and z must be floating or complex
//!======================================================================================================================

#include <kyosu/types/impl/detail/bessel_utils.hpp>

#include <kyosu/types/impl/detail/bessel_utils.hpp>
// These files only contain implementation details and do not contain any function
// belonging to the user interface
//
// bessels of integral order
#include <kyosu/types/impl/bessel/cb_jyn.hpp>
#include <kyosu/types/impl/bessel/sb_jyn.hpp>
Expand All @@ -112,5 +108,8 @@
#include <kyosu/types/impl/besselr/cb_jyr.hpp>
#include <kyosu/types/impl/besselr/cb_hr.hpp>
#include <kyosu/types/impl/besselr/cb_ikr.hpp>
#include <kyosu/types/impl/besselr/sb_jyr.hpp>
#include <kyosu/types/impl/besselr/sb_hr.hpp>
#include <kyosu/types/impl/besselr/sb_ikr.hpp>
// // airy functions
#include <kyosu/types/impl/besselr/airy.hpp>
8 changes: 8 additions & 0 deletions include/kyosu/details/decorators.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,13 +13,21 @@ namespace kyosu
struct assume_unitary_mode {};
struct intrinsic_mode {};
struct extrinsic_mode {};
struct kind_1_mode {};
struct kind_2_mode {};

[[maybe_unused]] inline constexpr auto assume_unitary = ::rbr::flag(assume_unitary_mode{});
[[maybe_unused]] inline constexpr auto intrinsic = ::rbr::flag(intrinsic_mode{});
[[maybe_unused]] inline constexpr auto extrinsic = ::rbr::flag(extrinsic_mode{});
[[maybe_unused]] inline constexpr auto kind_1 = ::rbr::flag(kind_1_mode{});
[[maybe_unused]] inline constexpr auto kind_2 = ::rbr::flag(kind_2_mode{});


struct assume_unitary_option : eve::detail::exact_option<assume_unitary> {};
struct intrinsic_option : eve::detail::exact_option<intrinsic> {};
struct extrinsic_option : eve::detail::exact_option<extrinsic> {};
struct kind_1_option : eve::detail::exact_option<kind_1> {};
struct kind_2_option : eve::detail::exact_option<kind_2> {};


}
1 change: 1 addition & 0 deletions include/kyosu/functions.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@
#include <kyosu/functions/atanh.hpp>
#include <kyosu/functions/atanpi.hpp>
#include <kyosu/functions/average.hpp>

#include <kyosu/functions/beta.hpp>
#include <kyosu/functions/ceil.hpp>
#include <kyosu/functions/chi.hpp>
Expand Down
7 changes: 4 additions & 3 deletions include/kyosu/functions/airy.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,8 @@
#pragma once
#include <kyosu/details/callable.hpp>
#include <kyosu/constants/wrapped.hpp>
#include <kyosu/functions/cyl_bessel_i.hpp>
#include <kyosu/types/impl/besselr/cb_ikr.hpp>


namespace kyosu
{
Expand All @@ -28,8 +29,8 @@ namespace kyosu
};

auto [sqzo3, zeta] = zet(z);
auto ip = cyl_bessel_i(eve::third(as<u_t>()), zeta);
auto im = cyl_bessel_i(-eve::third(as<u_t>()), zeta);
auto ip = _::cb_ir(eve::third(as<u_t>()), zeta);
auto im = _::cb_ir(-eve::third(as<u_t>()), zeta);
return kumi::tuple{inv_pi(as<u_t>())*sqzo3*(im-ip), sqzo3*(ip+im)};
}
else
Expand Down
151 changes: 151 additions & 0 deletions include/kyosu/functions/bessel_h.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
//======================================================================================================================
/*
Kyosu - Complex Without Complexes
Copyright: KYOSU Contributors & Maintainers
SPDX-License-Identifier: BSL-1.0
*/
//======================================================================================================================
#pragma once
#include "eve/traits/as_logical.hpp"
#include <kyosu/details/callable.hpp>
#include <kyosu/bessel.hpp>


namespace kyosu
{

template<typename Options>
struct bessel_h_t : eve::strict_elementwise_callable<bessel_h_t, Options, eve::spherical_option, eve::cylindrical_option, kind_1_option, kind_2_option>
{
template<eve::scalar_value N, typename Z>
requires(concepts::real<Z> || concepts::cayley_dickson<Z>)
KYOSU_FORCEINLINE constexpr auto operator()(N const& n, Z const & z) const noexcept
{
constexpr int Kind = Options::contains(kind_2) ? 2 : 1;
if constexpr(concepts::complex<Z> )
{
if constexpr(eve::integral_scalar_value<N>)
{
if constexpr(Options::contains(eve::spherical))
{
if (eve::is_ltz(n))
{
using u_t = eve::underlying_type_t<Z>;
return _::sb_hr<Kind>(u_t(n), z);
}
else
return _::sb_hn<Kind>(n, z);
}
else
return _::cb_hn<Kind>(n, z);
}
else if constexpr( eve::floating_scalar_value<N>)
{
if constexpr(Options::contains(eve::spherical))
return _::sb_hr<Kind>(n, z);
else
return _::cb_hr<Kind>(n, z);
}
}
else
return _::cayley_extend_rev(*this, n, z);
}

template<eve::scalar_value N, typename Z, std::size_t S>
requires(concepts::real<Z> || concepts::complex<Z>)
KYOSU_FORCEINLINE constexpr auto operator()(N const& n, Z const & z, std::span<Z, S> hs) const noexcept
{
constexpr int Kind = Options::contains(kind_2) ? 2 : 1;
if constexpr(eve::integral_scalar_value<N>)
{
if constexpr(Options::contains(eve::spherical))
{
if (eve::is_ltz(n))
{
using u_t = eve::underlying_type_t<Z>;
return _::sb_hr(u_t(n), z, hs);
}
else
return _::sb_hn<Kind>(n, z, hs);
}
else
return _::cb_hn<Kind>(n, z, hs);
}
else
{
if constexpr(Options::contains(eve::spherical))
return _::sb_hr<Kind>(n, z, hs);
else
return _::cb_hr<Kind>(n, z, hs);
}
}

KYOSU_CALLABLE_OBJECT(bessel_h_t, bessel_h_);
};

//======================================================================================================================
//! @addtogroup functions
//! @{
//! @var bessel_h
//! @brief Computes the spherical or cylindrical Hankel functions,
//! extended to the complex plane and cayley_dickson algebras.
//!
//! @code
//! #include <kyosu/functions.hpp>
//! @endcode
//!
//! @groupheader{Callable Signatures}
//!
//! @code
//! namespace kyosu
//! {
//! template<eve;scalar_value N, kyosu::concepts::cayley_dickson T> constexpr auto bessel_h(N n, T z) noexcept; //1
//! template<eve;scalar_value N, kyosu::concepts::real T> constexpr T bessel_h(N n, T z) noexcept; //1
//! template<eve;scalar_value N, kyosu::concepts::complex T, size_t S> constexpr auto bessel_h(N n, T z, std::span<Z, S> chs) noexcept; //2
//! template<eve;scalar_value N, kyosu::concepts::real T, size_t S> constexpr T bessel_h(N n, T z, std::span<Z, S> chs) noexcept; //2
//!
//! template<eve;scalar_value N, kyosu::concepts::cayley_dickson T> constexpr auto bessel_h[spherical](N n, T z) noexcept; //3
//! template<eve;scalar_value N, kyosu::concepts::real T> constexpr T bessel_h[spherical](N n, T z) noexcept; //3
//! template<eve;scalar_value N, kyosu::concepts::complex T, size_t S> constexpr auto bessel_h[spherical](N n, T z, std::span<Z, S> shs) noexcept; //4
//! template<eve;scalar_value N, kyosu::concepts::real T, size_t S> constexpr T bessel_h[spherical](N n, T z, std::span<Z, S> shs) noexcept; //4
//!
//! template<eve;scalar_value N, kyosu::concepts::cayley_dickson T> constexpr auto bessel_h[kind_2](/*any previous overloads*/) noexcept; //5
//! template<eve;scalar_value N, kyosu::concepts::cayley_dickson T> constexpr auto bessel_h[kind_2][spherical](/*any previous overloads*/) noexcept; //6
//! }
//! @endcode
//!
//! **Parameters**
//!
//! * `n`: scalar order (integral or floating)
//! * `z`: Value to process.
//! * `shs`, `chs` : std::span of T
//!
//! **Return value**
//!
//! 1. returns \f$H_n\f$(z) (cylindrical first kind).
//! 2. Same as 1, but chs is filled by the lesser orders.
//! 3. returns \f$h_n\f$(z) (spherical).
//! 4. Same as 3, but shs is filled by the lesser orders.
//! 5. Same as 1., 2., 3., 4. but returns second kind Hankel functions values.
//! 6. Same as 1., 2., 3., 4. but returns second kind spherical Hankel functions values.
//!
//! @note
//! * Let \f$ i = \lfloor |n| \rfloor \f$ and \f$ j=|n|-i\f$. If \f$n\f$ is
//! positive the lesser order are \f$(\pm j, \pm(j+1), \dots, \pm(j+i)) \f$
//! with \f$+\f$ sign if \f$n\f$ is positive and \f$-\f$ sign if \f$n\f$ is negative.
//! * The span parameters are filled according the minimum of their allocated size and \f$i\f$.
//! * `cylindical` and `kind_1` options can be used and their results are identical to the regular call (no option).
//!
//! @groupheader{External references}
//! * [Wolfram MathWorld: Hankel Function](https://mathworld.wolfram.com/HankelFunction.html)
//! * [Wikipedia: Bessel function](https://en.wikipedia.org/wiki/Bessel_function)
//! * [DLMF: Bessel functions](https://dlmf.nist.gov/10.2)
//!
//! @groupheader{Example}
//! @godbolt{doc/bessel_h.cpp}
//======================================================================================================================
inline constexpr auto bessel_h = eve::functor<bessel_h_t>;
//======================================================================================================================
//! @}
//======================================================================================================================
}
Loading

0 comments on commit 0189a47

Please sign in to comment.