Spherical Bessel functions

doc/index.md

@@ -106,7 +106,7 @@ Operators
 ---------
 
 Operators (as said before) `+`, `-`, `*` and `/` can be used in infix form and can mix cayley-dickson values of
- different dimensinalities. Of course the biggest dimensionlity is recovered in the output.
+ different dimensionalities. Of course the biggest dimensionlity is recovered in the output.
 
 Prefix forms are also provided as `add`, `sub`, `multiply` and `div`. Also plus and minus for unary versions.
 
@@ -115,10 +115,9 @@ The left division sometimes necessary if the dimensionality is greater than 2 is
 Functions
 ---------
 
-Most **KYOSU** callables are usable with all cayley_dickson types. The exception being mainly special
-complex functions and rotation related quaternion usage.
+Most **KYOSU** callables are usable with all cayley_dickson types. The exceptions mainly being rotation related quaternion usage.
 
-@warning: **EVE** callables that correspond to mathematical functions that
+@warning  **EVE** callables that correspond to mathematical functions that
           are only defined on a proper part of the real axis as, for example, `acos` DOES NOT ever provide the same result
           if called in **EVE** or **KYOSU** context.
 
@@ -163,8 +162,14 @@ complex functions and rotation related quaternion usage.
    | [sqr](@ref kyosu::sqr )                   | [sqr_abs](@ref kyosu::sqr_abs )         | [sqrt](@ref kyosu::sqrt )             | [tan](@ref kyosu::tan )        | [tanpi](@ref kyosu::tanpi )          |
    | [tanh](@ref kyosu::tanh )                 | [to_polar](@ref kyosu::to_polar )       | [tgamma](@ref kyosu::tgamma )         | [trunc](@ref kyosu::trunc )    | [zeta](@ref kyosu::zeta )            |
 
+  * Bessel functions
 
-  * callables usable with quaternion complex and real only
+     They deserve a separate list to avoid verbosity:
+
+     * cylindrical Bessel functions: cyl_bessel_xxx are all defined for xxx being j0, j1, jn, y0, y1, yn, i0, i1, in, h1_0, h1_1, h1n, h2_0, h2_1, h2n, k1, k0, kn.
+     * spherical   Bessel functions: sph_bessel_xxx are all defined for xxx being j0, j1, jn, y0, y1, yn, i1_0, i1_1, i1n, i2_0, i2_1, i2n, h1_0, h1_1, h1n, h2_0, h2_1 h2n, k1, k0, kn.
+
+  * Callables usable with quaternion complex and real only
 
     These functions are related to \f$\mathbb{R}^3\f$ rotations
 

include/kyosu/details/cayleyify.hpp

@@ -23,6 +23,7 @@ auto cayley_extend_rev(auto f, auto z1, auto z2)
   return kyosu::real(c) + kyosu::ipart(c)*kyosu::sign(p);
 }
 
+
 auto cayley_extend2(auto f, auto z, auto... rs)
 {
   auto p = kyosu::pure(z);

include/kyosu/types/impl/detail/besseli0.hpp → include/kyosu/details/with_alloca.hpp

@@ -6,17 +6,16 @@
 */
 //======================================================================================================================
 #pragma once
-#include <kyosu/types/impl/detail/besselj0.hpp>
+#include <concepts>
 
 namespace kyosu::_
 {
 
-  //===-------------------------------------------------------------------------------------------
-  //  cyl_bessel_i0
-  //===-------------------------------------------------------------------------------------------
-  template<typename Z>
-  auto dispatch(eve::tag_of<kyosu::cyl_bessel_i0>, Z z) noexcept
+  template<typename T, std::invocable<T*> F>
+  decltype(auto) with_alloca(auto size, F f)
   {
-    return cyl_bessel_j0(complex(-ipart(z), real(z)));
+    T* p = (T*)(__builtin_alloca_with_align(size*sizeof(T), 8*alignof(T)));
+    return f(p);
   }
+
 }

include/kyosu/functions.hpp

@@ -163,8 +163,38 @@
 #include <kyosu/functions/sinpi.hpp>
 #include <kyosu/functions/sinpicospi.hpp>
 #include <kyosu/functions/sinh.hpp>
+#include <kyosu/functions/sinhc.hpp>
 #include <kyosu/functions/sinhcosh.hpp>
 #include <kyosu/functions/slerp.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_j0.hpp>
+#include <kyosu/functions/sph_bessel_j1.hpp>
+#include <kyosu/functions/sph_bessel_jn.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_y0.hpp>
+#include <kyosu/functions/sph_bessel_y1.hpp>
+#include <kyosu/functions/sph_bessel_yn.hpp>
+
+#include <kyosu/functions/sph_bessel_k0.hpp>
+#include <kyosu/functions/sph_bessel_k1.hpp>
+#include <kyosu/functions/sph_bessel_kn.hpp>
+
 #include <kyosu/functions/sqr.hpp>
 #include <kyosu/functions/sqr_abs.hpp>
 #include <kyosu/functions/sqrt.hpp>

include/kyosu/functions/cyl_bessel_i1.hpp

@@ -40,7 +40,7 @@ namespace kyosu
 //! @{
 //!   @var cyl_bessel_i1
 //!   @brief Computes  the modified Bessel function of the first kind,
-//!   \f$ I_1(x)= _iJ_1(ix) \f$ extended to the complex plane and cayley_dickson algebras.
+//!   \f$ I_1(x)= iJ_1(ix) \f$ extended to the complex plane and cayley_dickson algebras.
 //!
 //!   It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 0\f$.
 //!

include/kyosu/functions/from_angle_axis.hpp

@@ -26,10 +26,6 @@ namespace kyosu::tags
     {
       auto fn = callable_from_angle_axis{};
       return fn(angle, axis, normalize);
-//       using e_t = decltype(axis[0]+angle);
-//       auto q =  quaternion(e_t(0), e_t(axis[0]), e_t(axis[1]), e_t(axis[2]));
-//       auto [s, c] = eve::sincos(angle*eve::half(eve::as(angle)));
-//       return c+s*q;
     }
 
     template<eve::floating_ordered_value V,  typename  U, bool n>

include/kyosu/functions/sinhc.hpp

@@ -0,0 +1,74 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright: KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#pragma once
+
+#include <kyosu/details/invoke.hpp>
+#include <eve/module/math.hpp>
+
+namespace kyosu::tags
+{
+  struct callable_sinhc: eve::elementwise
+  {
+    using callable_tag_type = callable_sinhc;
+
+    KYOSU_DEFERS_CALLABLE(sinhc_);
+
+    template<eve::floating_ordered_value T>
+    static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::sinhc(v); }
+
+    template<typename T>
+    KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target))
+    {
+      return eve::tag_invoke(*this, target);
+    }
+
+    template<typename... T>
+    eve::unsupported_call<callable_sinhc(T&&...)> operator()(T&&... x) const
+    requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete;
+  };
+}
+
+namespace kyosu
+{
+//======================================================================================================================
+//! @addtogroup functions
+//! @{
+//!   @var sinhc
+//!   @brief Computes the hyperbolic sine cardinal of the argument.
+//!
+//!   **Defined in Header**
+//!
+//!   @code
+//!   #include <kyosu/functions.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace kyosu
+//!   {
+//!      template<kyosu::concepts::cayley_dickson T> constexpr auto sinhc(T z) noexcept;
+//!      template<eve::floating_ordered_value T>     constexpr auto sinhc(T z) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `z`: Value to process.
+//!
+//!   **Return value**
+//!
+//!     Returns the hyperbolic sine cardinal of the argument.
+//!
+//!  @groupheader{Example}
+//!
+//!  @godbolt{doc/sinhc.cpp}
+//! @}
+//======================================================================================================================
+inline constexpr tags::callable_sinhc sinhc = {};
+}

include/kyosu/functions/sph_bessel_h1_0.hpp

@@ -0,0 +1,77 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright: KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#pragma once
+
+#include <kyosu/details/invoke.hpp>
+#include <eve/module/bessel.hpp>
+
+namespace kyosu::tags
+{
+  struct callable_sph_bessel_h1_0: eve::elementwise
+  {
+    using callable_tag_type = callable_sph_bessel_h1_0;
+
+    KYOSU_DEFERS_CALLABLE(sph_bessel_h1_0_);
+
+    template<eve::floating_ordered_value T>
+    static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept
+    {
+      return complex(eve::sph_bessel_j0(z), eve::sph_bessel_y0(z));
+    }
+
+    template<typename T>
+    KYOSU_FORCEINLINE auto operator()(T const& target1) const noexcept
+    -> decltype(eve::tag_invoke(*this, target1))
+    {
+      return eve::tag_invoke(*this, target1);
+    }
+
+    template<typename... T>
+    eve::unsupported_call<callable_sph_bessel_h1_0(T&&...)> operator()(T&&... x) const
+    requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete;
+  };
+}
+
+namespace kyosu
+{
+//======================================================================================================================
+//! @addtogroup functions
+//! @{
+//!   @var sph_bessel_h1_0
+//!   @brief Computes the spherical Bessel/hankel functions of the third kind,
+//!   \f$ h_0^{(1)}(z) = j_1(z)+iy_1(z)\f$.
+//!
+//!   @code
+//!   #include <kyosu/functions.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace kyosu
+//!   {
+//!      template<kyosu::concepts::cayley_dickson T> constexpr auto sph_bessel_h1_0(int n, T z) noexcept;
+//!      template<eve::floating_ordered_value T>     constexpr auto sph_bessel_h1_0(int n, T z) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `z`: Value to process.
+//!
+//!   **Return value**
+//!
+//!     * returns  \f$h_0^{(1)}(z)\f$.
+//!
+//!  @groupheader{Example}
+//!
+//!  @godbolt{doc/sph_bessel_h1_0.cpp}
+//! @}
+//======================================================================================================================
+inline constexpr tags::callable_sph_bessel_h1_0 sph_bessel_h1_0 = {};
+}

include/kyosu/functions/sph_bessel_h1_1.hpp

@@ -0,0 +1,77 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright: KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#pragma once
+
+#include <kyosu/details/invoke.hpp>
+#include <eve/module/bessel.hpp>
+
+namespace kyosu::tags
+{
+  struct callable_sph_bessel_h1_1: eve::elementwise
+  {
+    using callable_tag_type = callable_sph_bessel_h1_1;
+
+    KYOSU_DEFERS_CALLABLE(sph_bessel_h1_1_);
+
+    template<eve::floating_ordered_value T>
+    static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept
+    {
+      return complex(eve::sph_bessel_j1(z), eve::sph_bessel_y1(z));
+    }
+
+    template<typename T>
+    KYOSU_FORCEINLINE auto operator()(T const& target1) const noexcept
+    -> decltype(eve::tag_invoke(*this, target1))
+    {
+      return eve::tag_invoke(*this, target1);
+    }
+
+    template<typename... T>
+    eve::unsupported_call<callable_sph_bessel_h1_1(T&&...)> operator()(T&&... x) const
+    requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete;
+  };
+}
+
+namespace kyosu
+{
+//======================================================================================================================
+//! @addtogroup functions
+//! @{
+//!   @var sph_bessel_h1_1
+//!   @brief Computes the spherical Bessel/hankel functions of the third kind,
+//!   \f$ h_1^{(1)}(z) = j_1(z)+iy_1(z)\f$.
+//!
+//!   @code
+//!   #include <kyosu/functions.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace kyosu
+//!   {
+//!      template<kyosu::concepts::cayley_dickson T> constexpr auto sph_bessel_h1_1(int n, T z) noexcept;
+//!      template<eve::floating_ordered_value T>     constexpr auto sph_bessel_h1_1(int n, T z) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `z`: Value to process.
+//!
+//!   **Return value**
+//!
+//!     * returns  \f$h_1^{(1)}(z)\f$.
+//!
+//!  @groupheader{Example}
+//!
+//!  @godbolt{doc/sph_bessel_h1_1.cpp}
+//! @}
+//======================================================================================================================
+inline constexpr tags::callable_sph_bessel_h1_1 sph_bessel_h1_1 = {};
+}