besse_y with_alloca

include/kyosu/details/with_alloca.hpp

@@ -0,0 +1,20 @@
+//======================================================================================================================
+/*
+  Kyosu - Complex Without Complexes
+  Copyright : KYOSU Contributors & Maintainers
+  SPDX-License-Identifier: BSL-1.0
+*/
+//======================================================================================================================
+#pragma once
+
+namespace kyosu::_
+{
+
+  template<typename T, std::invocable<T*> F>
+  decltype(auto) with_alloca(auto size, F f)
+  {
+    T* p = (T*)(__builtin_alloca_with_align(size*sizeof(T), 8*alignof(T)));
+    return f(p);
+  }
+
+}

include/kyosu/types/impl/detail/bessel_y.hpp

@@ -7,6 +7,7 @@
 //======================================================================================================================
 #pragma once
 #include <kyosu/types/impl/detail/bessel_j.hpp>
+#include <kyosu/details/with_alloca.hpp>
 
 namespace kyosu::_
 {
@@ -80,15 +81,16 @@ namespace kyosu::_
     auto y = cyl_bessel_y0(z);
     if (nn != 0)
     {
-      auto n = eve::abs(nn);
-      auto js = Js(n, z);
-      using u_t   =  eve::underlying_type_t<Z>;
-      auto twoopi = eve::two_o_pi(eve::as<u_t>());
-      auto b = twoopi*rec(z);
-      for(int i=1; i <= n ; ++i)
-      {
-        y = fms(js[i], y, b)/js[i-1];
-      }
+      int n = eve::abs(nn);
+      auto do_it =  [n, z, &y](auto js){
+        mkjs jj(n, z);
+        for(int i=n; i >= 0 ; --i)  js[i] = jj();
+        using u_t   =  eve::underlying_type_t<Z>;
+        auto twoopi = eve::two_o_pi(eve::as<u_t>());
+        auto b = twoopi*rec(z);
+        for(int i=1; i <= n ; ++i) y = fms(js[i], y, b)/js[i-1];
+      };
+      kyosu::_::with_alloca<Z>(n+1, do_it);
       auto r = if_else(eve::is_gtz(real(z)) && is_real(z) && is_nan(y), complex(real(y)), y);
       r = if_else(is_eqz(z), complex(eve::minf(eve::as<u_t>())), r);
       return nn < 0 ? r*eve::sign_alternate(u_t(n)) : r;