bessel_jn bug

jfalcou · Dec 8, 2023 · e3f78f1 · e3f78f1
1 parent f5e9bcf
commit e3f78f1
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 8 deletions.
diff --git a/include/kyosu/types/impl/bessel/cb_jyn.hpp b/include/kyosu/types/impl/bessel/cb_jyn.hpp
@@ -477,20 +477,21 @@ namespace kyosu::_
       auto backwardj = [az, nn, n, &cjv](auto z){
         auto j0 = cjv[0];
         auto j1 = cjv[1];;
-        auto m1 = eve::maximum(ini_for_br_1(az, u_t(200)));
-        auto m2 = eve::maximum(ini_for_br_2(n, az, u_t(15)));
-        auto m = eve::if_else( m1 >= n && eve::is_not_nan(m2), m2, m1);
+        auto m1 = ini_for_br_1(az, u_t(200));
+        auto m2 = ini_for_br_2(n, az, u_t(15));
+        auto m0 = eve::if_else( m1 >= n && eve::is_not_nan(m2), m2, m1);
+        auto m = eve::maximum(m0);
         auto cf2 = Z(0);
         auto cf1 = complex(eve::sqrtsmallestposval(eve::as<Z>()));
         Z cf(cf2);
         int k = m;
         auto kgez = eve::is_gez(k);
         while (eve::any(kgez))
         {
-          cf  = kyosu::if_else(kgez,  2*inc(k)*cf1*rec(z)-cf2, cf);
+          cf  = kyosu::if_else(kgez&& k <= m0,  2*inc(k)*cf1*rec(z)-cf2, cf);
           if(k <= n && k > 1) cjv[k] = cf;
-          cf2 = kyosu::if_else(kgez, cf1, cf2);
-          cf1 = kyosu::if_else(kgez, cf, cf1);
+          cf2 = kyosu::if_else(kgez&& k <= m0, cf1, cf2);
+          cf1 = kyosu::if_else(kgez&& k <= m0, cf, cf1);
           k = dec(k);
           kgez = eve::is_gez(k);
         }

diff --git a/include/kyosu/types/impl/besselr/cb_jyr.hpp b/include/kyosu/types/impl/besselr/cb_jyr.hpp
@@ -64,7 +64,7 @@ namespace kyosu::_
         auto cone  = Z{1};
         auto cii   = i(as<Z>());
 
-        auto forward = [n, v0, rz, cjv0, cjv1, &cjv](){
+        auto forward = [n, v0, rz, cjv0 = cjv0, cjv1 = cjv1, &cjv](){
           // az large compared to n : 4*n <  az
           auto cf0 = cjv0;
           auto cf1 = cjv1;
@@ -78,7 +78,7 @@ namespace kyosu::_
           return cf;
         };
 
-        auto backward = [czero, v0, az, n, rz, cjv0, cjv1, &cjv](){
+        auto backward = [czero, v0, az, n, rz, cjv0 = cjv0, cjv1 = cjv1, &cjv](){
           auto m1 = eve::maximum(ini_for_br_1(az, u_t(200)));
           auto m2 = eve::maximum(ini_for_br_2(u_t(n), az, u_t(15)));
           auto m = eve::if_else( m1 >= n && eve::is_not_nan(m2), m2, m1);