Skip to content

Commit

Permalink
Simplify bessel functions interface
Browse files Browse the repository at this point in the history
  • Loading branch information
jtlap authored Dec 21, 2024
1 parent 7a90989 commit 7815fc0
Show file tree
Hide file tree
Showing 98 changed files with 1,741 additions and 3,919 deletions.
22 changes: 4 additions & 18 deletions include/eve/module/bessel.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,21 +23,7 @@
#include <eve/module/bessel/regular/airy_ai.hpp>
#include <eve/module/bessel/regular/airy_bi.hpp>
#include <eve/module/bessel/regular/airy.hpp>
#include <eve/module/bessel/regular/cyl_bessel_i0.hpp>
#include <eve/module/bessel/regular/cyl_bessel_i1.hpp>
#include <eve/module/bessel/regular/cyl_bessel_in.hpp>
#include <eve/module/bessel/regular/cyl_bessel_j0.hpp>
#include <eve/module/bessel/regular/cyl_bessel_j1.hpp>
#include <eve/module/bessel/regular/cyl_bessel_jn.hpp>
#include <eve/module/bessel/regular/cyl_bessel_k0.hpp>
#include <eve/module/bessel/regular/cyl_bessel_k1.hpp>
#include <eve/module/bessel/regular/cyl_bessel_kn.hpp>
#include <eve/module/bessel/regular/cyl_bessel_y0.hpp>
#include <eve/module/bessel/regular/cyl_bessel_y1.hpp>
#include <eve/module/bessel/regular/cyl_bessel_yn.hpp>
#include <eve/module/bessel/regular/sph_bessel_j0.hpp>
#include <eve/module/bessel/regular/sph_bessel_j1.hpp>
#include <eve/module/bessel/regular/sph_bessel_jn.hpp>
#include <eve/module/bessel/regular/sph_bessel_y0.hpp>
#include <eve/module/bessel/regular/sph_bessel_y1.hpp>
#include <eve/module/bessel/regular/sph_bessel_yn.hpp>
#include <eve/module/bessel/regular/bessel_j.hpp>
#include <eve/module/bessel/regular/bessel_y.hpp>
#include <eve/module/bessel/regular/bessel_i.hpp>
#include <eve/module/bessel/regular/bessel_k.hpp>
76 changes: 76 additions & 0 deletions include/eve/module/bessel/detail/cbs_y.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
//==================================================================================================
/*
EVE - Expressive Vector Engine
Copyright : EVE Project Contributors
SPDX-License-Identifier: BSL-1.0
*/
//==================================================================================================
#pragma once
#include <eve/module/bessel/detail/kernel_bessel_y0.hpp>
#include <eve/module/bessel/detail/kernel_bessel_y1.hpp>
#include <eve/module/bessel/detail/kernel_bessel_y.hpp>
#include <eve/module/core.hpp>
#include <eve/module/math.hpp>
#include <eve/as_element.hpp>

namespace eve::detail
{
// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T cb_y(I nu, T x) noexcept
{
if constexpr(floating_value<I>)
{
if constexpr(scalar_floating_value<I>)
{
auto flintx = is_flint(nu);
return flintx ? kernel_bessel_y_int(nu, x) : kernel_bessel_y_flt(T(nu), x);
}
else if constexpr( scalar_value<T> )
{
auto flintx = is_flint(nu);
return flintx ? kernel_bessel_y_int(nu, x) : kernel_bessel_y_flt(nu, x);
}
else // simd
{
auto flint_nu = is_flint(nu);
if( eve::all(flint_nu) ) return kernel_bessel_y_int(nu, x);
else if( eve::none(flint_nu) ) return kernel_bessel_y_flt(nu, x);
else
{
auto nu_int = if_else(flint_nu, nu, zero);
auto nu_flt = if_else(flint_nu, T(0.5), nu);
return if_else(flint_nu, kernel_bessel_y_int(nu_int, x), kernel_bessel_y_flt(nu_flt, x));
}
}
}
else constexpr(integral_value<I>)
{
if constexpr(simd_value<I> && scalar_value<T>)
{
using c_t = wide<T, cardinal_t<I>>;
return cb_yn(nu, c_t(x));
}
else if constexpr(scalar_value<I> && simd_value<T>)
{
using i_t = wide<I, cardinal_t<T>>;
return kernel_bessel_y_int(i_t(nu), x);
}
else if constexpr(simd_value<I> && simd_value<T>)
{
auto n = convert(nu, as_element(x));
return kernel_bessel_y_int(n, x);
}
else return kernel_bessel_y_int(T(nu), x);
}
}

// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T sb_y(I nu, T x) noexcept
{
using elt_t = element_type_t<T>;
if constexpr( integral_value<I> ) return sph_bessel_yn(convert(n, as<elt_t>()), x);
else return cyl_bessel_yn(n + half(as(x)), x) * rsqrt(2 * x * inv_pi(as(x)));
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -6,50 +6,62 @@
*/
//==================================================================================================
#pragma once

#include <eve/module/bessel/detail/kernel_bessel_i0.hpp>
#include <eve/module/bessel/detail/kernel_bessel_i1.hpp>
#include <eve/module/bessel/detail/kernel_bessel_i.hpp>
#include <eve/module/core.hpp>
#include <eve/module/math.hpp>
#include <eve/as_element.hpp>

namespace eve::detail
{

template<typename I, typename T, callable_options O> constexpr
EVE_FORCEINLINE as_wide_as_t<T, I>
cyl_bessel_in_(EVE_REQUIRES(cpu_), O const&, I nu, T x)
{
// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T cb_i(I nu, T x) noexcept
{
if constexpr( integral_value<I> )
{
nu = eve::abs(nu); //DLMF 10.27.1 I-n(z) = In(z),
if constexpr( simd_value<I> && scalar_value<T> )
{
using c_t = wide<T, cardinal_t<I>>;
return cyl_bessel_in(convert(nu, as(x)), c_t(x));
return cb_i(convert(nu, as(x)), c_t(x));
}
else if constexpr( simd_value<I> && simd_value<T> )
{
using elt_t = element_type_t<T>;
auto tnu = convert(nu, as(elt_t()));
return cyl_bessel_in(tnu, x);
return cb_i(tnu, x);
}
else if constexpr( integral_scalar_value<I> )
{
return cyl_bessel_in(T(nu), x);
return cb_i(T(nu), x);
}
}
else // I is floating
{
if constexpr( simd_value<I> && scalar_value<T> )
{
using c_t = wide<T, cardinal_t<I>>;
return cyl_bessel_in(nu, c_t(x));
return cb_i(nu, c_t(x));
}
else if constexpr( scalar_value<I> && simd_value<T> )
{
return cyl_bessel_in(T(nu), x);
return cb_i(T(nu), x);
}
else
{
return kernel_bessel_i(nu, x);
}
}
}

// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T sb_i(I n, T x) noexcept
{
using elt_t = element_type_t<T>;
if constexpr( integral_value<I> ) return sb_i(convert(n, as<elt_t>()), x);
else return cb_i(n + half(as(x)), x) * rsqrt(2 * x * inv_pi(as(x)));
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -6,17 +6,17 @@
*/
//==================================================================================================
#pragma once

#include <eve/module/bessel/detail/kernel_bessel_j0.hpp>
#include <eve/module/bessel/detail/kernel_bessel_j1.hpp>
#include <eve/module/bessel/detail/kernel_bessel_j.hpp>
#include <eve/module/core.hpp>
#include <eve/module/math.hpp>
#include <eve/as_element.hpp>

namespace eve::detail
{
template<typename I, typename T, callable_options O> constexpr
EVE_FORCEINLINE as_wide_as_t<T, I>
cyl_bessel_jn_(EVE_REQUIRES(cpu_), O const&, I nu, T x)
template < typename I, typename T > constexpr EVE_FORCEINLINE
T cb_j(I nu, T x) noexcept
{
using w_t = as_wide_as_t<T, I>;
if constexpr(floating_value<I>)
Expand Down Expand Up @@ -49,9 +49,20 @@ namespace eve::detail
}
else // nu is integral
{
if constexpr(simd_value<I> && scalar_value<T>) return cyl_bessel_jn(nu, w_t(x));
if constexpr(simd_value<I> && scalar_value<T>) return cb_j(nu, w_t(x));
else if constexpr(scalar_value<I> && scalar_value<T>) return kernel_bessel_j_int(nu, x);
else if constexpr(simd_value<T>) return kernel_bessel_j_int(convert(nu, as_element(x)), x);
}
}

template < typename I, typename T > constexpr EVE_FORCEINLINE
T sb_j(I n, T x) noexcept
{
using elt_t = element_type_t<T>;
if constexpr( integral_value<I> ) return sb_j(convert(n, as<elt_t>()), x);
else return if_else(abs(x) < eps(as(x)) && is_gez(n),
if_else(is_eqz(n), one(as(x)), zero),
detail::cb_j(n + half(as(n)), x) * rsqrt(2 * x * inv_pi(as(x))));

}
}
79 changes: 79 additions & 0 deletions include/eve/module/bessel/detail/csb_k.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
//==================================================================================================
/*
EVE - Expressive Vector Engine
Copyright : EVE Project Contributors
SPDX-License-Identifier: BSL-1.0
*/
//==================================================================================================
#pragma once
#include <eve/module/bessel/detail/kernel_bessel_k0.hpp>
#include <eve/module/bessel/detail/kernel_bessel_k1.hpp>
#include <eve/module/bessel/detail/kernel_bessel_k.hpp>
#include <eve/module/core.hpp>
#include <eve/module/math.hpp>
#include <eve/as_element.hpp>

namespace eve::detail
{
// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T cb_k(I nu, T x) noexcept
{
if constexpr(floating_value<I>)
{
if constexpr(scalar_value<I> && simd_value<T>)
{
auto flintx = is_flint(nu);
return flintx ? kernel_bessel_k_int(nu, x) : kernel_bessel_k_flt(T(nu), x);
}
else if constexpr( scalar_value<T> )
{
auto flintx = is_flint(nu);
return flintx ? kernel_bessel_k_int(nu, x) : kernel_bessel_k_flt(nu, x);
}
else // simd
{
auto flint_nu = is_flint(nu);
if( eve::all(flint_nu) ) return kernel_bessel_k_int(nu, x);
else if( eve::none(flint_nu) ) return kernel_bessel_k_flt(nu, x);
else
{
auto nu_int = if_else(flint_nu, nu, zero);
auto nu_flt = if_else(flint_nu, T(0.5), nu);
return if_else(flint_nu, kernel_bessel_k_int(nu_int, x), kernel_bessel_k_flt(nu_flt, x));
}
}
}
else //if constexpr(integral_value<I>)
{
if constexpr(simd_value<I> && scalar_value<T>)
{
using c_t = wide<T, cardinal_t<I>>;
return cyl_bessel_kn(nu, c_t(x));
}
else if constexpr(scalar_value<I> && scalar_value<T>)
{
return kernel_bessel_k_int(nu, x);
}
else if constexpr(scalar_value<I> && simd_value<T>)
{
using i_t = wide<I, cardinal_t<T>>;
return kernel_bessel_k_int(i_t(nu), x);
}
else if constexpr(simd_value<I> && simd_value<T>)
{
auto n = convert(nu, as_element(x));
return kernel_bessel_k_int(n, x);
}
}
}

// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T sb_k(I n, T x) noexcept
{
using elt_t = element_type_t<T>;
if constexpr( integral_value<I> ) return sb_k(convert(n, as<elt_t>()), x);
else return cb_k(n + half(as(x)), x) * rsqrt(2 * x * inv_pi(as(x)));
}
}
76 changes: 76 additions & 0 deletions include/eve/module/bessel/detail/csb_y.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
//==================================================================================================
/*
EVE - Expressive Vector Engine
Copyright : EVE Project Contributors
SPDX-License-Identifier: BSL-1.0
*/
//==================================================================================================
#pragma once
#include <eve/module/bessel/detail/kernel_bessel_y0.hpp>
#include <eve/module/bessel/detail/kernel_bessel_y1.hpp>
#include <eve/module/bessel/detail/kernel_bessel_y.hpp>
#include <eve/module/core.hpp>
#include <eve/module/math.hpp>
#include <eve/as_element.hpp>

namespace eve::detail
{
// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T cb_y(I nu, T x) noexcept
{
if constexpr(floating_value<I>)
{
if constexpr(floating_scalar_value<I>)
{
auto flintx = is_flint(nu);
return flintx ? kernel_bessel_y_int(nu, x) : kernel_bessel_y_flt(T(nu), x);
}
else if constexpr( scalar_value<T> )
{
auto flintx = is_flint(nu);
return flintx ? kernel_bessel_y_int(nu, x) : kernel_bessel_y_flt(nu, x);
}
else // simd
{
auto flint_nu = is_flint(nu);
if( eve::all(flint_nu) ) return kernel_bessel_y_int(nu, x);
else if( eve::none(flint_nu) ) return kernel_bessel_y_flt(nu, x);
else
{
auto nu_int = if_else(flint_nu, nu, zero);
auto nu_flt = if_else(flint_nu, T(0.5), nu);
return if_else(flint_nu, kernel_bessel_y_int(nu_int, x), kernel_bessel_y_flt(nu_flt, x));
}
}
}
else //if constexpr(integral_value<I>)
{
if constexpr(simd_value<I> && scalar_value<T>)
{
using c_t = wide<T, cardinal_t<I>>;
return cb_yn(nu, c_t(x));
}
else if constexpr(scalar_value<I> && simd_value<T>)
{
using i_t = wide<I, cardinal_t<T>>;
return kernel_bessel_y_int(i_t(nu), x);
}
else if constexpr(simd_value<I> && simd_value<T>)
{
auto n = convert(nu, as_element(x));
return kernel_bessel_y_int(n, x);
}
else return kernel_bessel_y_int(T(nu), x);
}
}

// T is always floating
template < typename I, typename T > constexpr EVE_FORCEINLINE
T sb_y(I n, T x) noexcept
{
using elt_t = element_type_t<T>;
if constexpr( integral_value<I> ) return sb_y(convert(n, as<elt_t>()), x);
else return cb_y(n + half(as(x)), x) * rsqrt(2 * x * inv_pi(as(x)));
}
}
Loading

0 comments on commit 7815fc0

Please sign in to comment.