Skip to content

Commit

Permalink
airy_ai airy_bi + bug bessel_j for v <0
Browse files Browse the repository at this point in the history
  • Loading branch information
@jtlap
jtlap committed Dec 16, 2023
1 parent f7e1fd0 commit c721ab9
Show file tree
Hide file tree
Showing 10 changed files with 363 additions and 181 deletions.
14 changes: 14 additions & 0 deletions include/kyosu/functions/cyl_bessel_y.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,13 +25,27 @@ 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))
{
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;
Expand Down
49 changes: 9 additions & 40 deletions include/kyosu/types/impl/besselr/airy.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,27 +46,24 @@ namespace kyosu::_
auto third = eve::third(as<u_t>());
auto invpi = eve::inv_pi(as<u_t>());
auto br_im_eq_0 = [](auto z){
// std::cout << " === 3 "<< std::endl;
// return Z(3.550280538878172e-01);
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;
// std::cout << " === 2 "<< z << "arg z " << arg(z) << " arg zet " << arg(zet) << std::endl;
auto argzetgtpio2 = arg(zet) > eve::pio_2(eve::as< u_t>());
auto sqtmz = sqrt(-z/3); //pow(z, u_t(0.5));
auto sqtmz = sqrt(-z/3);
auto r = mulmi(invpi*sqtmz*cyl_bessel_k(third, zet));
// return if_else(imag(z) > 0, r, -r);
r = if_else(izgt0, r, conj(r));
return if_else(argzetgtpio2, -r, r);
// return (cyl_bessel_i(third, zet)+cyl_bessel_i(-third, zet))*sqrt(-z)*third;
};
auto br_re_gt_0 = [third, invpi](auto z){
// std::cout << " === 1 "<< std::endl;
auto zet = 2*pow(z, u_t(1.5) )/3;
return invpi*sqrt(z/3)*cyl_bessel_k(third, zet);
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>());
Expand All @@ -84,14 +81,7 @@ namespace kyosu::_
}
}
}
return r;
// using u_t = eve::underlying_type_t<Z>;
// auto [sqzo3, zeta] = zet(z);
// auto r = eve::inv_pi(as(z))*sqzo3*cyl_bessel_k(eve::third(as<u_t>()), zeta);
// return r;
// return if_else(eve::is_ltz(arg(z)), -r, r);
// return sqzo3*(cyl_bessel_i(eve::third(as<u_t>()), zeta)
// +cyl_bessel_i(-eve::third(as<u_t>()), zeta));
return r;
}

//===-------------------------------------------------------------------------------------------
Expand All @@ -104,28 +94,15 @@ namespace kyosu::_
auto third = eve::third(as<u_t>());
auto invpi = eve::inv_pi(as<u_t>());
auto br_im_eq_0 = [](auto z){
// std::cout << " === 3 "<< std::endl;
// return Z(3.550280538878172e-01);
return Z(eve::airy_bi(real(z)));
};

auto br_re_lt_0 = [third, invpi](auto z){
std::cout << " === 1 "<< std::endl;
auto zet = 2*pow(-z, u_t(1.5) )/3;
std::cout << " === 2 "<< z << "arg z " << arg(z) << " arg zet " << arg(zet) << std::endl;
return (cyl_bessel_j(third, zet)+cyl_bessel_j(-third, zet))*sqrt(-z/3) ;
// 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); //pow(z, u_t(0.5));
// auto r = mulmi(invpi*sqtmz*cyl_bessel_k(third, zet));
// r = if_else(izgt0, r, conj(r));
// return if_else(argzetgtpio2, -r, r);

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){
std::cout << " === 1 "<< std::endl;
auto zet = 2*pow(z, u_t(1.5) )/3;
return sqrt(z/3)*(cyl_bessel_i(third, zet)+cyl_bessel_i(-third, zet));
};
Expand All @@ -145,15 +122,7 @@ namespace kyosu::_
}
}
}
return r;
// using u_t = eve::underlying_type_t<Z>;
// auto [sqzo3, zeta] = zet(z);
// // std::cout << "z" << z << " sqzo3 "<< sqzo3 << " zeta "<< zeta << std::endl;
// // std::cout << "third " << eve::third(as<u_t>()) << std::endl;
// // std::cout << "cyl_bessel_i(eve::third(as<u_t>()), zeta) " << cyl_bessel_i(eve::third(as<u_t>()), zeta) << std::endl;
// // std::cout << "cyl_bessel_i(-eve::third(as<u_t>()), zeta) " << cyl_bessel_i(-eve::third(as<u_t>()), zeta) << std::endl;
// return sqzo3*(cyl_bessel_i(-eve::third(as<u_t>()), zeta)
// -cyl_bessel_i(eve::third(as<u_t>()), zeta));
return r;
}

//===-------------------------------------------------------------------------------------------
Expand Down
43 changes: 26 additions & 17 deletions include/kyosu/types/impl/besselr/cb_jyr.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down Expand Up @@ -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]};
}
}

Expand All @@ -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);
}
Expand All @@ -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);
}
Expand Down Expand Up @@ -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){
Expand Down
13 changes: 9 additions & 4 deletions test/doc/cyl_bessel_j.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
34 changes: 2 additions & 32 deletions test/unit/function/airy_ai.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,45 +41,15 @@ TTS_CASE_WITH ( "Check kyosu::airy_ai over real"
, 0.0000000000000000e+00, -1.5188956587718155e-01, 2.1867029634299019e-02, 6.0030132776007987e-02
, 3.9632838829414818e-02, 2.2033585182583651e-02, 1.3313384345125566e-02, 8.4631437627202440e-03};

auto z = kyosu::complex(0.0);
std::cout << kyosu::airy_ai(z) << std::endl;
for(int i=0; i < 16; ++i)
{
auto c = kyosu::complex(re[i], im[i]);
auto res = kyosu::complex(reres[i], imres[i]);
// std::cout << "i " << i << " c " << c << " arg(c) " << kyosu::arg(c) << " sqrt(c) "<< kyosu::sqrt(c) << std::endl;
TTS_RELATIVE_EQUAL(kyosu::airy_ai(c), res, 1.0e-3) << i << " <- " << c << "arg " << kyosu::arg(c) << '\n';
// auto z = 2*kyosu::pow(c, T(1.5) )/3;
// auto kz = -kyosu::cyl_bessel_k(1/T(3), z)*kyosu::sqrt(c/3)/eve::pi(eve::as(c));
// std::cout << "kz " << kz << std::endl;
// TTS_RELATIVE_EQUAL(kz, res, 1.0e-4) << i << " <- " << 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';

}
}
};

// TTS_CASE_WITH ( "Check kyosu::airy_ai 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_ai;
// auto j0c = airy_ai(c);
// TTS_EQUAL(j0c, kyosu::conj(airy_ai(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_ai(z), kyosu::complex(airy_ai(T(0))));
// }
// };
4 changes: 2 additions & 2 deletions test/unit/function/airy_bi.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real"
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_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';
}
}
};
Expand Down
1 change: 0 additions & 1 deletion test/unit/function/cyl_bessel_i.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]);
Expand Down
Loading

0 comments on commit c721ab9

Please sign in to comment.