diff --git a/include/kyosu/details/with_alloca.hpp b/include/kyosu/details/with_alloca.hpp index 6eb422a9..1a939e25 100644 --- a/include/kyosu/details/with_alloca.hpp +++ b/include/kyosu/details/with_alloca.hpp @@ -11,19 +11,25 @@ namespace kyosu::_ { - template F> + template> 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); + auto s = std::span(p, size); + return f(s); } - - template - requires(std::invocable) - decltype(auto) with_alloca(auto size, F f) + + template, std::span> F> + decltype(auto) with_alloca(auto size, F f) { T* p1 = (T*)(__builtin_alloca_with_align(size*sizeof(T), 8*alignof(T))); T* p2 = (T*)(__builtin_alloca_with_align(size*sizeof(T), 8*alignof(T))); - return f(p1, p2); + auto s1 = std::span(p1, size); + auto s2 = std::span(p2, size); + return f(s1, s2); } + + + + } diff --git a/include/kyosu/types/impl/bessel/cb_hn.hpp b/include/kyosu/types/impl/bessel/cb_hn.hpp index bb8ee3f0..13e0c136 100644 --- a/include/kyosu/types/impl/bessel/cb_hn.hpp +++ b/include/kyosu/types/impl/bessel/cb_hn.hpp @@ -95,8 +95,11 @@ namespace kyosu::_ template KYOSU_FORCEINLINE auto cb_h12n(N n, Z z, R & h1s, R& h2s) noexcept { + auto an = abs(n); + EVE_ASSERT(std::size(cjv) > size_t(an), "not room enough in js"); + EVE_ASSERT(std::size(cyv) > size_t(an), "not room enough in ys"); cb_jyn(n, z, h1s, h2s); - for(int i=0; i <= eve::abs(n+1); ++i) + for(int i=0; i <= an; ++i) { auto miyst = muli(h2s[i]); h2s[i] = h1s[i]-miyst; diff --git a/include/kyosu/types/impl/bessel/cb_ikn.hpp b/include/kyosu/types/impl/bessel/cb_ikn.hpp index 043a63ad..10245700 100644 --- a/include/kyosu/types/impl/bessel/cb_ikn.hpp +++ b/include/kyosu/types/impl/bessel/cb_ikn.hpp @@ -135,16 +135,18 @@ namespace kyosu::_ //===------------------------------------------------------------------------------------------- // cb_ikn //===------------------------------------------------------------------------------------------- - template KYOSU_FORCEINLINE - auto cb_ikn(N n, Z z, R& is, R& ks) noexcept + template KYOSU_FORCEINLINE + auto cb_ikn(N n, Z z, R1& is, R2& ks) noexcept requires(concepts::complex || eve::floating_ordered_value) { using u_t = eve::underlying_type_t; auto nn = eve::abs(n); + EVE_ASSERT(N(size(is)) > nn, "not room enough in is"); + EVE_ASSERT(N(size(ks)) > nn, "not room enough in ks"); - for(int j=0; j <= nn; ++j) + for(N j=0; j <= nn; ++j) { - is[j] = cb_in(j, z); + is[j] = cb_in(j, z); ks[j] = cb_kn(j, z); } return kumi::tuple{is[n], ks[n]}; @@ -250,10 +252,12 @@ namespace kyosu::_ //===------------------------------------------------------------------------------------------- // cyl_bessel_ikn //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N n, Z z, R& js, R& ys) noexcept + template + auto dispatch(eve::tag_of, N n, Z z, R1& is, R2& ks) noexcept requires(concepts::complex) { - return cb_ikn(n, z, js, ys); + EVE_ASSERT(N(size(is)) > abs(n), "not room enough in is"); + EVE_ASSERT(N(size(ks)) > abs(n), "not room enough in ks"); + return cb_ikn(n, z, is, ks); } } diff --git a/include/kyosu/types/impl/bessel/cb_jyn.hpp b/include/kyosu/types/impl/bessel/cb_jyn.hpp index b5bc72b3..d82fc1a7 100644 --- a/include/kyosu/types/impl/bessel/cb_jyn.hpp +++ b/include/kyosu/types/impl/bessel/cb_jyn.hpp @@ -418,10 +418,10 @@ namespace kyosu::_ auto cb_jyn(N nn, Z z, R1& cjv, R2& cyv) noexcept requires(concepts::complex || eve::floating_ordered_value) { -// EVE_ASSERT(size(cjv) > nn, "not room enough in js"); -// EVE_ASSERT(size(cyv) > nn, "not room enough in ys"); - using u_t = eve::underlying_type_t; auto n = eve::abs(nn); + EVE_ASSERT(N(size(cjv)) > n, "not room enough in cjv"); + EVE_ASSERT(N(size(cyv)) > n, "not room enough in cyv"); + using u_t = eve::underlying_type_t; if (n <= 1) { cjv[0] = cyl_bessel_j0(z); diff --git a/include/kyosu/types/impl/bessel/sb_ikn.hpp b/include/kyosu/types/impl/bessel/sb_ikn.hpp index 83468212..93e436f5 100644 --- a/include/kyosu/types/impl/bessel/sb_ikn.hpp +++ b/include/kyosu/types/impl/bessel/sb_ikn.hpp @@ -93,6 +93,7 @@ namespace kyosu::_ template KYOSU_FORCEINLINE auto sb_i1n(N n, Z z, R & sys) noexcept { + EVE_ASSERT(N(size(sys)) > n, "not enough room in sys"); using e_t = as_real_type_t; return riton(n)*sph_bessel_yn(n,muli(z), sys); } @@ -115,9 +116,10 @@ namespace kyosu::_ return -eve::pio_2(eve::as())*iton(n)*sph_bessel_h1n(n, complex(-imag(z), real(z))); } - template KYOSU_FORCEINLINE - auto sb_kn(int n, Z z, R& sks) noexcept + template KYOSU_FORCEINLINE + auto sb_kn(N n, Z z, R& sks) noexcept { + EVE_ASSERT(N(size(sks)) > n, "not enough room in sks"); using u_t = eve::underlying_type_t; auto iton = [](int n){ if (n%4 == 0) return complex(eve::one(eve::as())); @@ -132,6 +134,8 @@ namespace kyosu::_ template KYOSU_FORCEINLINE auto sb_ikn(int n, Z z, R1& sis, R2& sks) noexcept { + EVE_ASSERT(N(size(sis)) > n, "not enough room in sis"); + EVE_ASSERT(N(size(sks)) > n, "not enough room in sks"); sb_i1n(n, z, sis); sb_kn(n, z, sks); return kumi::tuple{sis[n], sks[n]}; diff --git a/include/kyosu/types/impl/bessel/sb_jyn.hpp b/include/kyosu/types/impl/bessel/sb_jyn.hpp index 48215a13..6f665e7f 100644 --- a/include/kyosu/types/impl/bessel/sb_jyn.hpp +++ b/include/kyosu/types/impl/bessel/sb_jyn.hpp @@ -82,9 +82,10 @@ namespace kyosu::_ return if_else(is_eqz(z), eve::zero, res); } - template KYOSU_FORCEINLINE + template auto sb_jn(N n, Z z, R & sjs) noexcept { + EVE_ASSERT(N(size(sjs)) > n, "not enough room in js"); using u_t = eve::underlying_type_t; auto bd = [z](int n){ auto st = eve::abs(eve::sin(eve::abs(arg(z)))); diff --git a/include/kyosu/types/impl/besselr/cb_h.hpp b/include/kyosu/types/impl/besselr/cb_h.hpp index 874e9ecc..7048adb7 100644 --- a/include/kyosu/types/impl/besselr/cb_h.hpp +++ b/include/kyosu/types/impl/besselr/cb_h.hpp @@ -24,11 +24,14 @@ namespace kyosu::_ //===------------------------------------------------------------------------------------------- // cb_h12 //===------------------------------------------------------------------------------------------- - template KYOSU_FORCEINLINE - auto cb_h12(N n, Z z, R & h1s, R& h2s) noexcept + template KYOSU_FORCEINLINE + auto cb_h12(N n, Z z, R1 & h1s, R2 & h2s) noexcept { + auto an = int(eve::abs(n)); + EVE_ASSERT(int(size(h1s)) > an, "not enough room in h1s"); + EVE_ASSERT(int(size(h2s)) > an, "not enough room in h2s"); cb_jy(n, z, h1s, h2s); - for(int i=0; i <= eve::abs(n+1); ++i) + for(int i=0; i <= an; ++i) { auto miyst = muli(h2s[i]); h2s[i] = h1s[i]-miyst; @@ -72,8 +75,8 @@ namespace kyosu::_ //===------------------------------------------------------------------------------------------- // cyl_bessel_h12 //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N n, Z z, R& h1s, R& h2s) noexcept + template + auto dispatch(eve::tag_of, N n, Z z, R1& h1s, R2& h2s) noexcept requires(concepts::complex) { return cb_h12(n, z, h1s, h2s); diff --git a/include/kyosu/types/impl/besselr/cb_ikr.hpp b/include/kyosu/types/impl/besselr/cb_ikr.hpp index 973edb16..ca0d6581 100644 --- a/include/kyosu/types/impl/besselr/cb_ikr.hpp +++ b/include/kyosu/types/impl/besselr/cb_ikr.hpp @@ -58,6 +58,7 @@ namespace kyosu::_ v = eve::abs(v); //k(-v, z) == K(v, z) DLMF 10.27.1 using u_t = eve::underlying_type_t; auto n = int(v); //n>= 0 + EVE_ASSERT(int(size(ks)) > n, "not enough room in h1s"); auto v0 = v-n; auto vi = v0; for(int j=0; j <= n; ++j) diff --git a/include/kyosu/types/impl/besselr/cb_jyr.hpp b/include/kyosu/types/impl/besselr/cb_jyr.hpp index 7e7f4131..07f2e4bb 100644 --- a/include/kyosu/types/impl/besselr/cb_jyr.hpp +++ b/include/kyosu/types/impl/besselr/cb_jyr.hpp @@ -22,8 +22,8 @@ namespace kyosu::_ { int n = int(v); int an = eve::abs(n); -// EVE_ASSERT(size(cjv) > an, "not enough room in js"); -// EVE_ASSERT(size(cyv) > an, "not enough room in ys"); + EVE_ASSERT(int(size(cjv)) > an, "not enough room in cjv"); + EVE_ASSERT(int(size(cyv)) > an, "not enough room in cyv"); if (eve::is_flint(v)) { return cb_jyn(int(v), z, cjv, cyv); @@ -100,7 +100,7 @@ namespace kyosu::_ } return cjv[n]; }; - + // compute j_{v+i} i = 2...n Z r; auto notdone = eve::true_(as()); @@ -109,7 +109,7 @@ namespace kyosu::_ { last_interval(backward, notdone, r); } - + // compute y_{v+i} i = 2...n auto twoopiz = twoopi*rz; auto fouropisqz = 2*twoopiz*rz; @@ -142,8 +142,8 @@ namespace kyosu::_ //===------------------------------------------------------------------------------------------- // cyl_bessel_jy //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N v, Z z, R& js, R& ys) noexcept + template + auto dispatch(eve::tag_of, N v, Z z, R1& js, R2& ys) noexcept { return cb_jy(v, z, js, ys); } @@ -200,7 +200,7 @@ namespace kyosu::_ //===------------------------------------------------------------------------------------------- // cyl_bessel_y //===------------------------------------------------------------------------------------------- - template + template auto dispatch(eve::tag_of, N v, Z z) noexcept { auto n = int(abs(v))+1; diff --git a/include/kyosu/types/impl/detail/bessel_jr.hpp b/include/kyosu/types/impl/detail/bessel_jr.hpp deleted file mode 100644 index a5a62bf5..00000000 --- a/include/kyosu/types/impl/detail/bessel_jr.hpp +++ /dev/null @@ -1,60 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include -#include -#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_bessel_jnu - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N nu, Z z) noexcept - { - if constexpr(concepts::complex ) - { - auto doit = [nu, z](auto js, auto ys){ - auto [j, y] = cb_jyr(eve::abs(nu), z, js, ys); - return kumi::tuple{j, y}; - }; - auto [j, y] = with_alloca(int(nu)+1, doit); - - if(eve::is_ltz(nu)) - { - if (eve::is_flint(nu)) - { - return cyl_bessel_jn(int(nu), z); - } - else - { - auto [s, c] = sinpicospi(z); - return j*s-y*c; - } - } - else - { - if (eve::is_flint(nu)) - { - return cyl_bessel_jn(int(nu), z); - } - else - { - return j; - } - } - } - else - { - return cayley_extend_rev(cyl_bessel_jnu, nu, z); - } - } - -} diff --git a/include/kyosu/types/impl/detail/bessel_jyr.hpp b/include/kyosu/types/impl/detail/bessel_jyr.hpp deleted file mode 100644 index 6bd9d7bd..00000000 --- a/include/kyosu/types/impl/detail/bessel_jyr.hpp +++ /dev/null @@ -1,132 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include -#include -#include -#include -#include - -namespace kyosu::_ -{ - //===------------------------------------------------------------------------------------------- - // cb_jyr - //===------------------------------------------------------------------------------------------- - // here nu is always > 0 and z is at most complex - template - auto cb_jyr( N v, Z z, R& cjv, R& cyv) noexcept - { - EVE_ASSERT(is_not_flint(v), "v must not be a flint"); - int n = int(v); -// EVE_ASSERT(size(cjv) > n, "not enough room in js"); -// EVE_ASSERT(size(cyv) > n, "not enough room in ys"); - auto v0 = frac(v); // v = n+v0 0 < v0 < 1 as v is not a flint - // compute cjv0, cjv1, cyv0, cyv1; - auto [cjv0, cyv0, cjv1, cyv1] = kyosu::_::cb_jyr01(v0, z); - cjv[0] = cjv0; - cyv[0] = cyv0; - if (n > 0) - { - cjv[1] = cjv1; - cyv[1] = cyv1; - } - if(n == 0) { - return kumi::tuple{cjv0, cyv0}; - } - else if (n == 1) - { - return kumi::tuple{cjv1, cyv1}; - } - else if (n >= 2) - { - // now jv0 jv1 yv0 and yv1 had been computed - // we make forward or backward recurrence to compute the others - using u_t = eve::underlying_type_t; - auto eps = 16*eve::eps(as()); - auto twoopi = eve::two_o_pi(as()); - auto az = abs(z); - auto z1 = z; - auto z2 = sqr(z); - auto v1 = inc(v0); // 1 < v1 < 2 - auto rz = rec(z); - - auto czero = Z{0}; - auto cone = Z{1}; - auto cii = i(as()); - - auto forward = [n, v0, rz, cjv0, cjv1, &cjv](){ - // az large compared to n : 4*n < az - auto cf0 = cjv0; - auto cf1 = cjv1; - auto cf = cf1; - for (int k=2; k <= n; ++k) { - cf = 2*dec(k+v0)*cf1*rz-cf0; - cjv[k] = cf; - cf0 = cf1; - cf1 = cf; - } - return cf; - }; - - auto backward = [czero, v0, az, n, rz, cjv0, 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); - auto cf2 = czero; - auto cf1 = eve::sqrtsmallestposval(eve::as()); - auto cf = cf1; - auto bn(cf2); - - for (int k=m; k>=0; --k) - { - cf = u_t(2)*inc(v0+u_t(k))*cf1*rz-cf2; - if (k <= n) cjv[k] = cf; - cf2 = cf1; - cf1 = cf; - } - auto cs = if_else(abs(cjv0) > abs(cjv1), cjv0/cf, cjv1/cf2); - for (int k=0; k<=n; ++k) { - cjv[k] *= cs; ; - } - return cjv[n]; - }; - // compute j_{v+i} i = 2...n - Z r; - auto notdone = eve::true_(as()); - notdone = next_interval(forward, notdone, 4*n <= az, r); - if( eve::any(notdone) ) - { - last_interval(backward, notdone, r); - } - // compute y_{v+i} i = 2...n - auto twoopiz = twoopi*rz; - auto fouropisqz = 2*twoopiz*rz; - for (int k=2; k<=n; ++k) - { - cyv[k] = if_else(kyosu::abs(cjv[k-1]) > kyosu::abs(cjv[k-2]) - , (cjv[k]*cyv[k-1]-twoopiz)/cjv[k-1] - , (cjv[k]*cyv[k-2]-dec(v0+k)*fouropisqz)/cjv[k-2] - ); - } - return kumi::tuple{cjv[n], cyv[n]}; - } - } - - //===------------------------------------------------------------------------------------------- - // cyl_bessel_jy - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N v, Z z, R& js, R& ys) noexcept - { - if (is_flint(v)) - return cb_jyn(int(v), z, js, ys); - else - return cb_jyr(v, z, js, ys); - } - -} diff --git a/include/kyosu/types/impl/detail/bessel_jyr01.hpp b/include/kyosu/types/impl/detail/bessel_jyr01.hpp deleted file mode 100644 index 5701326f..00000000 --- a/include/kyosu/types/impl/detail/bessel_jyr01.hpp +++ /dev/null @@ -1,154 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include -#include -#include - -namespace kyosu::_ -{ - //===------------------------------------------------------------------------------------------- - // cb_jyr01 - // this internal routine computes jv0 jv1 yv0 and yv1 with v0 in ]0, 1[ and v1 = v0+1 - //===------------------------------------------------------------------------------------------- - // here nu is always > 0 and z is at most complex - template - auto cb_jyr01( N v0, Z z) noexcept - { - EVE_ASSERT(eve::is_gtz(v0) && eve::is_ltz(oneminus(v0)), "v0 is not in ]0, 1["); - auto v1 = inc(v0); // 1 < v1 < 2 - - using u_t = eve::underlying_type_t; - const auto eps = 16*eve::eps(as()); - const auto hlf = eve::half(as()); - const auto two = u_t(2); - const auto three = u_t(3); - const auto four = u_t(4); - const auto quarter = u_t(0.25); - const auto twoopi = eve::two_o_pi(as()); - const auto invpi = eve::inv_pi(as()); - - auto az = abs(z); - auto isltzrz = eve::is_ltz(real(z)); - auto isltziz = eve::is_ltz(imag(z)); - auto z1 = if_else(isltzrz, -z, z); - auto z2 = sqr(z); - - const Z cone = eve::one(as()); - const Z cnan = eve::nan(as()); - - Z cjv0(cnan), cjv1(cnan), cyv0(cnan), cyv1(cnan); - - - auto br_lt12 = [&](){ - auto cjv01 = [quarter, hlf, cone, eps, z1, z2](auto vl){ - auto cjvl = cone; - auto cr = cone; - for (int k=1;k<=40;++k) - { - cr *= -quarter*z2/(k*(k+vl)); - cjvl += cr; - if (eve::all(kyosu::abs(cr) < kyosu::abs(cjvl)*eps)) break; - } - auto vg = inc(vl); - auto ga = eve::tgamma(vg); - auto ca = pow(hlf*z1,vl)/ga; - return ca*cjvl; - }; - cjv0 = cjv01(v0); - cjv1 = cjv01(v1); - - auto cju01 = [&](auto vl){ - auto cjvl = cone; - auto cr = cone; - for (int k=1;k<=40;++k) { - cr *= -quarter*z2/(k*(k-vl)); - cjvl += cr; - if (eve::all(kyosu::abs(cr) < kyosu::abs(cjvl)*eps)) break; - } - auto vg = eve::oneminus(vl); - auto gb = eve::tgamma(vg); - auto cb = pow(two/z1,vl)/gb; - return cjvl*cb; - }; - auto cju0 = cju01(v0); - auto cju1 = cju01(v1); - - auto [s0, c0] = sinpicospi(v0); - auto [s1, c1] = sinpicospi(v1); - cyv0 = fms(cjv0, c0, cju0)/s0; - cyv1 = fms(cjv1, c1, cju1)/s1; - return kumi::tuple{cjv0, cyv0, cjv1, cyv1}; - }; - - auto br_gt12 = [&](){ - auto bds = eve::if_else(az >= 50.0, 8, eve::if_else(az >= 35.0, 10, 11)); - auto cjyv01 = [&](auto v0p){ - auto vv = sqr(two*v0p); - auto cpz = cone; - auto crp = cone; - size_t k = 1; - auto bound_not_reached = u_t(k) <= bds; - while (eve::any(bound_not_reached)) - { - crp = u_t(-0.78125e-2)*crp*(vv-sqr(four*k-three))* - (vv-sqr(dec(four*k)))/(k*dec(two*k)*z2); - cpz = if_else(bound_not_reached, cpz+ crp, cpz); - bound_not_reached = u_t(++k) <= bds; - } - auto cqz = cone; - auto crq = cone; - k = 1; - bound_not_reached = u_t(k) <= bds; - while (eve::any(bound_not_reached)) - { - crq = u_t(-0.78125e-2)*crq*(vv-sqr(dec(four*k)))* - (vv-sqr(inc(four*k)))/(k*inc(two*k)*z2); - cqz += crq; - cqz = if_else(bound_not_reached, cqz+crp, cqz); - bound_not_reached = u_t(++k) <= bds; - } - cqz *= u_t(0.125)*dec(vv)/z1; - auto z1opi = z1*invpi; - auto zk = z1opi-(hlf*(v0p+hlf)); - auto ca0 = sqrt(twoopi/z1); - auto [ csk, cck] = sinpicospi(zk); - return kumi::tuple{ca0*(cpz*cck-cqz*csk) - , ca0*(cpz*csk+cqz*cck)}; - }; - kumi::tie(cjv0, cyv0) = cjyv01(v0); - kumi::tie(cjv1, cyv1) = cjyv01(v1); - return kumi::tuple{cjv0, cyv0, cjv1, cyv1}; - }; - - auto r = kumi::tuple{cjv0, cyv0, cyv1, cjv1}; - auto notdone = eve::true_(as()); - if( eve::any(notdone) ) - { - notdone = next_interval(br_lt12, notdone, az <= u_t(12), r); - if( eve::any(notdone) ) - { - last_interval(br_gt12, notdone, r); - } - } - kumi::tie(cjv0, cyv0, cjv1, cyv1) = r; - - if( eve::any(isltzrz) ) //treating real(z) < 0 - { - auto eipiv0 = exp_ipi(v0); - auto eipimv0 = rec(eipiv0); - auto eipiv1 = -eipiv0; - auto eipimv1 = -eipimv0; - cyv0 = if_else(isltzrz, cyv0*eipimv0+2*muli(real(eipimv0))*cjv0, cyv0); - cyv1 = if_else(isltzrz, cyv1*eipimv1+2*muli(real(eipimv1))*cjv1, cyv1); - cjv0 = if_else(isltzrz, cjv0*eipiv0, cjv0); - cjv1 = if_else(isltzrz, cjv1*eipiv1, cjv1); - } - return kumi::tuple{cjv0, cyv0, cjv1, cyv1}; - } -} diff --git a/include/kyosu/types/impl/detail/bessel_jyr_old.hpp b/include/kyosu/types/impl/detail/bessel_jyr_old.hpp deleted file mode 100644 index 507c87e3..00000000 --- a/include/kyosu/types/impl/detail/bessel_jyr_old.hpp +++ /dev/null @@ -1,248 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include - -//#include -#include -#include - -namespace kyosu::_ -{ - //===------------------------------------------------------------------------------------------- - // cyl_bessel_jyr - //===------------------------------------------------------------------------------------------- - // here nu is always > 0 and z is at most complex - template - auto cb_jy( N v, Z z, R& cjv, R& cyv) noexcept - { - EVE_ASSERT(is_not_flint(v), "v must not be a flint"); - int n = int(v); -// EVE_ASSERT(size(js) > n, "not enough room in js"); -// EVE_ASSERT(size(ys) > n, "not enough room in ys"); - using u_t = eve::underlying_type_t; - auto el = u_t(0.5772156649015329); - auto rp2 = u_t( 0.63661977236758); - auto eps = 16*eve::eps(as()); - auto hlf = eve::half(as()); - auto two = u_t(2); - auto three = u_t(3); - auto four = u_t(4); - auto quarter = u_t(0.25); - auto pi = eve::pi(as()); - auto twopi = eve::two_pi(as()); - auto twoopi = eve::two_o_pi(as()); - auto invpi = eve::inv_pi(as()); - auto cspv = cbrt(eve::smallestposval(as())); - auto az = abs(z); - auto z1 = z; - auto z2 = sqr(z); - auto v0 = frac(v); // v = n+v0 0 < v0 < 1 as v is not a flint - auto v1 = inc(v0); // 1 < v1 < 2 - - auto czero = Z{0}; - auto cone = Z{1}; - auto cii = i(as()); - - auto smallz = [z, v0, v1, n, invpi, hlf, &cjv, &cyv](){ - // DLMF 10.7.3 and 10.7.4 - auto lzh = kyosu::log(z*hlf); - auto vk = v0; - auto lgvk = log_gamma(v0); - for(int k=0; k <= n; ++k) - { - auto lzhtovk = lzh*vk; - auto lgvkp1 = lgvk+log(vk); - cjv[k] = kyosu::exp(lzhtovk-lgvkp1); - cyv[k] = -invpi*kyosu::exp(lgvk-lzhtovk); - lgvk = lgvkp1; - vk = inc(vk); - } - return kumi::tuple{cjv[n], cyv[n]}; - }; - - if constexpr(eve::scalar_value) - { - if (az < cspv) - { - return smallz(); - } - -// if(eve::is_ltz(real(z1))) z1 = -z; - - Z cjv0, cjv1, cyv0, cyv1; - if (az <= 12) - { - auto cjv01 = [quarter, hlf, cone, eps, z1, z2](auto vl){ - auto cjvl = cone; - auto cr = cone; - for (int k=1;k<=40;++k) - { - cr *= -quarter*z2/(k*(k+vl)); - cjvl += cr; - if (abs(cr) < abs(cjvl)*eps) break; - } - auto vg = inc(vl); - auto ga = eve::tgamma(vg); - auto ca = pow(hlf*z1,vl)/ga; - return ca*cjvl; - }; - cjv0 = cjv01(v0); - cjv1 = cjv01(v1); - - auto cju01 = [two, cone, quarter, eps, z1, z2](auto vl){ - auto cjvl = cone; - auto cr = cone; - for (int k=1;k<=40;++k) { - cr *= -quarter*z2/(k*(k-vl)); - cjvl += cr; - if (kyosu::abs(cr) < kyosu::abs(cjvl)*eps) break; - } - auto vg = eve::oneminus(vl); - auto gb = eve::tgamma(vg); - auto cb = pow(two/z1,vl)/gb; - return cjvl*cb; - }; - auto cju0 = cju01(v0); - auto cju1 = cju01(v1); - - auto [s0, c0] = sinpicospi(v0); - auto [s1, c1] = sinpicospi(v1); - cyv0 = fms(cjv0, c0, cju0)/s0; - cyv1 = fms(cjv1, c1, cju1)/s1; - } - else - { - auto kz = 11; - if (az >= 35.0) kz = 10; - else if (az >= 50.0) kz = 8; - - auto cjyv01 = [twoopi, pi, quarter, hlf, three, two, cone, four, invpi, kz, z1, z2](auto v0p){ - auto vv = sqr(two*v0p); - auto cpz = cone; - auto crp = cone; - for (int k=1;k<=kz;++k) - { - crp = u_t(-0.78125e-2)*crp*(vv-sqr(four*k-three))* - (vv-sqr(dec(four*k)))/(k*dec(two*k)*z2); - cpz += crp; - } - auto cqz = cone; - auto crq = cone; - for (int k=1;k<=kz;++k) - { - crq = u_t(-0.78125e-2)*crq*(vv-sqr(dec(four*k)))* - (vv-sqr(inc(four*k)))/(k*inc(two*k)*z2); - cqz += crq; - } - cqz *= u_t(0.125)*dec(vv)/z1; - auto z1opi = z1*invpi; - auto zk = z1opi-(hlf*(v0p+hlf)); - auto ca0 = sqrt(twoopi/z1); - auto [ csk, cck] = sinpicospi(zk); - return kumi::tuple{ca0*(cpz*cck-cqz*csk) - , ca0*(cpz*csk+cqz*cck)}; - }; - kumi::tie(cjv0, cyv0) = cjyv01(v0); - kumi::tie(cjv1, cyv1) = cjyv01(v1); - } - -// if (real(z) < 0.0) { -// std::cout << "icitte" << std::endl; -// auto cfac0 = exp_ipi(v0); -// auto cfac1 = exp_ipi(v1); -// auto rcfac0= rec(cfac0); -// auto rcfac1= rec(cfac1); -// auto icospv0 = 2*muli(cospi(v0)); -// auto icospv1 = 2*muli(cospi(v1)); -// if (imag(z) < 0.0) { -// cyv0 = cfac0*cyv0-icospv0*cjv0; -// cyv1 = cfac1*cyv1-icospv1*cjv1; -// cjv0 *= rcfac0; -// cjv1 *= rcfac1; -// } -// else if (imag(z) > 0.0) { -// cyv0 = cyv0*rcfac0+icospv0*cjv0; -// cyv1 = cyv1*rcfac1+icospv1*cjv1; -// cjv0 *= cfac0; -// cjv1 *= cfac1; -// } -// } - -// std::cout<< std::endl << "besselj(" << v0 << ") = " << cjv0 << std::endl; -// std::cout<< std::endl << "besselj(" << v0+1 << ") = " << cjv1 << std::endl; - - // now jv0 jv1 yv0 and yv1 had been computed - // we make forward or backward recurrence - cjv[0] = cjv0; - cyv[0] = cyv0; - cjv[1] = cjv1; - cyv[1] = cyv1; - - if(n == 0) return kumi::tuple{cjv0, cyv0}; - else if (n == 1) return kumi::tuple{cjv1, cyv1}; - else if (n >= 2) - { - auto rz = rec(z); - - - auto forward = [n, v0, rz, cjv0, cjv1, &cjv](){ - // az large compared to n : 4*n < az - auto cf0 = cjv0; - auto cf1 = cjv1; - auto cf = cf1; - for (int k=2; k<= n; ++k) { - cf = 2*dec(k+v0)*cf1*rz-cf0; - if (k <= n) cjv[k] = cf; - cf0 = cf1; - cf1 = cf; - } - return cf; - }; - - auto backward = [v0, az, n, rz, cjv0, cjv1, &cjv](){ - auto m1 = ini_for_br_1(az, u_t(200)); - auto m2 = 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); - auto cf2 = Z(0); - auto cf1 = complex(eve::sqrtsmallestposval(eve::as< u_t>())); - auto cf = cf1; - auto bn(cf2); - - for (int k=m; k>=0; --k) - { - cf = 2*inc(v0+k)*cf1*rz-cf2; - if (k <= n) cjv[k] = cf; - cf2 = cf1; - cf1 = cf; - } - auto cs = (abs(cjv0) > abs(cjv1)) ? cjv0/cf : cjv1/cf2; - for (int k=0; k<=n; ++k) { - cjv[k] *= cs; ; - } - return cjv[n]; - }; - auto res_cjv = (n <= (int)(quarter*az)) ? forward() : backward(); - Z res_cyv; - - auto ya0 = kyosu::abs(cyv0); - int lb = 0; - auto cg0 = cyv0; - auto cg1 = cyv1; - for (int k=2; k<=n; ++k) - { - if(kyosu::abs(cjv[k-1]) > kyosu::abs(cjv[k-2])) - cyv[k] = (cjv[k]*cyv[k-1]-2/(pi*z))/cjv[k-1]; - else - cyv[k] = (cjv[k]*cyv[k-2]-4*dec(v0+k)/(pi*sqr(z)))/cjv[k-2]; - } - return kumi::tuple{res_cjv, cyv[n]}; - } - } - } -} diff --git a/include/kyosu/types/impl/detail/bessel_yr.hpp b/include/kyosu/types/impl/detail/bessel_yr.hpp deleted file mode 100644 index 38abe6a3..00000000 --- a/include/kyosu/types/impl/detail/bessel_yr.hpp +++ /dev/null @@ -1,58 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include -#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_bessel_ynu - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N nu, Z z) noexcept - { - if constexpr(concepts::complex ) - { - auto doit = [nu, z](auto js, auto ys){ - auto [j, y] = cb_jyr(eve::abs(nu), z, js, ys); - return kumi::tuple{j, y}; - }; - auto [j, y] = with_alloca(int(nu)+1, doit); - - if(eve::is_ltz(nu)) - { - if (eve::is_flint(nu)) - { - return cyl_bessel_yn(int(nu), z); - } - else - { - auto [s, c] = sinpicospi(z); - return j*s+y*c; - } - } - else - { - if (eve::is_flint(nu)) - { - return cyl_bessel_yn(int(nu), z); - } - else - { - return y; - } - } - } - else - { - return cayley_extend_rev(cyl_bessel_ynu, nu, z); - } - } -}