diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp
index 8d6c61b5..355dfcce 100644
--- a/include/kyosu/functions.hpp
+++ b/include/kyosu/functions.hpp
@@ -17,6 +17,11 @@
 #include <kyosu/functions/acospi.hpp>
 #include <kyosu/functions/acot.hpp>
 #include <kyosu/functions/acotpi.hpp>
+
+#include <kyosu/functions/airy.hpp>
+#include <kyosu/functions/airy_ai.hpp>
+#include <kyosu/functions/airy_bi.hpp>
+
 #include <kyosu/functions/atanpi.hpp>
 #include <kyosu/functions/acoth.hpp>
 #include <kyosu/functions/arg.hpp>
diff --git a/include/kyosu/functions/airy.hpp b/include/kyosu/functions/airy.hpp
new file mode 100644
index 00000000..fc777f30
--- /dev/null
+++ b/include/kyosu/functions/airy.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_airy: eve::elementwise
+  {
+    using callable_tag_type = callable_airy;
+
+    KYOSU_DEFERS_CALLABLE(airy_);
+
+    template<eve::floating_ordered_value T>
+    static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept
+    {
+      return eve::airy(v);
+    }
+
+    template<typename T>
+    KYOSU_FORCEINLINE auto operator()(T const & target0) const noexcept
+    -> decltype(eve::tag_invoke(*this, target0))
+    {
+      return eve::tag_invoke(*this, target0);
+    }
+
+    template<typename... T>
+    eve::unsupported_call<callable_airy(T&&...)> operator()(T&&... x) const
+    requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete;
+  };
+}
+
+namespace kyosu
+{
+//======================================================================================================================
+//! @addtogroup functions
+//! @{
+//!   @var airy
+//!   @brief Computes simultaneously the airy functions \f$Ai\f$ and \f$Bi\f$.
+//!
+//!   @code
+//!   #include <kyosu/functions.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace kyosu
+//!   {
+//!      template<eve::floating_ordered_value T> constexpr auto airy(T z) noexcept;
+//!
+//!      template<kyosu::concepts::cayley_dickson T>  constexpr auto airy(T z) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `z`: Value to process.
+//!
+//!   **Return value**
+//!
+//!     * returns  a kumi pair containing \f$\{Ai(z)\}f$ and \f$\{Bi(z)\}f$.
+//!
+//!  @groupheader{Example}
+//!
+//!  @godbolt{doc/airy.cpp}
+//! @}
+//======================================================================================================================
+inline constexpr tags::callable_airy airy = {};
+}
diff --git a/include/kyosu/functions/airy_ai.hpp b/include/kyosu/functions/airy_ai.hpp
new file mode 100644
index 00000000..d4d4a630
--- /dev/null
+++ b/include/kyosu/functions/airy_ai.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_airy_ai: eve::elementwise
+  {
+    using callable_tag_type = callable_airy_ai;
+
+    KYOSU_DEFERS_CALLABLE(airy_ai_);
+
+    template<eve::floating_ordered_value T>
+    static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept
+    {
+      return eve::airy_ai(v);
+    }
+
+    template<typename T>
+    KYOSU_FORCEINLINE auto operator()(T const & target0) const noexcept
+    -> decltype(eve::tag_invoke(*this, target0))
+    {
+      return eve::tag_invoke(*this, target0);
+    }
+
+    template<typename... T>
+    eve::unsupported_call<callable_airy_ai(T&&...)> operator()(T&&... x) const
+    requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete;
+  };
+}
+
+namespace kyosu
+{
+//======================================================================================================================
+//! @addtogroup functions
+//! @{
+//!   @var airy_ai
+//!   @brief Computes the airy function \f$Ai\f$,
+//!
+//!   @code
+//!   #include <kyosu/functions.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace kyosu
+//!   {
+//!      template<eve::floating_ordered_value T> constexpr auto airy_ai(T z) noexcept;
+//!
+//!      template<kyosu::concepts::cayley_dickson T> constexpr auto airy_ai(T z) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `z`: Value to process.
+//!
+//!   **Return value**
+//!
+//!     * returns  \f$\{Ai(z)\}f$.
+//!
+//!  @groupheader{Example}
+//!
+//!  @godbolt{doc/airy_ai.cpp}
+//! @}
+//======================================================================================================================
+inline constexpr tags::callable_airy_ai airy_ai = {};
+}
diff --git a/include/kyosu/functions/airy_bi.hpp b/include/kyosu/functions/airy_bi.hpp
new file mode 100644
index 00000000..d3575bd1
--- /dev/null
+++ b/include/kyosu/functions/airy_bi.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_airy_bi: eve::elementwise
+  {
+    using callable_tag_type = callable_airy_bi;
+
+    KYOSU_DEFERS_CALLABLE(airy_bi_);
+
+    template<eve::floating_ordered_value T>
+    static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept
+    {
+      return eve::airy_bi(v);
+    }
+
+    template<typename T>
+    KYOSU_FORCEINLINE auto operator()(T const & target0) const noexcept
+    -> decltype(eve::tag_invoke(*this, target0))
+    {
+      return eve::tag_invoke(*this, target0);
+    }
+
+    template<typename... T>
+    eve::unsupported_call<callable_airy_bi(T&&...)> operator()(T&&... x) const
+    requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete;
+  };
+}
+
+namespace kyosu
+{
+//======================================================================================================================
+//! @addtogroup functions
+//! @{
+//!   @var airy_bi
+//!   @brief Computes the airy function \f$Bi\f$,
+//!
+//!   @code
+//!   #include <kyosu/functions.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace kyosu
+//!   {
+//!      template<eve::floating_ordered_value T> constexpr auto airy_bi(T z) noexcept;
+//!
+//!      template<kyosu::concepts::cayley_dickson T> constexpr auto airy_bi(T z) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `z`: Value to process.
+//!
+//!   **Return value**
+//!
+//!     * returns  \f$\{Bi(z)\}f$.
+//!
+//!  @groupheader{Example}
+//!
+//!  @godbolt{doc/airy_bi.cpp}
+//! @}
+//======================================================================================================================
+inline constexpr tags::callable_airy_bi airy_bi = {};
+}
diff --git a/include/kyosu/functions/cyl_bessel_y.hpp b/include/kyosu/functions/cyl_bessel_y.hpp
index 2b4c0e0e..33df0eea 100644
--- a/include/kyosu/functions/cyl_bessel_y.hpp
+++ b/include/kyosu/functions/cyl_bessel_y.hpp
@@ -25,6 +25,13 @@ namespace kyosu::tags
       return fnu(nu, complex(v));
     }
 
+    template<eve::floating_scalar_value N, eve::floating_ordered_value T, typename R>
+    static KYOSU_FORCEINLINE auto deferred_call(auto, N nu, T const& v, R& ys) noexcept
+    {
+      auto fnu = callable_cyl_bessel_y{};
+      return fnu(nu, complex(v), ys);
+    }
+
     template<typename N, typename T>
     KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept
     -> decltype(eve::tag_invoke(*this, target0, target1))
@@ -32,6 +39,13 @@ namespace kyosu::tags
       return eve::tag_invoke(*this, target0, target1);
     }
 
+    template<typename N, typename T, typename R>
+    KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1, R & out0) const noexcept
+    -> decltype(eve::tag_invoke(*this, target0, target1, out0))
+    {
+      return eve::tag_invoke(*this, target0, target1, out0);
+    }
+
     template<typename... T>
     eve::unsupported_call<callable_cyl_bessel_y(T&&...)> operator()(T&&... x) const
     requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete;
diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp
index 2d29e9b2..71866a6f 100644
--- a/include/kyosu/types/impl/bessel.hpp
+++ b/include/kyosu/types/impl/bessel.hpp
@@ -111,3 +111,5 @@
 #include <kyosu/types/impl/besselr/cb_jyr.hpp>
 #include <kyosu/types/impl/besselr/cb_h.hpp>
 #include <kyosu/types/impl/besselr/cb_ikr.hpp>
+// airy functions
+#include <kyosu/types/impl/besselr/airy.hpp>
diff --git a/include/kyosu/types/impl/besselr/airy.hpp b/include/kyosu/types/impl/besselr/airy.hpp
new file mode 100644
index 00000000..60001def
--- /dev/null
+++ b/include/kyosu/types/impl/besselr/airy.hpp
@@ -0,0 +1,189 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright : KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#pragma once
+
+namespace kyosu::_
+{
+  /////////////////////////////////
+  // contains implementations of
+  // airy_ai
+  // airy_bi
+  // airy
+  /////////////////////////////////
+  // utilities
+  // zet
+  // ai
+  // bi
+  // aibi
+  /////////////////////////////////
+  // implementation is done using j and y functions
+
+  //===-------------------------------------------------------------------------------------------
+  // zet
+  //===-------------------------------------------------------------------------------------------
+  template<typename Z> KYOSU_FORCEINLINE
+  auto zet(Z z) noexcept
+  {
+    using u_t = eve::underlying_type_t<Z>;
+     auto sqz = sqrt(z);
+//    auto zeta = 2*sqz*z/3;
+     auto zeta = (pow(z, u_t(1.5))*2)/3;
+    return kumi::tuple{sqz/eve::sqrt_3(as(z)), zeta};
+  }
+
+  //===-------------------------------------------------------------------------------------------
+  // ai
+  //===-------------------------------------------------------------------------------------------
+  template<typename Z> KYOSU_FORCEINLINE
+  auto ai(Z z) noexcept
+  {
+    using u_t =  eve::underlying_type_t<Z>;
+    auto third =  eve::third(as<u_t>());
+    auto invpi =  eve::inv_pi(as<u_t>());
+    auto br_im_eq_0 = [](auto z){
+      return Z(eve::airy_ai(real(z)));
+    };
+
+    auto br_re_lt_0 = [third, invpi](auto z){
+      auto izgt0 = eve::is_gez(imag(z));
+      z = if_else(izgt0, z, conj(z));
+      auto zet = 2*pow(z, u_t(1.5) )/3;
+      auto argzetgtpio2 =  arg(zet) > eve::pio_2(eve::as< u_t>());
+      auto sqtmz = sqrt(-z/3);
+      auto r = mulmi(invpi*sqtmz*cyl_bessel_k(third, zet));
+      r =  if_else(izgt0, r, conj(r));
+      return if_else(argzetgtpio2, -r, r);
+    };
+    auto br_re_gt_0 = [third, invpi](auto z){
+      auto zet = 2*pow(z, u_t(1.5) )/3;
+      std::array<Z, 1> ks;
+      auto r = invpi*sqrt(z/3)*cyl_bessel_k(third, zet, ks);
+      return r;
+    };
+
+    auto notdone = eve::true_(as<Z>());
+    Z r = eve::nan(as<Z>());
+    auto re = real(z);
+    if( eve::any(notdone) )
+    {
+      notdone = next_interval(br_im_eq_0, notdone, is_real(z), r, z);
+      if( eve::any(notdone) )
+      {
+        notdone = next_interval(br_re_gt_0, notdone, eve::is_gtz(re), r, z);
+        if( eve::any(notdone) )
+        {
+          notdone = last_interval(br_re_lt_0, notdone, r, z);
+        }
+      }
+    }
+    return r;
+  }
+
+  //===-------------------------------------------------------------------------------------------
+  // bi
+  //===-------------------------------------------------------------------------------------------
+  template<typename Z> KYOSU_FORCEINLINE
+  auto bi(Z z) noexcept
+  {
+    using u_t =  eve::underlying_type_t<Z>;
+    auto third =  eve::third(as<u_t>());
+    auto invpi =  eve::inv_pi(as<u_t>());
+    auto br_im_eq_0 = [](auto z){
+      return Z(eve::airy_bi(real(z)));
+    };
+
+    auto br_re_lt_0 = [third, invpi](auto z){
+      auto zet = 2*pow(-z, u_t(1.5) )/3;
+      return (kyosu::cyl_bessel_j(-third, zet)-kyosu::cyl_bessel_j(third, zet))*sqrt(-z/3) ;
+    };
+
+    auto br_re_gt_0 = [third, invpi](auto z){
+      auto zet = 2*pow(z, u_t(1.5) )/3;
+      return sqrt(z/3)*(cyl_bessel_i(third, zet)+cyl_bessel_i(-third, zet));
+    };
+
+    auto notdone = eve::true_(as<Z>());
+    Z r = eve::nan(as<Z>());
+    auto re = real(z);
+    if( eve::any(notdone) )
+    {
+      notdone = next_interval(br_im_eq_0, notdone, is_real(z), r, z);
+      if( eve::any(notdone) )
+      {
+        notdone = next_interval(br_re_gt_0, notdone, eve::is_gtz(re), r, z);
+        if( eve::any(notdone) )
+        {
+          notdone = last_interval(br_re_lt_0, notdone, r, z);
+        }
+      }
+    }
+    return r;
+  }
+
+  //===-------------------------------------------------------------------------------------------
+  // aibi
+  //===-------------------------------------------------------------------------------------------
+  template<typename Z> KYOSU_FORCEINLINE
+  auto aibi(Z z) noexcept
+  {
+    using u_t = eve::underlying_type_t<Z>;
+    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);
+
+    return kumi::tuple{ invpi(as<u_t>())*sqzo3*(im-ip),  sqzo3*(ip+im)};
+  }
+
+  //===-------------------------------------------------------------------------------------------
+  //  airy_ai
+  //===-------------------------------------------------------------------------------------------
+  template<typename Z> KYOSU_FORCEINLINE
+  auto dispatch(eve::tag_of<kyosu::airy_ai>, Z z) noexcept
+  {
+    if constexpr(kyosu::concepts::complex<Z> )
+    {
+      return ai(z);
+    }
+    else
+    {
+      return cayley_extend(airy_ai, z);
+    }
+  }
+
+  //===-------------------------------------------------------------------------------------------
+  //  airy_bi
+  //===-------------------------------------------------------------------------------------------
+  template<typename Z> KYOSU_FORCEINLINE
+  auto dispatch(eve::tag_of<kyosu::airy_bi>, Z z) noexcept
+  {
+    if constexpr(concepts::complex<Z> )
+    {
+      return bi(z);
+    }
+    else
+    {
+      return cayley_extend(airy_bi, z);
+    }
+  }
+
+  //===-------------------------------------------------------------------------------------------
+  //  airy
+  //===-------------------------------------------------------------------------------------------
+  template<typename Z> KYOSU_FORCEINLINE
+  auto dispatch(eve::tag_of<kyosu::airy>, Z z) noexcept
+  {
+    if constexpr(concepts::complex<Z> )
+    {
+      return aibi(z);
+    }
+    else
+    {
+      return cayley_extend2(airy, z);
+    }
+  }
+}
diff --git a/include/kyosu/types/impl/besselr/cb_jyr.hpp b/include/kyosu/types/impl/besselr/cb_jyr.hpp
index 0494dec2..e38f3cac 100644
--- a/include/kyosu/types/impl/besselr/cb_jyr.hpp
+++ b/include/kyosu/types/impl/besselr/cb_jyr.hpp
@@ -40,12 +40,13 @@ namespace kyosu::_
         cjv[1] = cjv1;
         cyv[1] = cyv1;
       }
-      if(n == 0)  {
-        return kumi::tuple{cjv0, cyv0};
+      if(n == 0)
+      {
+        return kumi::tuple{cjv[0], cyv[0]};
       }
       else if (n == 1)
       {
-        return kumi::tuple{cjv1, cyv1};
+        return kumi::tuple{cjv[1], cyv[1]};
       }
       else if (n >= 2)
       {
@@ -126,16 +127,24 @@ namespace kyosu::_
     else
     {
       // v is not flint and is less than 0;
+      // we use:
+      // J_{-v] = J_v(z) cos(v pi) - Y_v *sin(v pi)
+      // y_{-v] = J_v(z) sin(v pi) + Y_v *cos(v pi)
       cb_jy(-v, z, cjv, cyv);
-      for(int i=0; i < an ; ++i)
+      auto v0 = frac(-v);
+      auto [s, c] = sinpicospi(v0);
+      for(int i=0; i <= an ; ++i)
       {
-        auto [s, c] = sinpicospi(v);
-        auto a = cjv[i]*s;
-        auto b = cyv[i]*c;
-        cjv[i] = a-b;
-        cyv[i] = a+b;
+        auto aa = cjv[i]*c;
+        auto bb = cyv[i]*s;
+        auto cc = cjv[i]*s;
+        auto dd = cyv[i]*c;
+        cjv[i] = aa-bb; //mulmi(a-b);
+        cyv[i] = cc+dd;
+        s = -s;
+        c = -c;
       }
-      return kumi::tuple{cjv[n], cyv[n]};
+      return kumi::tuple{cjv[an], cyv[an]};
     }
   }
 
@@ -154,10 +163,10 @@ namespace kyosu::_
   template<eve::floating_scalar_value N, kyosu::concepts::complex Z, typename R>
   auto dispatch(eve::tag_of<kyosu::cyl_bessel_j>, N v, Z z, R & js) noexcept
   {
-    auto n = int(abs(v))+1;
-    auto doit = [v, z, &js](auto ys){
+    auto n = int(abs(v));
+    auto doit = [n, v, z, &js](auto ys){
       auto [jvn, _] = cb_jy(v, z, js, ys);
-      return jvn;
+      return js[n];
     };
     return with_alloca<Z>(n+1, doit);
   }
@@ -168,12 +177,12 @@ namespace kyosu::_
   template<eve::floating_scalar_value N, typename Z>
   auto dispatch(eve::tag_of<kyosu::cyl_bessel_j>, N v, Z z) noexcept
   {
-    auto n = int(abs(v))+1;
+    auto n = int(abs(v));
     if constexpr(concepts::complex<Z> )
     {
-      auto doit = [v, z](auto js, auto ys){
+      auto doit = [n, v, z](auto js, auto ys){
         auto [jv, _] = cb_jy(v, z, js, ys);
-        return jv;
+        return js[n];
       };
       return with_alloca<Z>(n+1, doit);
     }
@@ -203,7 +212,7 @@ namespace kyosu::_
   template<eve::floating_scalar_value N, typename Z>
   auto dispatch(eve::tag_of<kyosu::cyl_bessel_y>, N v, Z z) noexcept
   {
-    auto n = int(abs(v))+1;
+    auto n = int(abs(v));
     if constexpr(concepts::complex<Z> )
     {
       auto doit = [v, z](auto js, auto ys){
diff --git a/test/doc/cyl_bessel_j.cpp b/test/doc/cyl_bessel_j.cpp
index 4200b279..0d86ba39 100644
--- a/test/doc/cyl_bessel_j.cpp
+++ b/test/doc/cyl_bessel_j.cpp
@@ -7,13 +7,18 @@ int main()
   std::cout.precision(16);
   using w_t = eve::wide<double, eve::fixed<2>>;
   auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5));
-  auto v = 3.3;
+  auto v = -(2+1/3.0);
+  int nb = int(eve::abs(v)+1);
   std::cout << "z                    " << z << std::endl;
-  std::vector<decltype(z)> js(4);
+  std::vector<decltype(z)> js(nb);
   kyosu::cyl_bessel_j(v, z, js);
-  for(int n=0; n <= 3; ++n)
+  auto inc = eve::sign(v);
+  auto v0 =  eve::frac(v);
+  for(int n=0; n < nb; ++n)
   {
-    std::cout << "js[" << n << "] = " << js[n] << std::endl;
+   std::cout << "js[" << n << "]                    = " << js[n] << std::endl;
+    std::cout << "j(" << v0 << ",  z) " << kyosu::cyl_bessel_j(v0, z) << std::endl;
+    v0 += inc;
   }
   return 0;
 }
diff --git a/test/unit/function/airy_ai.cpp b/test/unit/function/airy_ai.cpp
new file mode 100644
index 00000000..98d10292
--- /dev/null
+++ b/test/unit/function/airy_ai.cpp
@@ -0,0 +1,55 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright : KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#include <kyosu/kyosu.hpp>
+#include <test.hpp>
+
+TTS_CASE_WITH ( "Check kyosu::airy_ai over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+<typename T>(T )
+{
+  using e_t = eve::element_type_t<T>;
+  if constexpr(sizeof(e_t) == 8)
+  {
+    std::array<T, 16> re{
+      -5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00
+        };
+    std::array<T, 16> im{
+      8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00
+        };
+    std::array<T, 16> reres{
+      1.6763717755962929e+07, 1.0671669307692060e+06, 5.1770135828865488e+04, 2.7209127311234188e+03
+        , 1.7464957947925365e+02, 1.5043848907202149e+01, 1.9864751505771967e+00, 5.2597973488369920e-01
+        , 2.7880648195500490e-01, 6.0458308371838042e-02, -1.0036228447036250e-01, 2.4722152686112525e-02
+        , 8.1141306614663667e-03, -1.1116711994387154e-02, 6.2362124668067535e-03, -1.1845517847442273e-03
+        };
+    std::array<T, 16> imres{
+      5.7419608040318936e+07, -1.2409242416987466e+06, 2.6427659633621028e+04, -3.4216030566322274e+02
+        , -3.0804431341868103e+01, 6.6862293791164378e+00, -1.3118365573923054e+00, 3.2173225694177510e-01
+        , 0.0000000000000000e+00, -1.5188956587718155e-01, 2.1867029634299019e-02, 6.0030132776007987e-02
+        , 3.9632838829414818e-02, 2.2033585182583651e-02, 1.3313384345125566e-02, 8.4631437627202440e-03};
+
+    for(int i=0; i < 16; ++i)
+    {
+     auto c = kyosu::complex(re[i], im[i]);
+     auto res = kyosu::complex(reres[i], imres[i]);
+     TTS_RELATIVE_EQUAL(kyosu::airy_ai(c), res, 1.0e-3) << i <<  " <- " << c << "arg " << kyosu::arg(c) << '\n';
+     TTS_RELATIVE_EQUAL(kyosu::airy_ai(re[i]), kyosu::real(kyosu::airy_ai(kyosu::complex(re[i], e_t(0.0)))), 1.0e-4)<< i <<  " -> " << re[i] << '\n';
+     TTS_RELATIVE_EQUAL(kyosu::airy_ai(im[i]), kyosu::real(kyosu::airy_ai(kyosu::complex(im[i], e_t(0.0)))), 1.0e-4)<< i <<  " -> " << im[i] << '\n';
+     TTS_RELATIVE_EQUAL(kyosu::airy_ai(kyosu::conj(c)), kyosu::conj(res), 1.0e-3) << i <<  " <- " << c << "arg " << kyosu::arg(c) << '\n';
+
+    }
+  }
+};
diff --git a/test/unit/function/airy_bi.cpp b/test/unit/function/airy_bi.cpp
new file mode 100644
index 00000000..f346e4a5
--- /dev/null
+++ b/test/unit/function/airy_bi.cpp
@@ -0,0 +1,62 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright : KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#include <kyosu/kyosu.hpp>
+#include <test.hpp>
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+<typename T>(T )
+{
+  using e_t = eve::element_type_t<T>;
+  if constexpr(sizeof(e_t) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00, -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01, 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00, 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00, 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00, -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00, -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+
+    std::array<T, 16> reres{-5.7419608040318385e+07, -1.2409242416987633e+06, -2.6427659633595144e+04, -3.4216029692014837e+02, 3.0804691728218572e+01, 6.6903996777286583e+00, 1.3523669469741970e+00, 5.5870307892357129e-01, 7.5248558508731545e-01, 7.1665807338276855e-01, -7.3369695701857451e-01, -1.3293492209488911e-01, 1.1081359835209212e+00, -2.0805776844472024e+00, 3.1132412017028535e+00, -3.5952187284708805e+00};
+
+    std::array<T, 16> imres{1.6763717755961375e+07, -1.0671669307691315e+06, 5.1770135827841987e+04, -2.7209127089468357e+03, 1.7464925781992241e+02, -1.5040759898113270e+01, 1.9670200948712224e+00, -4.4241768244532709e-01, 0.0000000000000000e+00, 6.1988929040084473e-01, -4.7545485913890873e-01, -1.2208933006131721e+00, -1.3901622231770903e+00, -1.5528897868070297e+00, -2.4569086817625085e+00, -5.1871698800009600e+00};
+
+    for(int i=0; i < 16; ++i)
+    {
+      auto c = kyosu::complex(re[i], im[i]);
+      auto res = kyosu::complex(reres[i], imres[i]);
+      TTS_RELATIVE_EQUAL(kyosu::airy_bi(c), res, 1.0e-3) << i <<  " <- " << c << '\n';
+      TTS_RELATIVE_EQUAL(kyosu::airy_bi(re[i]), kyosu::real(kyosu::airy_bi(kyosu::complex(re[i], e_t(0.0)))), 1.0e-4)<< re[i] << '\n';
+      TTS_RELATIVE_EQUAL(kyosu::airy_bi(im[i]), kyosu::real(kyosu::airy_bi(kyosu::complex(im[i], e_t(0.0)))), 1.0e-4)<< im[i] << '\n';
+    }
+  }
+};
+
+// TTS_CASE_WITH ( "Check kyosu::airy_bi over real"
+//               , kyosu::real_types
+//               , tts::generate(tts::randoms(-10,10),
+//                               tts::randoms(-10,10)
+//                              )
+//               )
+// <typename T>(T a0, T a1)
+// {
+//   using u_t = eve::underlying_type_t<T>;
+//   if constexpr(sizeof(u_t) == 8)
+//   {
+//     auto c =  kyosu::complex(a0, a1);
+//     auto cb=  kyosu::conj(c);
+// //     auto cr=  kyosu::complex(a0);
+// //     auto ci=  kyosu::complex(T(0), a1);
+//     using kyosu::airy_bi;
+//     auto j0c = airy_bi(c);
+//     TTS_EQUAL(j0c, kyosu::conj(airy_bi(cb))) << "c = " << c << '\n';
+// //     TTS_EXPECT(eve::all(kyosu::is_real(cr)));
+// //     TTS_EXPECT(eve::all(kyosu::is_pure(ci)));
+// //     auto z =   kyosu::complex(T(0), T(0));
+// //     TTS_IEEE_EQUAL(airy_bi(z), kyosu::complex(airy_bi(T(0))));
+//   }
+// };
diff --git a/test/unit/function/cyl_bessel_i.cpp b/test/unit/function/cyl_bessel_i.cpp
index 33b35b11..f552918b 100644
--- a/test/unit/function/cyl_bessel_i.cpp
+++ b/test/unit/function/cyl_bessel_i.cpp
@@ -65,7 +65,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i over real"
       auto vv = T(0.3);
       for (n = 0; n <= N; ++n)
       {
-       std::cout << "n = "<< n << " k = "<< k << std::endl;
         auto c = kyosu::complex(re[k], im[k]);
         kyosu::cyl_bessel_i(v1, c, is_03);
         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
diff --git a/test/unit/function/cyl_bessel_j.cpp b/test/unit/function/cyl_bessel_j.cpp
new file mode 100644
index 00000000..01ba2710
--- /dev/null
+++ b/test/unit/function/cyl_bessel_j.cpp
@@ -0,0 +1,119 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright : KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#include <kyosu/kyosu.hpp>
+#include <test.hpp>
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+  <typename T>(T )
+{
+  if constexpr(sizeof(T) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+    constexpr int N = 0;
+    T v1 = 1.0/3;
+    std::cout << "v1 " << v1 << std::endl;
+    auto n = int(v1);
+    using vc_t = std::array<T, 16>;
+    using a_t = std::array<vc_t, N+1>;
+    a_t reresi_03{
+      vc_t{2.8090373106452773e+02, 1.7536493250871054e+01, -3.3400462158788272e+01, -2.3524142528283910e+01
+          , -9.3274540332656208e+00, -1.9656734427524503e+00, 4.3413771695853465e-01, 7.4917804310983560e-01
+          , 5.8501480583737409e-01, 1.0213243250718731e+00, 1.3227852389571668e+00, 1.9818890844578396e-01
+          , -6.0472153767481238e+00, -2.3374276615174153e+01, -5.1769656721627001e+01, -5.2626473557055462e+01
+          }
+
+    };
+    a_t imresi_03{
+      vc_t{-2.6513665958745401e+02, 1.5122859663021515e+02, -5.0648785179729913e+01, 7.0192828558358924e+00
+          , 4.0222943595848140e+00, -3.8995769510156886e+00, 1.9902620949926999e+00, -7.5756275269309670e-01
+          , 0.0000000000000000e+00, -1.1158707376004939e-01, 1.3330326497253218e+00, -4.1355677304829221e+00
+          , 7.6943469735731851e+00, -4.6838788665148172e+00, -2.8742293351241671e+01, 1.3952363114976131e+02
+          }
+    };
+
+    std::vector<kyosu::complex_t<T>> is_03(N+2);
+    for (int k = 0; k < 16; ++k)
+    {
+      auto c = kyosu::complex(re[k], im[k]);
+      kyosu::cyl_bessel_j(v1, c, is_03);
+      auto vv = T(N+1.0/3);
+      for (n = 0; n <= N; ++n)
+      {
+         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
+        TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_j(vv, c), 1.e-2) << "n " << n  << " k " << k  << " c "<< c <<'\n';
+      }
+    }
+  }
+};
+
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+  <typename T>(T )
+{
+  if constexpr(sizeof(T) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+    constexpr int N = 0;
+    T v1 = -1.0/3;
+    std::cout << "v1 " << v1 << std::endl;
+    auto n = int(v1);
+    using vc_t = std::array<T, 16>;
+    using a_t = std::array<vc_t, N+1>;
+    a_t reresi_03{
+      vc_t{-8.9163143837610050e+01, -1.2219942089268812e+02, -6.0563369116851803e+01, -1.7842181975937766e+01
+          , -1.1860080210885884e+00, 2.3789187887624883e+00, 1.9232059348460351e+00, 1.1368224659946016e+00
+          , 1.3432948699326057e+00, 5.5139655214064687e-01, -4.8990051277299940e-01, -3.4926635951494012e+00
+          , -9.6925116789155155e+00, -1.5745152038061331e+01, -9.9352493046059109e-01, 9.4517830645215696e+01
+          }
+    };
+    a_t imresi_03{
+      vc_t{-3.7583808427575212e+02, 9.0801181249849066e+01, 3.6018913915299962e+00, -1.6864299541716473e+01
+          , 1.0089846321224439e+01, -3.6425542565379301e+00, 5.6128150602785309e-01, 4.7851915025080155e-01
+          , 0.0000000000000000e+00, -7.8534132361994757e-01, 1.7555311861860581e+00, -2.2254640329062294e+00
+          , -1.3912074018251386e+00, 1.7899970399493917e+01, -5.9204417397002132e+01, 1.1533747474469016e+02
+          }
+    };
+
+    std::vector<kyosu::complex_t<T>> is_03(N+2);
+    for (int k = 0; k < 16; ++k)
+    {
+      auto c = kyosu::complex(re[k], im[k]);
+      kyosu::cyl_bessel_j(v1, c, is_03);
+      auto vv = -T(N+1.0/3);
+      for (n = 0; n <= N; ++n)
+      {
+         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
+        TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_j(vv, c), 1.e-2);
+      }
+    }
+  }
+};
diff --git a/test/unit/function/cyl_bessel_k.cpp b/test/unit/function/cyl_bessel_k.cpp
index cb9a6ef2..99cfc8de 100644
--- a/test/unit/function/cyl_bessel_k.cpp
+++ b/test/unit/function/cyl_bessel_k.cpp
@@ -9,7 +9,7 @@
 #include <test.hpp>
 
 
-TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real v =  2.3"
               , kyosu::scalar_real_types
               , tts::generate(tts::randoms(-10,10))
               )
@@ -60,19 +60,173 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
           , -3.4235642328623889e-02, 7.3102426796495309e-03, 2.5947374159400196e-03, -2.7525951662690611e-03}
     };
 
-    std::vector<kyosu::complex_t<T>> is_03(N+2);
+    std::vector<kyosu::complex_t<T>> ks_03(N+2);
     for (int k = 0; k < 16; ++k)
     {
       auto c = kyosu::complex(re[k], im[k]);
-      kyosu::cyl_bessel_k(v1, c, is_03);
+      kyosu::cyl_bessel_k(v1, c, ks_03);
       auto vv = T(-0.3);
       for (n = 0; n <= N; ++n)
       {
          auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
-        TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-6) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-6) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
         TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_k(vv, c), 1.e-6);
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_k(-vv, c), 1.e-6);
         vv = eve::dec(vv);
       }
     }
   }
 };
+
+
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+  <typename T>(T )
+{
+  if constexpr(sizeof(T) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+    constexpr int N = 0;
+    T v1 = N+1.0/3;
+    std::cout << "v1 " << v1 << std::endl;
+    auto n = int(v1);
+    using vc_t = std::array<T, 16>;
+    using a_t = std::array<vc_t, N+1>;
+    a_t reresi_03{
+      vc_t{-7.6339268874816966e+01, -9.2192136865580512e+00, 1.6508508118018792e+01, 1.2375096096000567e+01
+          , 2.3503760861380147e+00, -2.5566294609268230e+00, -2.4913144969412047e+00, -6.0742679340229055e-01
+          , 1.5091129245821360e+00, 7.5263258886450732e-02, -1.0402938701852732e-01, -5.4882419994366216e-02
+          , -6.5488350696814700e-03, 7.4931851019796377e-03, 4.9792792996237569e-03, 8.9451308743070568e-04},
+    };
+    a_t imresi_03{
+      vc_t{-2.8143596488672866e+01, 4.2273285530950140e+01, -1.6385214857913763e+01, -2.8174548075889092e+00
+          , 6.6761389661094128e+00, -3.1878291715173703e+00, -3.3868350891838195e-01, 1.6945714309783351e+00
+          , 0.0000000000000000e+00, -3.6796936276059572e-01, 9.2097126162947351e-02, 1.6635565810036414e-02
+          , -2.4007330768135104e-02, 8.2054910688589972e-03, 8.7325515356663085e-04, -2.1515168333318974e-03},
+    };
+
+    std::vector<kyosu::complex_t<T>> ks_03(N+2);
+    for (int k = 0; k < 16; ++k)
+    {
+      auto c = kyosu::complex(re[k], im[k]);
+      kyosu::cyl_bessel_k(v1, c, ks_03);
+      auto vv = T(N+1.0/3);
+      for (n = 0; n <= N; ++n)
+      {
+         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
+        TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-6) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_k(vv, c), 1.e-6);
+
+//         auto z = 2*kyosu::pow(c, T(1.5) )/3;
+//         std::cout << "c " << c << " z " << z << " -> " << kyosu::cyl_bessel_k(vv, z) << std::endl;
+//         vv = eve::dec(vv);
+      }
+    }
+  }
+};
+
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+  <typename T>(T )
+{
+  if constexpr(sizeof(T) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+    constexpr int N = 0;
+    T v1 = N+1.0/3;
+    std::cout << "v1 " << v1 << std::endl;
+    auto n = int(v1);
+    using vc_t = std::array<T, 16>;
+    using a_t = std::array<vc_t, N+1>;
+    a_t reresi_03{
+      vc_t{1.1900413043351881e-03, 1.5495899619039162e-03, 9.2313922815213754e-03, 1.4640077997088552e-02, -1.0136111628580130e-02, -1.0257686840586897e-01, -2.1965377035090466e-01, 1.3644035629075407e-02, 7.5455646229106821e-01, -1.5515848980450035e+00, -4.2927833713114696e+00, -4.0556204704598322e+00, 4.7439243146286634e+00, 2.2123510870256215e+01, 2.8376292912481247e+01, -1.8674126596346056e+01},
+        };
+    a_t imresi_03{
+      vc_t{1.6213718228760029e-03, -4.0459476346553763e-03, 1.9336808997634971e-03, 1.4818266297842454e-02, -4.5874488719344970e-02, 3.8478715861459151e-02, 1.5990522343845420e-01, -7.8491087169627471e-01, -3.2078969068337546e+00, 2.4846307644453542e+00, -1.2426224367177374e-01, -5.8516695697616443e+00, 1.1525345800362158e+01, -4.0835724859787863e+00, -3.0165300008239711e+01, 7.4997855427520179e+01}
+    };
+
+    std::vector<kyosu::complex_t<T>> ks_03(N+2);
+    for (int k = 0; k < 16; ++k)
+    {
+      auto c = -kyosu::complex(re[k], im[k]);
+      kyosu::cyl_bessel_k(v1, c, ks_03);
+      auto vv = T(N+1.0/3);
+      for (n = 0; n <= N; ++n)
+      {
+         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
+        TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-2) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_k(vv, c), 1.e-2);
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_k(-vv,c), 1.e-2);
+      }
+    }
+  }
+};
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+  <typename T>(T )
+{
+  if constexpr(sizeof(T) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+    constexpr int N = 0;
+    T v1 = N-1.0/3;
+    std::cout << "v1 " << v1 << std::endl;
+    auto n = int(v1);
+    using vc_t = std::array<T, 16>;
+    using a_t = std::array<vc_t, N+1>;
+    a_t reresi_03{
+      vc_t{-7.6339268874816966e+01, -9.2192136865580512e+00, 1.6508508118018792e+01, 1.2375096096000567e+01, 2.3503760861380147e+00, -2.5566294609268230e+00, -2.4913144969412047e+00, -6.0742679340229055e-01, 1.5091129245821360e+00, 7.5263258886450732e-02, -1.0402938701852732e-01, -5.4882419994366216e-02, -6.5488350696814700e-03, 7.4931851019796377e-03, 4.9792792996237569e-03, 8.9451308743070568e-04},
+        };
+    a_t imresi_03{
+      vc_t{2.8143596488672866e+01, -4.2273285530950140e+01, 1.6385214857913763e+01, 2.8174548075889092e+00, -6.6761389661094128e+00, 3.1878291715173703e+00, 3.3868350891838195e-01, -1.6945714309783351e+00, 0.0000000000000000e+00, 3.6796936276059572e-01, -9.2097126162947351e-02, -1.6635565810036414e-02, 2.4007330768135104e-02, -8.2054910688589972e-03, -8.7325515356663085e-04, 2.1515168333318974e-03}
+    };
+
+    std::vector<kyosu::complex_t<T>> ks_03(N+2);
+    for (int k = 0; k < 16; ++k)
+    {
+      auto c = kyosu::complex(re[k], -im[k]);
+      kyosu::cyl_bessel_k(v1, c, ks_03);
+      auto vv = T(N+1.0/3);
+      for (n = 0; n <= N; ++n)
+      {
+         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
+        TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-2) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_k(vv, c), 1.e-2);
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_k(-vv, c), 1.e-2);
+      }
+    }
+  }
+};
diff --git a/test/unit/function/cyl_bessel_y.cpp b/test/unit/function/cyl_bessel_y.cpp
new file mode 100644
index 00000000..8ccbbc20
--- /dev/null
+++ b/test/unit/function/cyl_bessel_y.cpp
@@ -0,0 +1,102 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright : KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#include <kyosu/kyosu.hpp>
+#include <test.hpp>
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+  <typename T>(T )
+{
+  if constexpr(sizeof(T) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+    constexpr int N = 0;
+    T v1 = 1.0/3;
+    std::cout << "v1 " << v1 << std::endl;
+    auto n = int(v1);
+    using vc_t = std::array<T, 16>;
+    using a_t = std::array<vc_t, N+1>;
+    a_t reresi_03{
+      vc_t{2.6513657493935034e+02, 1.5122843619345218e+02, 5.0648789106855794e+01, 7.0207071123160771e+00, -4.0157239964866642e+00, -3.8818208974565880e+00, -1.9700774006295092e+00, -8.8015136866517407e-01, -1.2133448538831417e+00, -4.7036021607109069e-02, 1.3293987996432377e+00, 4.1474049533381976e+00, 7.7005870282777540e+00, 4.6857906393290154e+00, -2.8742001471989823e+01, -1.3952369860713620e+02}
+    };
+    a_t imresi_03{
+      vc_t{2.8090371647178273e+02, -1.7536301901048642e+01, -3.3401195686627879e+01, 2.3525800606543957e+01, -9.3284782480155641e+00, 1.9546375587054139e+00, 5.0096629910926505e-01, -9.8992537961477600e-01, 0.0000000000000000e+00, 8.4240922212198033e-01, -1.2574860466731328e+00, 1.8207337449423835e-01, 6.0487612323155489e+00, -2.3373344181702222e+01, 5.1769001839281714e+01, -5.2626238180368290e+01}
+    };
+
+    std::vector<kyosu::complex_t<T>> is_03(N+2);
+    for (int k = 0; k < 16; ++k)
+    {
+      auto c = kyosu::complex(re[k], im[k]);
+      kyosu::cyl_bessel_y(v1, c, is_03);
+      auto vv = T(N+1.0/3);
+      for (n = 0; n <= N; ++n)
+      {
+         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
+        TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_y(vv, c), 1.e-2) << "n " << n  << " k " << k  << " c "<< c <<'\n';
+      }
+    }
+  }
+};
+
+
+TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real"
+              , kyosu::scalar_real_types
+              , tts::generate(tts::randoms(-10,10))
+              )
+  <typename T>(T )
+{
+  if constexpr(sizeof(T) == 8)
+  {
+    std::array<T, 16> re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00
+        , -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01
+        , 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00
+        , 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
+    std::array<T, 16> im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00
+        , 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00
+        , -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00
+        , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00};
+
+    constexpr int N = 0;
+    T v1 = -1.0/3;
+    std::cout << "v1 " << v1 << std::endl;
+    auto n = int(v1);
+    using vc_t = std::array<T, 16>;
+    using a_t = std::array<vc_t, N+1>;
+    a_t reresi_03{
+      vc_t{3.7583805458938821e+02, 9.0801266745274788e+01, -3.6012541742235697e+00, -1.6862151475581715e+01, -1.0085674143682983e+01, -3.6432335856963327e+00, -6.0906440868768552e-01, 2.0873153295804381e-01, -1.0003474349638408e-01, 8.6097480021168360e-01, 1.8102650205095943e+00, 2.2453391061314565e+00, -1.3867486242808820e+00, -1.7899822024160851e+01, -5.9204838602123715e+01, -1.1533771231556815e+02}
+    };
+    a_t imresi_03{
+      vc_t{-8.9163224441390696e+01, 1.2219965550991174e+02, -6.0563732479780832e+01, 1.7841777572774447e+01, -1.1808300271084744e+00, -2.3998139242391443e+00, 1.9741006840075483e+00, -1.1510312786004779e+00, 0.0000000000000000e+00, 3.2456737045081957e-01, 5.2569711539964548e-01, -3.4904700264222477e+00, 9.6878805608040661e+00, -1.5743030177502007e+01, 9.9294471444101262e-01, 9.4517889913758978e+01}
+    };
+
+    std::vector<kyosu::complex_t<T>> is_03(N+2);
+    for (int k = 0; k < 16; ++k)
+    {
+      auto c = kyosu::complex(re[k], im[k]);
+      kyosu::cyl_bessel_y(v1, c, is_03);
+      auto vv = -T(N+1.0/3);
+      for (n = 0; n <= N; ++n)
+      {
+         auto refi_03=  kyosu::complex_t<T>(reresi_03[n][k], imresi_03[n][k]);
+        TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n  << " k " << k  << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n';
+        TTS_RELATIVE_EQUAL(refi_03,  kyosu::cyl_bessel_y(vv, c), 1.e-2);
+      }
+    }
+  }
+};