diff --git a/doc/index.md b/doc/index.md index 79a448e4..87031da1 100644 --- a/doc/index.md +++ b/doc/index.md @@ -129,7 +129,8 @@ complex functions and rotation related quaternion usage. * callables usable with all cayley_dickson types - Most **EVE** arithmetic and math functions are provided. + Most **EVE** arithmetic and math functions are provided. In particular, real analytic functions of one variable + can be quite naturally extended using the polar form described above from a mere complex implementation. | | | | | | |--------------|---------------|-----------------|------------|-----------------| @@ -137,38 +138,31 @@ complex functions and rotation related quaternion usage. | [acotpi](@ref kyosu::acotpi ) | [atanpi](@ref kyosu::atanpi ) | [acoth](@ref kyosu::acoth ) | [arg](@ref kyosu::arg ) | [acsc](@ref kyosu::acsc ) | | [acscpi](@ref kyosu::acscpi ) | [acsch](@ref kyosu::acsch ) | [asec](@ref kyosu::asec ) | [asecpi](@ref kyosu::asecpi ) | [asech](@ref kyosu::asech ) | | [asin](@ref kyosu::asin ) | [asinpi](@ref kyosu::asinpi ) | [asinh](@ref kyosu::asinh ) | [atan](@ref kyosu::atan ) | [atanh](@ref kyosu::atanh ) | - | [average](@ref kyosu::average ) | [ceil](@ref kyosu::ceil ) | [conj](@ref kyosu::conj ) | [cos](@ref kyosu::cos ) | [cosh](@ref kyosu::cosh ) | - | [cospi](@ref kyosu::cospi ) | [cot](@ref kyosu::cot ) | [cotpi](@ref kyosu::cotpi ) | [coth](@ref kyosu::coth ) | [convert](@ref kyosu::convert ) | - | [csc](@ref kyosu::csc ) | [cscpi](@ref kyosu::cscpi ) | [csch](@ref kyosu::csch ) | [dec](@ref kyosu::dec ) | [dist](@ref kyosu::dist ) | - | [dot](@ref kyosu::dot ) | [exp](@ref kyosu::exp ) | [exp10](@ref kyosu::exp10 ) | [exp2](@ref kyosu::exp2 ) | [exp_i](@ref kyosu::exp_i ) | - | [exp_ipi](@ref kyosu::exp_ipi ) | [expm1](@ref kyosu::expm1 ) | [expmx2](@ref kyosu::expmx2 ) | [expx2](@ref kyosu::expx2 ) | [fam](@ref kyosu::fam ) | - | [floor](@ref kyosu::floor ) | [fma](@ref kyosu::fma ) | [fms](@ref kyosu::fms ) | [fnma](@ref kyosu::fnma ) | [fnms](@ref kyosu::fnms ) | - | [frac](@ref kyosu::frac ) | [fsm](@ref kyosu::fsm ) | [from_polar](@ref kyosu::from_polar ) | [horner](@ref kyosu::horner ) | [hypot](@ref kyosu::hypot ) | - | [if_else](@ref kyosu::if_else ) | [inc](@ref kyosu::inc ) |||| + | [average](@ref kyosu::average ) | [beta](@ref kyosu::beta ) | [ceil](@ref kyosu::ceil ) | [conj](@ref kyosu::conj ) | [cos](@ref kyosu::cos ) | + | [cosh](@ref kyosu::cosh ) | [cospi](@ref kyosu::cospi ) | [cot](@ref kyosu::cot ) | [cotpi](@ref kyosu::cotpi ) | [coth](@ref kyosu::coth ) | + | [convert](@ref kyosu::convert ) | [csc](@ref kyosu::csc ) | [cscpi](@ref kyosu::cscpi ) | [csch](@ref kyosu::csch ) | [dec](@ref kyosu::dec ) | + | [deta](@ref kyosu::deta ) | [digamma](@ref kyosu::digamma ) | [dist](@ref kyosu::dist ) | [dot](@ref kyosu::dot ) | [erf](@ref kyosu::erf ) | + | [erfcx](@ref kyosu::erfcx ) | [erfi](@ref kyosu::erfi ) | [eta](@ref kyosu::eta ) | [exp](@ref kyosu::exp ) | [exp10](@ref kyosu::exp10 ) | + | [exp2](@ref kyosu::exp2 ) | [exp_i](@ref kyosu::exp_i ) | [exp_ipi](@ref kyosu::exp_ipi ) | [expm1](@ref kyosu::expm1 ) | [expmx2](@ref kyosu::expmx2 ) | + | [expx2](@ref kyosu::expx2 ) | [faddeeva](@ref kyosu::faddeeva ) | [fam](@ref kyosu::fam ) | [floor](@ref kyosu::floor ) | [fma](@ref kyosu::fma ) | + | [fms](@ref kyosu::fms ) | [fnma](@ref kyosu::fnma ) | [fnms](@ref kyosu::fnms ) | [frac](@ref kyosu::frac ) | [fsm](@ref kyosu::fsm ) | + | [from_polar](@ref kyosu::from_polar ) | [horner](@ref kyosu::horner ) | [hypot](@ref kyosu::hypot ) | [if_else](@ref kyosu::if_else )| [inc](@ref kyosu::inc ) | | [ipart](@ref kyosu::ipart ) | [is_denormal](@ref kyosu::is_denormal ) | [is_equal](@ref kyosu::is_equal ) | [is_eqz](@ref kyosu::is_eqz ) | [is_finite](@ref kyosu::is_finite ) | | [is_infinite](@ref kyosu::is_infinite ) | [is_imag](@ref kyosu::is_imag ) | [is_nan](@ref kyosu::is_nan ) | [is_nez](@ref kyosu::is_nez ) | [is_not_denormal](@ref kyosu::is_not_denormal ) | - | [is_not_equal](@ref kyosu::is_not_equal ) | [is_not_finite](@ref kyosu::is_not_finite ) | [is_not_infinite](@ref kyosu::is_not_finite ) | [is_not_nan](@ref kyosu::is_not_nan ) | [is_not_real](@ref kyosu::is_not_real ) | + | [is_not_equal](@ref kyosu::is_not_equal ) | [is_not_finite](@ref kyosu::is_not_finite ) | [is_not_infinite](@ref kyosu::is_not_finite ) | [is_not_nan](@ref kyosu::is_not_nan ) | [is_not_real](@ref kyosu::is_not_real ) | | [is_pure](@ref kyosu::is_pure ) | [is_real](@ref kyosu::is_real ) | [is_unitary](@ref kyosu::is_unitary ) | [jpart](@ref kyosu::jpart ) | [kpart](@ref kyosu::kpart ) | - | [ldiv](@ref kyosu::ldiv ) | [lerp](@ref kyosu::lerp ) | [log](@ref kyosu::log ) | [log10](@ref kyosu::log10 ) | [log1p](@ref kyosu::log1p ) | - | [log_abs](@ref kyosu::log_abs ) | [log2](@ref kyosu::log2 ) | [lpnorm](@ref kyosu::lpnorm ) | [manhattan](@ref kyosu::manhattan ) | [minus](@ref kyosu::minus ) | - | [lpart](@ref kyosu::lpart ) | [lipart](@ref kyosu::lipart ) | [ljpart](@ref kyosu::ljpart ) | [lkpart](@ref kyosu::lkpart ) | | - | [nearest](@ref kyosu::nearest ) | [oneminus](@ref kyosu::oneminus ) | [pow](@ref kyosu::pow ) | [pow1p](@ref kyosu::pow1p ) | [pow_abs](@ref kyosu::pow_abs ) | - | [powm1](@ref kyosu::powm1 ) | [proj](@ref kyosu::proj ) | [pure](@ref kyosu::imag ) | [radinpi](@ref kyosu::radinpi )| [real](@ref kyosu::real ) | - | [rec](@ref kyosu::rec ) | [reldist](@ref kyosu::reldist ) | [reverse_horner](@ref kyosu::reverse_horner ) | [right_horner](@ref kyosu::right_horner ) | [right_reverse_horner](@ref kyosu::right_reverse_horner ) | - | [sec](@ref kyosu::sec ) | [secpi](@ref kyosu::secpi ) | [sech](@ref kyosu::sech ) | | | - | [sign](@ref kyosu::sign ) | [sin](@ref kyosu::sin ) | [sinc](@ref kyosu::sinc ) | [sincos](@ref kyosu::sincos ) | [sinpi](@ref kyosu::sinpi ) | - | [sinpicospi](@ref kyosu::sinpicospi ) | [sinh](@ref kyosu::sinh ) | [sinhcosh](@ref kyosu::sinhcosh ) | [slerp](@ref kyosu::slerp ) | [sqr](@ref kyosu::sqr ) | - | [sqr_abs](@ref kyosu::sqr_abs ) | [sqrt](@ref kyosu::sqrt ) | [tan](@ref kyosu::tan ) | [tanpi](@ref kyosu::tanpi ) | [tanh](@ref kyosu::tanh ) | - | [to_polar](@ref kyosu::to_polar ) | [trunc](@ref kyosu::trunc ) | | | | - - * callables usable with complex or real only. These are mainly implementation of some classical meromorphic functions. - - | | | | | | - |------------------|------------------|--------------------|------------------|---------------| - | [beta](@ref kyosu::beta ) | [deta](@ref kyosu::deta ) | [digamma](@ref kyosu::digamma ) | [erf](@ref kyosu::erf ) | [erfcx](@ref kyosu::erfcx ) | - | [erfi](@ref kyosu::erfi ) | [eta](@ref kyosu::eta ) | [faddeeva](@ref kyosu::faddeeva ) | [lambda](@ref kyosu::lambda ) | [lbeta](@ref kyosu::lbeta ) | - | [log_abs_gamma](@ref kyosu::log_abs_gamma ) | [log_gamma](@ref kyosu::log_gamma ) | [lrising_factorial](@ref kyosu::lrising_factorial ) | [rising_factorial](@ref kyosu::rising_factorial ) | [tgamma](@ref kyosu::tgamma ) | - | [zeta](@ref kyosu::zeta ) | | | | | + | [ldiv](@ref kyosu::ldiv ) | [lambda](@ref kyosu::lambda ) | [lbeta](@ref kyosu::lbeta ) | [lerp](@ref kyosu::lerp ) | [lipart](@ref kyosu::lipart ) | + | [ljpart](@ref kyosu::ljpart ) | [lkpart](@ref kyosu::lkpart ) | [log](@ref kyosu::log ) | [log2](@ref kyosu::log2 ) | [lpart](@ref kyosu::lpart ) | + | [lipart](@ref kyosu::lipart ) | [ljpart](@ref kyosu::ljpart ) | [lkpart](@ref kyosu::lkpart ) | [lpnorm](@ref kyosu::lpnorm ) | [lrising_factorial](@ref kyosu::lrising_factorial ) | + | [manhattan](@ref kyosu::manhattan ) | [minus](@ref kyosu::minus ) | [nearest](@ref kyosu::nearest ) | [oneminus](@ref kyosu::oneminus )| [pow](@ref kyosu::pow ) | + | [pow1p](@ref kyosu::pow1p ) | [pow_abs](@ref kyosu::pow_abs ) | [powm1](@ref kyosu::powm1 ) | [proj](@ref kyosu::proj ) | [pure](@ref kyosu::imag ) | + | [radinpi](@ref kyosu::radinpi ) | [real](@ref kyosu::real ) | [rec](@ref kyosu::rec ) | [reldist](@ref kyosu::reldist )| [reverse_horner](@ref kyosu::reverse_horner ) | + | [right_horner](@ref kyosu::right_horner ) | [right_reverse_horner](@ref kyosu::right_reverse_horner ) | [rising_factorial](@ref kyosu::rising_factorial ) | [sec](@ref kyosu::sec )| [secpi](@ref kyosu::secpi ) | + | [sech](@ref kyosu::sech ) | [sign](@ref kyosu::sign ) | [sin](@ref kyosu::sin ) | [sinc](@ref kyosu::sinc ) | [sincos](@ref kyosu::sincos ) | + | [sinpi](@ref kyosu::sinpi ) | [sinpicospi](@ref kyosu::sinpicospi ) | [sinh](@ref kyosu::sinh ) | [sinhcosh](@ref kyosu::sinhcosh ) | [slerp](@ref kyosu::slerp ) | + | [sqr](@ref kyosu::sqr ) | [sqr_abs](@ref kyosu::sqr_abs ) | [sqrt](@ref kyosu::sqrt ) | [tan](@ref kyosu::tan ) | [tanpi](@ref kyosu::tanpi ) | + | [tanh](@ref kyosu::tanh ) | [to_polar](@ref kyosu::to_polar ) | [tgamma](@ref kyosu::tgamma ) | [trunc](@ref kyosu::trunc ) | [zeta](@ref kyosu::zeta ) | + * callables usable with quaternion complex and real only diff --git a/include/kyosu/details/cayleyify.hpp b/include/kyosu/details/cayleyify.hpp new file mode 100644 index 00000000..8ae31006 --- /dev/null +++ b/include/kyosu/details/cayleyify.hpp @@ -0,0 +1,34 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +auto cayley_extend(auto f, auto z, auto... rs) +{ + auto p = kyosu::pure(z); + auto az = kyosu::abs(p); + auto c = f(kyosu::complex(kyosu::real(z), az), rs...); + return kyosu::real(c) + kyosu::ipart(c)*kyosu::sign(p); +} + +auto cayley_extend_rev(auto f, auto z1, auto z2) +{ + auto p = kyosu::pure(z2); + auto az2 = kyosu::abs(p); + auto c = f(z1, kyosu::complex(kyosu::real(z2), az2)); + return kyosu::real(c) + kyosu::ipart(c)*kyosu::sign(p); +} + +auto cayley_extend2(auto f, auto z, auto... rs) +{ + auto p = kyosu::pure(z); + auto az = kyosu::abs(p); + auto [c1, c2] = f(kyosu::complex(kyosu::real(z), az), rs...); + auto s = kyosu::sign(p); + return kumi::tuple{ kyosu::real(c1) + kyosu::ipart(c1)*s + , kyosu::real(c2) + kyosu::ipart(c2)*s}; +} diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 4a3cfbeb..5ab94fe5 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -47,8 +48,14 @@ #include #include #include +#include +#include #include #include +#include +#include +#include +#include #include #include #include @@ -57,6 +64,8 @@ #include #include #include +#include + #include #include #include @@ -90,18 +99,24 @@ #include #include #include +#include +#include #include #include #include #include #include #include +#include +#include #include #include #include #include #include #include +#include +#include #include #include #include @@ -143,35 +158,15 @@ #include #include #include +#include #include #include - -//====================================================================================================================== -//! @brief Functions performing computations over complex or real elements only. -//====================================================================================================================== - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include //====================================================================================================================== //! @brief Functions performing computations over quaternion complex or real elements only. //====================================================================================================================== - #include #include #include diff --git a/include/kyosu/complex/beta.hpp b/include/kyosu/functions/beta.hpp similarity index 90% rename from include/kyosu/complex/beta.hpp rename to include/kyosu/functions/beta.hpp index 12813f9a..51ac0869 100644 --- a/include/kyosu/complex/beta.hpp +++ b/include/kyosu/functions/beta.hpp @@ -61,12 +61,12 @@ namespace kyosu //! //! **Parameters** //! -//! * `x`,`y` : Values to process. Can be a mix of complex and real floating values and complex values. +//! * `x`,`y` : Values to process. Can be a mix of cayley_dickson and real floating values. //! //! **Return value** //! //! 1. If x and y are real typed values returns \f$\displaystyle \mathbf{B}(x,y) = \int_0^1 t^{x-1}(1-t)^{y-1}\mbox{d}t\f$ -//! 2. if x or y is complex the value \f$\displaystyle \mathbf{B}(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\f$ is returned. +//! 2. \f$\displaystyle \mathbf{B}(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\f$ is returned. //! //! @groupheader{Example} //! diff --git a/include/kyosu/complex/deta.hpp b/include/kyosu/functions/deta.hpp similarity index 90% rename from include/kyosu/complex/deta.hpp rename to include/kyosu/functions/deta.hpp index 4768962f..532ebf6e 100644 --- a/include/kyosu/complex/deta.hpp +++ b/include/kyosu/functions/deta.hpp @@ -56,8 +56,7 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto deta(K k, T z) noexcept; //1 -//! template constexpr auto deta(K k, T z) noexcept; //2 +//! constexpr auto deta(unsigned_scalar_value auto k, auto z) noexcept; //! } //! @endcode //! diff --git a/include/kyosu/complex/digamma.hpp b/include/kyosu/functions/digamma.hpp similarity index 90% rename from include/kyosu/complex/digamma.hpp rename to include/kyosu/functions/digamma.hpp index eea3ecd6..a5436ad5 100644 --- a/include/kyosu/complex/digamma.hpp +++ b/include/kyosu/functions/digamma.hpp @@ -52,14 +52,13 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr T digamma(T z) noexcept; -//! template constexpr T digamma(T z) noexcept; +//! constexpr auto digamma(auto z) noexcept; //! } //! @endcode //! //! **Parameters** //! -//! * `z` : real or complex value to process. +//! * `z` : value to process. //! //! **Return value** //! diff --git a/include/kyosu/complex/erf.hpp b/include/kyosu/functions/erf.hpp similarity index 85% rename from include/kyosu/complex/erf.hpp rename to include/kyosu/functions/erf.hpp index fd88ff58..cc22ce24 100644 --- a/include/kyosu/complex/erf.hpp +++ b/include/kyosu/functions/erf.hpp @@ -44,7 +44,7 @@ namespace kyosu //! @var erf //! @brief Computes the error function: \f$ \displaystyle //! \mbox{erf}(x)=\frac{2}{\sqrt\pi}\int_0^{x} e^{-t^2}\mbox{d}t\f$ or -//! its analytic continuation in the complex plane +//! its extension to complex and general cayley-dickson values //! //! **Defined in Header** //! @@ -57,8 +57,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto erf(T z) noexcept; //1 -//! template constexpr auto erf(T z) noexcept; //2 +//! { template constexpr T erf(T z) noexcept; +//! template constexpr T erf(T z) noexcept; //! } //! @endcode //! @@ -68,9 +68,7 @@ namespace kyosu //! //! **Return value** //! -//! 1. a real input z returns eve::erf(z). -//! -//! 2. The value of the error function in the complex plane is returned +//! * The value of the error function is returned. //! //! @groupheader{Example} //! diff --git a/include/kyosu/complex/erfcx.hpp b/include/kyosu/functions/erfcx.hpp similarity index 91% rename from include/kyosu/complex/erfcx.hpp rename to include/kyosu/functions/erfcx.hpp index abdd8f01..da3cb493 100644 --- a/include/kyosu/complex/erfcx.hpp +++ b/include/kyosu/functions/erfcx.hpp @@ -56,8 +56,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto erfcx(T z) noexcept; //1 -//! template constexpr auto erfcx(T z) noexcept; //2 +//! template constexpr auto erfcx(T z) noexcept; //1 +//! template constexpr auto erfcx(T z) noexcept; //2 //! } //! @endcode //! diff --git a/include/kyosu/complex/erfi.hpp b/include/kyosu/functions/erfi.hpp similarity index 88% rename from include/kyosu/complex/erfi.hpp rename to include/kyosu/functions/erfi.hpp index 5390a945..6a53ab03 100644 --- a/include/kyosu/complex/erfi.hpp +++ b/include/kyosu/functions/erfi.hpp @@ -23,7 +23,7 @@ namespace kyosu::tags { auto over = eve::sqr(v) > 720; auto r = eve::inf(eve::as(v))*eve::sign(v); - r = eve::if_else(over, r, -kyosu::imag(kyosu::erf(complex(eve::zero(eve::as(v)), -v)))); + r = eve::if_else(over, r, -get<1>(kyosu::erf(complex(eve::zero(eve::as(v)), -v)))); return complex(r); } @@ -58,8 +58,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto erfi(T z) noexcept; //1 -//! template constexpr auto erfi(T z) noexcept; //2 +//! template constexpr auto erfi(T z) noexcept; //1 +//! template constexpr auto erfi(T z) noexcept; //2 //! } //! @endcode //! diff --git a/include/kyosu/complex/eta.hpp b/include/kyosu/functions/eta.hpp similarity index 91% rename from include/kyosu/complex/eta.hpp rename to include/kyosu/functions/eta.hpp index 1ddff424..de593cd2 100644 --- a/include/kyosu/complex/eta.hpp +++ b/include/kyosu/functions/eta.hpp @@ -57,14 +57,14 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto eta(K, k, T z) noexcept; //1 -//! template constexpr auto eta(K, k, T z) noexcept; //2 +//! template constexpr auto eta(T z) noexcept; //1 +//! template constexpr auto eta(T z) noexcept; //2 //! } //! @endcode //! //! **Parameters** //! -//! * `z` : complex or real value to process. +//! * `z` : value to process. //! //! **Return value** //! diff --git a/include/kyosu/complex/faddeeva.hpp b/include/kyosu/functions/faddeeva.hpp similarity index 92% rename from include/kyosu/complex/faddeeva.hpp rename to include/kyosu/functions/faddeeva.hpp index d26f81d5..754ff9e3 100644 --- a/include/kyosu/complex/faddeeva.hpp +++ b/include/kyosu/functions/faddeeva.hpp @@ -56,8 +56,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto faddeeva(T z) noexcept; //1 -//! template constexpr auto faddeeva(T z) noexcept; //2 +//! template constexpr auto faddeeva(T z) noexcept; +//! template constexpr auto faddeeva(T z) noexcept; //! } //! @endcode //! diff --git a/include/kyosu/functions/is_imag.hpp b/include/kyosu/functions/is_imag.hpp index 3cd375d6..97a8f2e5 100644 --- a/include/kyosu/functions/is_imag.hpp +++ b/include/kyosu/functions/is_imag.hpp @@ -28,8 +28,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto is_imag(T z) noexcept; -//! template constexpr auto is_imag(T z) noexcept; +//! template constexpr auto is_imag(T z) noexcept; +//! template constexpr auto is_imag(T z) noexcept; //! } //! @endcode //! diff --git a/include/kyosu/complex/lambda.hpp b/include/kyosu/functions/lambda.hpp similarity index 90% rename from include/kyosu/complex/lambda.hpp rename to include/kyosu/functions/lambda.hpp index ce0dd606..649124d6 100644 --- a/include/kyosu/complex/lambda.hpp +++ b/include/kyosu/functions/lambda.hpp @@ -48,6 +48,7 @@ namespace kyosu //! This function can be extended to the whole complex plane as \f$\lambda(z) = \zeta(z)(1-2^{-x})\f$ //! (where \f$\zeta\f$ is the Riemann zeta function). It coincides with the serie where the serie converges. //! However for `z = 1` the result is \f$\infty\f$. +//! The usual extension mechanism is used for general Cayleyâ-dickson values. //! //! **Defined in Header** //! @@ -60,8 +61,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto lambda(K, k, T z) noexcept; //1 -//! template constexpr auto lambda(K, k, T z) noexcept; //2 +//! template constexpr auto lambda(T z) noexcept; //1 +//! template constexpr auto lambda(T z) noexcept; //2 //! } //! @endcode //! diff --git a/include/kyosu/complex/lbeta.hpp b/include/kyosu/functions/lbeta.hpp similarity index 77% rename from include/kyosu/complex/lbeta.hpp rename to include/kyosu/functions/lbeta.hpp index f343bd3a..c67fe289 100644 --- a/include/kyosu/complex/lbeta.hpp +++ b/include/kyosu/functions/lbeta.hpp @@ -53,17 +53,7 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template< eve::floating_ordered_value T, eve::floating_ordered_value U > -//! auto lbeta(T x,U y) noexcept; //1 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto lbeta(eve::complex_t x, U y) noexcept; //2 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto lbeta(T x, eve::complex_t y) noexcept; //2 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto lbeta(eve::complex_t x, eve::complex_t y) noexcept; //2 +//! auto lbeta(auto x,auto y) noexcept; //! } //! @endcode //! diff --git a/include/kyosu/complex/log_abs_gamma.hpp b/include/kyosu/functions/log_abs_gamma.hpp similarity index 91% rename from include/kyosu/complex/log_abs_gamma.hpp rename to include/kyosu/functions/log_abs_gamma.hpp index 5e6fc7c1..08ebcc6c 100644 --- a/include/kyosu/complex/log_abs_gamma.hpp +++ b/include/kyosu/functions/log_abs_gamma.hpp @@ -55,8 +55,7 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr T log_abs_gamma(T z) noexcept; -//! template constexpr as_real_t log_abs_gamma(T z) noexcept; +//! constexpr auto log_abs_gamma(auto z) noexcept; //! } //! @endcode //! diff --git a/include/kyosu/complex/log_gamma.hpp b/include/kyosu/functions/log_gamma.hpp similarity index 88% rename from include/kyosu/complex/log_gamma.hpp rename to include/kyosu/functions/log_gamma.hpp index b9955865..817f15ca 100644 --- a/include/kyosu/complex/log_gamma.hpp +++ b/include/kyosu/functions/log_gamma.hpp @@ -40,7 +40,7 @@ namespace kyosu //! @addtogroup functions //! @{ //! @var log_gamma -//! @brief Computes the log of the \f$\Gamma\f$ function of the parameter. +//! @brief Computes the log of the \f$\Gamma\f$ function. //! //! **Defined in Header** //! @@ -53,8 +53,7 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto log_gamma(T z) noexcept; -//! template constexpr auto log_gamma(T z) noexcept; +//! constexpr auto log_gamma(T z) noexcept; //! } //! @endcode //! diff --git a/include/kyosu/complex/lrising_factorial.hpp b/include/kyosu/functions/lrising_factorial.hpp similarity index 74% rename from include/kyosu/complex/lrising_factorial.hpp rename to include/kyosu/functions/lrising_factorial.hpp index 287beec5..008cb69f 100644 --- a/include/kyosu/complex/lrising_factorial.hpp +++ b/include/kyosu/functions/lrising_factorial.hpp @@ -53,17 +53,7 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template< eve::floating_ordered_value T, eve::floating_ordered_value U > -//! auto lrising_factorial(T x,U y) noexcept; //1 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto lrising_factorial(eve::complex_t x, U y) noexcept; //2 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto lrising_factorial(T x, eve::complex_t y) noexcept; //2 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto lrising_factorial(eve::complex_t x, eve::complex_t y) noexcept; //2 +//! auto lrising_factorial(auto x,auto y) noexcept; //! } //! @endcode //! @@ -73,7 +63,7 @@ namespace kyosu //! //! **Return value** //! -//! @brief Computes the logarithm Rising Factorial function i.e. \f$\log\frac{\Gamma(x+y)}{\Gamma(x)}\f$. +//! @brief Computes the logarithm of rising Factorial i.e. \f$\log\frac{\Gamma(x+y)}{\Gamma(x)}\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/complex/rising_factorial.hpp b/include/kyosu/functions/rising_factorial.hpp similarity index 74% rename from include/kyosu/complex/rising_factorial.hpp rename to include/kyosu/functions/rising_factorial.hpp index d3b9d044..4ece00a4 100644 --- a/include/kyosu/complex/rising_factorial.hpp +++ b/include/kyosu/functions/rising_factorial.hpp @@ -53,17 +53,7 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template< eve::floating_ordered_value T, eve::floating_ordered_value U > -//! auto rising_factorial(T x,U y) noexcept; //1 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto rising_factorial(eve::complex_t x, U y) noexcept; //2 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto rising_factorial(T x, eve::complex_t y) noexcept; //2 -//! -//! template< eve::floating_value T, eve::floating_value U > -//! auto rising_factorial(eve::complex_t x, eve::complex_t y) noexcept; //2 +//! constexpr auto rising_factorial(auto x,auto y) noexcept; //! } //! @endcode //! @@ -73,7 +63,7 @@ namespace kyosu //! //! **Return value** //! -//! @brief Computes the Rising Factorial function i.e. \f$\frac{\Gamma(x+y)}{\Gamma(x)}\f$. +//! @brief Computes the rising Factorial i.e. \f$\frac{\Gamma(x+y)}{\Gamma(x)}\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/complex/tgamma.hpp b/include/kyosu/functions/tgamma.hpp similarity index 92% rename from include/kyosu/complex/tgamma.hpp rename to include/kyosu/functions/tgamma.hpp index ac51cedf..234fb25a 100644 --- a/include/kyosu/complex/tgamma.hpp +++ b/include/kyosu/functions/tgamma.hpp @@ -54,8 +54,7 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr T tgamma(T z) noexcept; -//! template constexpr T tgamma(T z) noexcept; +//! constexpr T tgamma(T z) noexcept; //! } //! @endcode //! diff --git a/include/kyosu/complex/zeta.hpp b/include/kyosu/functions/zeta.hpp similarity index 89% rename from include/kyosu/complex/zeta.hpp rename to include/kyosu/functions/zeta.hpp index 2115450f..e7886e85 100644 --- a/include/kyosu/complex/zeta.hpp +++ b/include/kyosu/functions/zeta.hpp @@ -57,17 +57,17 @@ namespace kyosu //! namespace kyosu //! { //! template constexpr auto zeta(T z) noexcept; -//! template constexpr auto zeta(T z) noexcept; +//! template constexpr auto zeta(T z) noexcept; //! } //! @endcode //! //! **Parameters** //! -//! * `z` : complex or real value to process. +//! * `z` : value to process. //! //! **Return value** //! -//! Returns the Dirichlet zeta function: \f$ \displaystyle \sum_0^\infty \frac{1}{(n+1)^z}\f$ +//! Returns the Dirichlet zeta sum: \f$ \displaystyle \sum_0^\infty \frac{1}{(n+1)^z}\f$ //! //! @groupheader{Example} //! diff --git a/include/kyosu/types/cayley_dickson.hpp b/include/kyosu/types/cayley_dickson.hpp index c336d7a6..fff70f4c 100644 --- a/include/kyosu/types/cayley_dickson.hpp +++ b/include/kyosu/types/cayley_dickson.hpp @@ -16,8 +16,8 @@ #include #include #include -#include -#include +#include +#include #include #include diff --git a/include/kyosu/types/impl/complex/erf.hpp b/include/kyosu/types/impl/complex/erf.hpp deleted file mode 100644 index f1b8bb10..00000000 --- a/include/kyosu/types/impl/complex/erf.hpp +++ /dev/null @@ -1,290 +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::_ -{ - template - auto dispatch(eve::tag_of, Z z) noexcept - { - auto numeric_sqr = [](auto z){ // UNTIL DECORATORS ARE AT HAND - auto [zr, zi] = z; - return complex((zr-zi)*(zi+zr), 2 * zr * zi); - }; - auto x = real(z); - auto y = imag(z); - using real_t = decltype(x); - using r_t = eve::element_type_t; - auto mz2 = -numeric_sqr(z);// -z^2, being careful of overflow - using A9 = kumi::result::generate_t<9, r_t>; - static constexpr std::array< A9, 97> coefs = - { - A9{0.28351593328822191546e-2,0.28494783221378400759e-2,0.14427470563276734183e-4,0.10939723080231588129e-6,0.92474307943275042045e-9,0.89128907666450075245e-11,0.92974121935111111110e-13,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.85927161243940350562e-2,0.29085312941641339862e-2,0.15106783707725582090e-4,0.11716709978531327367e-6,0.10197387816021040024e-8,0.10122678863073360769e-10,0.10917479678400000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.14471159831187703054e-1,0.29703978970263836210e-2,0.15835096760173030976e-4,0.12574803383199211596e-6,0.11278672159518415848e-8,0.11547462300333495797e-10,0.12894535335111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.20476320420324610618e-1,0.30352843012898665856e-2,0.16617609387003727409e-4,0.13525429711163116103e-6,0.12515095552507169013e-8,0.13235687543603382345e-10,0.15326595042666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.26614461952489004566e-1,0.31034189276234947088e-2,0.17460268109986214274e-4,0.14582130824485709573e-6,0.13935959083809746345e-8,0.15249438072998932900e-10,0.18344741882133333333e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.32892330248093586215e-1,0.31750557067975068584e-2,0.18369907582308672632e-4,0.15761063702089457882e-6,0.15577638230480894382e-8,0.17663868462699097951e-10,0.22126732680711111111e-12,0.30273474177737853668e-14,0.00000000000000000000e+00}, - A9{0.39317207681134336024e-1,0.32504779701937539333e-2,0.19354426046513400534e-4,0.17081646971321290539e-6,0.17485733959327106250e-8,0.20593687304921961410e-10,0.26917401949155555556e-12,0.38562123837725712270e-14,0.00000000000000000000e+00}, - A9{0.45896976511367738235e-1,0.33300031273110976165e-2,0.20423005398039037313e-4,0.18567412470376467303e-6,0.19718038363586588213e-8,0.24175006536781219807e-10,0.33059982791466666666e-12,0.49756574284439426165e-14,0.00000000000000000000e+00}, - A9{0.52640192524848962855e-1,0.34139883358846720806e-2,0.21586390240603337337e-4,0.20247136501568904646e-6,0.22348696948197102935e-8,0.28597516301950162548e-10,0.41045502119111111110e-12,0.65151614515238361946e-14,0.00000000000000000000e+00}, - A9{0.59556171228656770456e-1,0.35028374386648914444e-2,0.22857246150998562824e-4,0.22156372146525190679e-6,0.25474171590893813583e-8,0.34122390890697400584e-10,0.51593189879111111110e-12,0.86775076853908006938e-14,0.00000000000000000000e+00}, - A9{0.66655089485108212551e-1,0.35970095381271285568e-2,0.24250626164318672928e-4,0.24339561521785040536e-6,0.29221990406518411415e-8,0.41117013527967776467e-10,0.65786450716444444445e-12,0.11791885745450623331e-13,0.00000000000000000000e+00}, - A9{0.73948106345519174661e-1,0.36970297216569341748e-2,0.25784588137312868792e-4,0.26853012002366752770e-6,0.33763958861206729592e-8,0.50111549981376976397e-10,0.85313857496888888890e-12,0.16417079927706899860e-13,0.00000000000000000000e+00}, - A9{0.81447508065002963203e-1,0.38035026606492705117e-2,0.27481027572231851896e-4,0.29769200731832331364e-6,0.39336816287457655076e-8,0.61895471132038157624e-10,0.11292303213511111111e-11,0.23558532213703884304e-13,0.00000000000000000000e+00}, - A9{0.89166884027582716628e-1,0.39171301322438946014e-2,0.29366827260422311668e-4,0.33183204390350724895e-6,0.46276006281647330524e-8,0.77692631378169813324e-10,0.15335153258844444444e-11,0.35183103415916026911e-13,0.00000000000000000000e+00}, - A9{0.97121342888032322019e-1,0.40387340353207909514e-2,0.31475490395950776930e-4,0.37222714227125135042e-6,0.55074373178613809996e-8,0.99509175283990337944e-10,0.21552645758222222222e-11,0.55728651431872687605e-13,0.00000000000000000000e+00}, - A9{0.10532778218603311137e00,0.41692873614065380607e-2,0.33849549774889456984e-4,0.42064596193692630143e-6,0.66494579697622432987e-8,0.13094103581931802337e-09,0.31896187409777777778e-11,0.97271974184476560742e-13,0.00000000000000000000e+00}, - A9{0.11380523107427108222e00,0.43099572287871821013e-2,0.36544324341565929930e-4,0.47965044028581857764e-6,0.81819034238463698796e-8,0.17934133239549647357e-09,0.50956666166186293627e-11,0.18850487318190638010e-12,0.79697813173519853340e-14}, - A9{0.12257529703447467345e00,0.44621675710026986366e-2,0.39634304721292440285e-4,0.55321553769873381819e-6,0.10343619428848520870e-7,0.26033830170470368088e-09,0.87743837749108025357e-11,0.34427092430230063401e-12,0.10205506615709843189e-13}, - A9{0.13166276955656699478e00,0.46276970481783001803e-2,0.43225026380496399310e-4,0.64799164020016902656e-6,0.13580082794704641782e-7,0.39839800853954313927e-09,0.14431142411840000000e-10,0.42193457308830027541e-12,0.00000000000000000000e+00}, - A9{0.14109647869803356475e00,0.48088424418545347758e-2,0.47474504753352150205e-4,0.77509866468724360352e-6,0.18536851570794291724e-7,0.60146623257887570439e-09,0.18533978397305276318e-10,0.41033845938901048380e-13,-0.46160680279304825485e-13}, - A9{0.15091057940548936603e00,0.50086864672004685703e-2,0.52622482832192230762e-4,0.95034664722040355212e-6,0.25614261331144718769e-7,0.80183196716888606252e-09,0.12282524750534352272e-10,-0.10531774117332273617e-11,-0.86157181395039646412e-13}, - A9{0.16114648116017010770e00,0.52314661581655369795e-2,0.59005534545908331315e-4,0.11885518333915387760e-5,0.33975801443239949256e-7,0.82111547144080388610e-09,-0.12357674017312854138e-10,-0.24355112256914479176e-11,-0.75155506863572930844e-13}, - A9{0.17185551279680451144e00,0.54829002967599420860e-2,0.67013226658738082118e-4,0.14897400671425088807e-5,0.40690283917126153701e-7,0.44060872913473778318e-09,-0.52641873433280000000e-10,-0.30940587864543343124e-11,0.00000000000000000000e+00}, - A9{0.18310194559815257381e00,0.57701559375966953174e-2,0.76948789401735193483e-4,0.18227569842290822512e-5,0.41092208344387212276e-7,-0.44009499965694442143e-09,-0.92195414685628803451e-10,-0.22657389705721753299e-11,0.10004784908106839254e-12}, - A9{0.19496527191546630345e00,0.61010853144364724856e-2,0.88812881056342004864e-4,0.21180686746360261031e-5,0.30652145555130049203e-7,-0.16841328574105890409e-08,-0.11008129460612823934e-09,-0.12180794204544515779e-12,0.15703325634590334097e-12}, - A9{0.20754006813966575720e00,0.64825787724922073908e-2,0.10209599627522311893e-3,0.22785233392557600468e-5,0.73495224449907568402e-8,-0.29442705974150112783e-08,-0.94082603434315016546e-10,0.23609990400179321267e-11,0.14141908654269023788e-12}, - A9{0.22093185554845172146e00,0.69182878150187964499e-2,0.11568723331156335712e-3,0.22060577946323627739e-5,-0.26929730679360840096e-7,-0.38176506152362058013e-08,-0.47399503861054459243e-10,0.40953700187172127264e-11,0.69157730376118511127e-13}, - A9{0.23524827304057813918e00,0.74063350762008734520e-2,0.12796333874615790348e-3,0.18327267316171054273e-5,-0.66742910737957100098e-7,-0.40204740975496797870e-08,0.14515984139495745330e-10,0.44921608954536047975e-11,-0.18583341338983776219e-13}, - A9{0.25058626331812744775e00,0.79377285151602061328e-2,0.13704268650417478346e-3,0.11427511739544695861e-5,-0.10485442447768377485e-6,-0.34850364756499369763e-08,0.72656453829502179208e-10,0.36195460197779299406e-11,-0.84882136022200714710e-13}, - A9{0.26701724900280689785e00,0.84959936119625864274e-2,0.14112359443938883232e-3,0.17800427288596909634e-6,-0.13443492107643109071e-6,-0.23512456315677680293e-08,0.11245846264695936769e-09,0.19850501334649565404e-11,-0.11284666134635050832e-12}, - A9{0.28457293586253654144e00,0.90581563892650431899e-2,0.13880520331140646738e-3,-0.97262302362522896157e-6,-0.15077100040254187366e-6,-0.88574317464577116689e-09,0.12760311125637474581e-09,0.20155151018282695055e-12,-0.10514169375181734921e-12}, - A9{0.30323425595617385705e00,0.95968346790597422934e-2,0.12931067776725883939e-3,-0.21938741702795543986e-5,-0.15202888584907373963e-6,0.61788350541116331411e-09,0.11957835742791248256e-09,-0.12598179834007710908e-11,-0.75151817129574614194e-13}, - A9{0.32292521181517384379e00,0.10082957727001199408e-1,0.11257589426154962226e-3,-0.33670890319327881129e-5,-0.13910529040004008158e-6,0.19170714373047512945e-08,0.94840222377720494290e-10,-0.21650018351795353201e-11,-0.37875211678024922689e-13}, - A9{0.34351233557911753862e00,0.10488575435572745309e-1,0.89209444197248726614e-4,-0.43893459576483345364e-5,-0.11488595830450424419e-6,0.28599494117122464806e-08,0.61537542799857777779e-10,-0.24935749227658002212e-11,0.00000000000000000000e+00}, - A9{0.36480946642143669093e00,0.10789304203431861366e-1,0.60357993745283076834e-4,-0.51855862174130669389e-5,-0.83291664087289801313e-7,0.33898011178582671546e-08,0.27082948188277716482e-10,-0.23603379397408694974e-11,0.19328087692252869842e-13}, - A9{0.38658679935694939199e00,0.10966119158288804999e-1,0.27521612041849561426e-4,-0.57132774537670953638e-5,-0.48404772799207914899e-7,0.35268354132474570493e-08,-0.32383477652514618094e-11,-0.19334202915190442501e-11,0.32333189861286460270e-13}, - A9{0.40858275583808707870e00,0.11006378016848466550e-1,-0.76396376685213286033e-5,-0.59609835484245791439e-5,-0.13834610033859313213e-7,0.33406952974861448790e-08,-0.26474915974296612559e-10,-0.13750229270354351983e-11,0.36169366979417390637e-13}, - A9{0.43051714914006682977e00,0.10904106549500816155e-1,-0.43477527256787216909e-4,-0.59429739547798343948e-5,0.17639200194091885949e-7,0.29235991689639918688e-08,-0.41718791216277812879e-10,-0.81023337739508049606e-12,0.33618915934461994428e-13}, - A9{0.45210428135559607406e00,0.10659670756384400554e-1,-0.78488639913256978087e-4,-0.56919860886214735936e-5,0.44181850467477733407e-7,0.23694306174312688151e-08,-0.49492621596685443247e-10,-0.31827275712126287222e-12,0.27494438742721623654e-13}, - A9{0.47306491195005224077e00,0.10279006119745977570e-1,-0.11140268171830478306e-3,-0.52518035247451432069e-5,0.64846898158889479518e-7,0.17603624837787337662e-08,-0.51129481592926104316e-10,0.62674584974141049511e-13,0.20055478560829935356e-13}, - A9{0.49313638965719857647e00,0.97725799114772017662e-2,-0.14122854267291533334e-3,-0.46707252568834951907e-5,0.79421347979319449524e-7,0.11603027184324708643e-08,-0.48269605844397175946e-10,0.32477251431748571219e-12,0.12831052634143527985e-13}, - A9{0.51208057433416004042e00,0.91542422354009224951e-2,-0.16726530230228647275e-3,-0.39964621752527649409e-5,0.88232252903213171454e-7,0.61343113364949928501e-09,-0.42516755603130443051e-10,0.47910437172240209262e-12,0.66784341874437478953e-14}, - A9{0.52968945458607484524e00,0.84400880445116786088e-2,-0.18908729783854258774e-3,-0.32725905467782951931e-5,0.91956190588652090659e-7,0.14593989152420122909e-09,-0.35239490687644444445e-10,0.54613829888448694898e-12,0.00000000000000000000e+00}, - A9{0.54578857454330070965e00,0.76474155195880295311e-2,-0.20651230590808213884e-3,-0.25364339140543131706e-5,0.91455367999510681979e-7,-0.23061359005297528898e-09,-0.27512928625244444444e-10,0.54895806008493285579e-12,0.00000000000000000000e+00}, - A9{0.56023851910298493910e00,0.67938321739997196804e-2,-0.21956066613331411760e-3,-0.18181127670443266395e-5,0.87650335075416845987e-7,-0.51548062050366615977e-09,-0.20068462174044444444e-10,0.50912654909758187264e-12,0.00000000000000000000e+00}, - A9{0.57293478057455721150e00,0.58965321010394044087e-2,-0.22841145229276575597e-3,-0.11404605562013443659e-5,0.81430290992322326296e-7,-0.71512447242755357629e-09,-0.13372664928000000000e-10,0.44461498336689298148e-12,0.00000000000000000000e+00}, - A9{0.58380635448407827360e00,0.49717469530842831182e-2,-0.23336001540009645365e-3,-0.51952064448608850822e-6,0.73596577815411080511e-7,-0.84020916763091566035e-09,-0.76700972702222222221e-11,0.36914462807972467044e-12,0.00000000000000000000e+00}, - A9{0.59281340237769489597e00,0.40343592069379730568e-2,-0.23477963738658326185e-3,0.34615944987790224234e-7,0.64832803248395814574e-7,-0.90329163587627007971e-09,-0.30421940400000000000e-11,0.29237386653743536669e-12,0.00000000000000000000e+00}, - A9{0.59994428743114271918e00,0.30976579788271744329e-2,-0.23308875765700082835e-3,0.51681681023846925160e-6,0.55694594264948268169e-7,-0.91719117313243464652e-09,0.53982743680000000000e-12,0.22050829296187771142e-12,0.00000000000000000000e+00}, - A9{0.60521224471819875444e00,0.21732138012345456060e-2,-0.22872428969625997456e-3,0.92588959922653404233e-6,0.46612665806531930684e-7,-0.89393722514414153351e-09,0.31718550353777777778e-11,0.15705458816080549117e-12,0.00000000000000000000e+00}, - A9{0.60865189969791123620e00,0.12708480848877451719e-2,-0.22212090111534847166e-3,0.12636236031532793467e-5,0.37904037100232937574e-7,-0.84417089968101223519e-09,0.49843180828444444445e-11,0.10355439441049048273e-12,0.00000000000000000000e+00}, - A9{0.61031580103499200191e00,0.39867436055861038223e-3,-0.21369573439579869291e-3,0.15339402129026183670e-5,0.29787479206646594442e-7,-0.77687792914228632974e-09,0.61192452741333333334e-11,0.60216691829459295780e-13,0.00000000000000000000e+00}, - A9{0.61027109047879835868e00,-0.43680904508059878254e-3,-0.20383783788303894442e-3,0.17421743090883439959e-5,0.22400425572175715576e-7,-0.69934719320045128997e-09,0.67152759655111111110e-11,0.26419960042578359995e-13,0.00000000000000000000e+00}, - A9{0.60859639489217430521e00,-0.12305921390962936873e-2,-0.19290150253894682629e-3,0.18944904654478310128e-5,0.15815530398618149110e-7,-0.61726850580964876070e-09,0.68987888999111111110e-11,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.60537899426486075181e00,-0.19790062241395705751e-2,-0.18120271393047062253e-3,0.19974264162313241405e-5,0.10055795094298172492e-7,-0.53491997919318263593e-09,0.67794550295111111110e-11,-0.17059208095741511603e-13,0.00000000000000000000e+00}, - A9{0.60071229457904110537e00,-0.26795676776166354354e-2,-0.16901799553627508781e-3,0.20575498324332621581e-5,0.51077165074461745053e-8,-0.45536079828057221858e-09,0.64488005516444444445e-11,-0.29311677573152766338e-13,0.00000000000000000000e+00}, - A9{0.59469361520112714738e00,-0.33308208190600993470e-2,-0.15658501295912405679e-3,0.20812116912895417272e-5,0.93227468760614182021e-9,-0.38066673740116080415e-09,0.59806790359111111110e-11,-0.36887077278950440597e-13,0.00000000000000000000e+00}, - A9{0.58742228631775388268e00,-0.39321858196059227251e-2,-0.14410441141450122535e-3,0.20743790018404020716e-5,-0.25261903811221913762e-8,-0.31212416519526924318e-09,0.54328422462222222221e-11,-0.40864152484979815972e-13,0.00000000000000000000e+00}, - A9{0.57899804200033018447e00,-0.44838157005618913447e-2,-0.13174245966501437965e-3,0.20425306888294362674e-5,-0.53330296023875447782e-8,-0.25041289435539821014e-09,0.48490437205333333334e-11,-0.42162206939169045177e-13,0.00000000000000000000e+00}, - A9{0.56951968796931245974e00,-0.49864649488074868952e-2,-0.11963416583477567125e-3,0.19906021780991036425e-5,-0.75580140299436494248e-8,-0.19576060961919820491e-09,0.42613011928888888890e-11,-0.41539443304115604377e-13,0.00000000000000000000e+00}, - A9{0.55908401930063918964e00,-0.54413711036826877753e-2,-0.10788661102511914628e-3,0.19229663322982839331e-5,-0.92714731195118129616e-8,-0.14807038677197394186e-09,0.36920870298666666666e-11,-0.39603726688419162617e-13,0.00000000000000000000e+00}, - A9{0.54778496152925675315e00,-0.58501497933213396670e-2,-0.96582314317855227421e-4,0.18434405235069270228e-5,-0.10541580254317078711e-7,-0.10702303407788943498e-09,0.31563175582222222222e-11,-0.36829748079110481422e-13,0.00000000000000000000e+00}, - A9{0.53571290831682823999e00,-0.62147030670760791791e-2,-0.85782497917111760790e-4,0.17553116363443470478e-5,-0.11432547349815541084e-7,-0.72157091369041330520e-10,0.26630811607111111111e-11,-0.33578660425893164084e-13,0.00000000000000000000e+00}, - A9{0.52295422962048434978e00,-0.65371404367776320720e-2,-0.75530164941473343780e-4,0.16613725797181276790e-5,-0.12003521296598910761e-7,-0.42929753689181106171e-10,0.22170894940444444444e-11,-0.30117697501065110505e-13,0.00000000000000000000e+00}, - A9{0.50959092577577886140e00,-0.68197117603118591766e-2,-0.65852936198953623307e-4,0.15639654113906716939e-5,-0.12308007991056524902e-7,-0.18761997536910939570e-10,0.18198628922666666667e-11,-0.26638355362285200932e-13,0.00000000000000000000e+00}, - A9{0.49570040481823167970e00,-0.70647509397614398066e-2,-0.56765617728962588218e-4,0.14650274449141448497e-5,-0.12393681471984051132e-7,0.92904351801168955424e-12,0.14706755960177777778e-11,-0.23272455351266325318e-13,0.00000000000000000000e+00}, - A9{0.48135536250935238066e00,-0.72746293327402359783e-2,-0.48272489495730030780e-4,0.13661377309113939689e-5,-0.12302464447599382189e-7,0.16707760028737074907e-10,0.11672928324444444444e-11,-0.20105801424709924499e-13,0.00000000000000000000e+00}, - A9{0.46662374675511439448e00,-0.74517177649528487002e-2,-0.40369318744279128718e-4,0.12685621118898535407e-5,-0.12070791463315156250e-7,0.29105507892605823871e-10,0.90653314645333333334e-12,-0.17189503312102982646e-13,0.00000000000000000000e+00}, - A9{0.45156879030168268778e00,-0.75983560650033817497e-2,-0.33045110380705139759e-4,0.11732956732035040896e-5,-0.11729986947158201869e-7,0.38611905704166441308e-10,0.68468768305777777779e-12,-0.14549134330396754575e-13,0.00000000000000000000e+00}, - A9{0.43624909769330896904e00,-0.77168291040309554679e-2,-0.26283612321339907756e-4,0.10811018836893550820e-5,-0.11306707563739851552e-7,0.45670446788529607380e-10,0.49782492549333333334e-12,-0.12191983967561779442e-13,0.00000000000000000000e+00}, - A9{0.42071877443548481181e00,-0.78093484015052730097e-2,-0.20064596897224934705e-4,0.99254806680671890766e-6,-0.10823412088884741451e-7,0.50677203326904716247e-10,0.34200547594666666666e-12,-0.10112698698356194618e-13,0.00000000000000000000e+00}, - A9{0.40502758809710844280e00,-0.78780384460872937555e-2,-0.14364940764532853112e-4,0.90803709228265217384e-6,-0.10298832847014466907e-7,0.53981671221969478551e-10,0.21342751381333333333e-12,-0.82975901848387729274e-14,0.00000000000000000000e+00}, - A9{0.38922115269731446690e00,-0.79249269708242064120e-2,-0.91595258799106970453e-5,0.82783535102217576495e-6,-0.97484311059617744437e-8,0.55889029041660225629e-10,0.10851981336888888889e-12,-0.67278553237853459757e-14,0.00000000000000000000e+00}, - A9{0.37334112915460307335e00,-0.79519385109223148791e-2,-0.44219833548840469752e-5,0.75209719038240314732e-6,-0.91848251458553190451e-8,0.56663266668051433844e-10,0.23995894257777777778e-13,-0.53819475285389344313e-14,0.00000000000000000000e+00}, - A9{0.35742543583374223085e00,-0.79608906571527956177e-2,-0.12530071050975781198e-6,0.68088605744900552505e-6,-0.86181844090844164075e-8,0.56530784203816176153e-10,-0.43120012248888888890e-13,-0.42372603392496813810e-14,0.00000000000000000000e+00}, - A9{0.34150846431979618536e00,-0.79534924968773806029e-2,0.37576885610891515813e-5,0.61419263633090524326e-6,-0.80565865409945960125e-8,0.55684175248749269411e-10,-0.95486860764444444445e-13,-0.32712946432984510595e-14,0.00000000000000000000e+00}, - A9{0.32562129649136346824e00,-0.79313448067948884309e-2,0.72539159933545300034e-5,0.55195028297415503083e-6,-0.75063365335570475258e-8,0.54281686749699595941e-10,-0.13545424295111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.30979191977078391864e00,-0.78959416264207333695e-2,0.10389774377677210794e-4,0.49404804463196316464e-6,-0.69722488229411164685e-8,0.52469254655951393842e-10,-0.16507860650666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.29404543811214459904e00,-0.78486728990364155356e-2,0.13190885683106990459e-4,0.44034158861387909694e-6,-0.64578942561562616481e-8,0.50354306498006928984e-10,-0.18614473550222222222e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.27840427686253660515e00,-0.77908279176252742013e-2,0.15681928798708548349e-4,0.39066226205099807573e-6,-0.59658144820660420814e-8,0.48030086420373141763e-10,-0.20018995173333333333e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.26288838011163800908e00,-0.77235993576119469018e-2,0.17886516796198660969e-4,0.34482457073472497720e-6,-0.54977066551955420066e-8,0.45572749379147269213e-10,-0.20852924954666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.24751539954181029717e00,-0.76480877165290370975e-2,0.19827114835033977049e-4,0.30263228619976332110e-6,-0.50545814570120129947e-8,0.43043879374212005966e-10,-0.21228012028444444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.23230087411688914593e00,-0.75653060136384041587e-2,0.21524991113020016415e-4,0.26388338542539382413e-6,-0.46368974069671446622e-8,0.40492715758206515307e-10,-0.21238627815111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.21725840021297341931e00,-0.74761846305979730439e-2,0.23000194404129495243e-4,0.22837400135642906796e-6,-0.42446743058417541277e-8,0.37958104071765923728e-10,-0.20963978568888888889e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.20239979200788191491e00,-0.73815761980493466516e-2,0.24271552727631854013e-4,0.19590154043390012843e-6,-0.38775884642456551753e-8,0.35470192372162901168e-10,-0.20470131678222222222e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.18773523211558098962e00,-0.72822604530339834448e-2,0.25356688567841293697e-4,0.16626710297744290016e-6,-0.35350521468015310830e-8,0.33051896213898864306e-10,-0.19811844544000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.17327341258479649442e00,-0.71789490089142761950e-2,0.26272046822383820476e-4,0.13927732375657362345e-6,-0.32162794266956859603e-8,0.30720156036105652035e-10,-0.19034196304000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.15902166648328672043e00,-0.70722899934245504034e-2,0.27032932310132226025e-4,0.11474573347816568279e-6,-0.29203404091754665063e-8,0.28487010262547971859e-10,-0.18174029063111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.14498609036610283865e00,-0.69628725220045029273e-2,0.27653554229160596221e-4,0.92493727167393036470e-7,-0.26462055548683583849e-8,0.26360506250989943739e-10,-0.17261211260444444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.13117165798208050667e00,-0.68512309830281084723e-2,0.28147075431133863774e-4,0.72351212437979583441e-7,-0.23927816200314358570e-8,0.24345469651209833155e-10,-0.16319736960000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.11758232561160626306e00,-0.67378491192463392927e-2,0.28525664781722907847e-4,0.54156999310046790024e-7,-0.21589405340123827823e-8,0.22444150951727334619e-10,-0.15368675584000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.10422112945361673560e00,-0.66231638959845581564e-2,0.28800551216363918088e-4,0.37758983397952149613e-7,-0.19435423557038933431e-8,0.20656766125421362458e-10,-0.14422990012444444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.91090275493541084785e-1,-0.65075691516115160062e-2,0.28982078385527224867e-4,0.23014165807643012781e-7,-0.17454532910249875958e-8,0.18981946442680092373e-10,-0.13494234691555555556e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.78191222288771379358e-1,-0.63914190297303976434e-2,0.29079759021299682675e-4,0.97885458059415717014e-8,-0.15635596116134296819e-8,0.17417110744051331974e-10,-0.12591151763555555556e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.65524757106147402224e-1,-0.62750311956082444159e-2,0.29102328354323449795e-4,-0.20430838882727954582e-8,-0.13967781903855367270e-8,0.15958771833747057569e-10,-0.11720175765333333333e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.53091065838453612773e-1,-0.61586898417077043662e-2,0.29057796072960100710e-4,-0.12597414620517987536e-7,-0.12440642607426861943e-8,0.14602787128447932137e-10,-0.10885859114666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, - A9{0.40889797115352738582e-1,-0.60426484889413678200e-2,0.28953496450191694606e-4,-0.21982952021823718400e-7,-0.11044169117553026211e-8,0.13344562332430552171e-10,-0.10091231402844444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00} - }; - - auto w_im_y100 = [](T y100, T x) //purely scalar - { - if (eve::is_nan(x)) return x; - auto t = eve::dec(2*eve::frac(y100)); - int index = ((int) y100); - if (index < 97) return eve::reverse_horner(t, coefs[index]); - else - { - // use Taylor expansion for small x (|x| <= 0.0309...) - // (2/sqrt(pi)) * (x -2/3 x^3 , 4/15 x^5 -8/105 x^7, 16/945 x^9) - auto mx2(-eve::sqr(x)); - kumi::result::generate_t<5, r_t> a{1.1283791670955125739, 0.75225277806367504925, 0.30090111122547001970 - , 0.085971746064420005629, 0.016931216931216931217}; - return x*eve::reverse_horner(mx2, a); - } - }; - - auto taylor = [mz2, z](){ - using r_t = eve::element_type_t; - kumi::result::generate_t<5, r_t> a{1.1283791670955125739, 0.37612638903183752464, 0.11283791670955125739 - , 0.026866170645131251760, 0.0052239776254421878422}; - return eve::reverse_horner(mz2, a)*z; - }; - - if constexpr(eve::scalar_value) - { - auto w_im = [w_im_y100](real_t xx){ - const r_t ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) - auto signxx = sign(xx); - xx = eve::abs(xx); - if (xx > 45) { // continued-fraction expansion is faster - if (xx > real_t(5e7)) // 1-term expansion, important to avoid overflow - return ispi / (xx*signxx); - /* 5-term expansion (rely on compiler for CSE), simplified from: - ispi / (xx-0.5/(xx-1/(xx-1.5/(xx-2/xx)))) */ - return ispi*((xx*xx) * (xx*xx-real_t(4.5)) + 2) / (xx * ((xx*xx) * (xx*xx-5) + real_t(3.75)))*signxx; - } - return w_im_y100(100/(1+xx), xx)*signxx; - }; - - if(eve::is_eqz(y)) - { - return complex(eve::erf(x), y); //call to real implementation - } - else if(eve::is_eqz(x)) - { - return complex(x, (eve::sqr(y) > real_t(720) ? eve::inf(eve::as(y))*eve::sign(y) : eve::expx2(y) * w_im(y))); - } - if (real(mz2) < -750) - { - return complex(eve::signnz(x)); // underflow - } - /* Handle positive and negative x via different formulas, - using the mirror symmetries of w, to avoid overflow/underflow - problems from multiplying exponentially large and small quantities. */ - auto signx = eve::sign(x); - auto ax = eve::abs(x); - auto smallx = ax < 8e-2; - auto smally = eve::abs(y) < 1e-2; - if (smallx && smally) return taylor(); - /* don't use complex exp function, since that will produce spurious NaN - values when multiplying w in an overflow situation. */ - auto [mRe_z2, mIm_z2] = mz2; - auto [s, c] = eve::sincos(mIm_z2); - return oneminus(eve::exp(mRe_z2)*(complex(c, s)*faddeeva(complex(-y, x)*signx)))*signx; - } - else //simd - { - auto signx = eve::signnz(x); - auto w_im = [signx, w_im_y100](real_t x){ - const real_t ispi = real_t(0.56418958354775628694807945156); // 1 / sqrt(pi) - auto greater5e7 = [ispi, signx](auto xx){ - return ispi / (xx*signx); - }; - auto greater45 = [ispi, signx](auto xx){ - auto xx2 = xx*xx; - return signx*ispi*((xx2) * (xx2-real_t(4.5)) + 2) / (xx * ((xx2) * (xx2-5) + real_t(3.75))); - }; - auto remain = [w_im_y100, signx](auto xx){ - - if constexpr(eve::scalar_value) - { - return w_im_y100(100/(1+xx), xx)*signx; - } - else - { - return eve::detail::map(w_im_y100, 100/(1+xx), xx)*signx; - } - }; - - auto ax = eve::abs(x); - auto notdone = eve::true_(eve::as(ax)); - auto r = eve::nan(eve::as(x)); - if( eve::any(notdone) ) - { - notdone = next_interval(greater5e7, notdone, ax > real_t(5e7), r, ax); - if( eve::any(notdone) ) - { - notdone = next_interval(greater45, notdone, ax > real_t(45), r, ax); - if( eve::any(notdone) ) - { - if (eve::any(notdone)) r = remain(ax); - notdone = eve::false_(eve::as(r)); - } - } - } - return r; - }; - - auto signy = eve::signnz(y); - auto no_underflow = real(mz2) >= -750; - auto nfin = eve::is_not_finite(y); //|| is_nan(z); - auto r = if_else(no_underflow || is_nan( real(mz2)) - , complex(eve::nan(eve::as(real_t(0))), eve::nan(eve::as(real_t(0)))) - , complex(signx, real_t(0))); // treat underflow and nan - auto notdone = (no_underflow && !nfin) || eve::is_eqz(y); // no underflow - r = if_else(nfin, complex(eve::nan(eve::as(y)), eve::nan(eve::as(y))), r); - r = if_else(eve::is_eqz(x) && nfin, complex(real_t(0), y), r); - auto ax = eve::abs(x); - auto xsmall = ax < 8e-2; - auto ysmall = eve::abs(y) < 1e-2; - auto nully = [](auto x, auto y){ - return complex(eve::erf(x), y); - }; - auto nullx = [signy, w_im](auto x, auto y){ - return complex(x, eve::if_else(eve::sqr(y) > real_t(720), eve::inf(eve::as(y)), eve::expx2(y) * w_im(y)))*signy; - }; - - auto smallxy = [taylor](){ - return taylor(); - }; - - auto remain = [mz2, signx](auto x, auto y){ - auto [mRe_z2, mIm_z2] = mz2; - auto [s, c] = eve::sincos(mIm_z2); - auto zz = oneminus(eve::exp(mRe_z2)*(complex(c, s)*faddeeva(complex(-y, x)*signx)))*signx; - return zz; - }; - - if( eve::any(notdone) ) - { - notdone = next_interval(nully, notdone, eve::is_eqz(y), r, x, y); - if( eve::any(notdone) ) - { - notdone = next_interval(nullx, notdone, eve::is_eqz(x), r, x, y); - if( eve::any(notdone) ) - { - notdone = next_interval(smallxy, notdone, xsmall && ysmall, r); - if( eve::any(notdone) ) - { - notdone = last_interval(remain, notdone, r, x, y); - } - } - } - } - return r; - } - } -} diff --git a/include/kyosu/types/impl/complex/faddeeva.hpp b/include/kyosu/types/impl/complex/faddeeva.hpp deleted file mode 100644 index 15a0a6b0..00000000 --- a/include/kyosu/types/impl/complex/faddeeva.hpp +++ /dev/null @@ -1,115 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include -#include - -namespace kyosu::_ -{ - template - auto dispatch(eve::tag_of, Z z) noexcept - { - using v_t = as_real_t; - using real_t = eve::element_type_t; - auto const sqrtpi = eve::sqrt_pi(eve::as()); - auto const iosqrtpi = complex(real_t(0), eve::rec(sqrtpi)); - - auto fexp = [iosqrtpi, sqrtpi](auto z){//Fourier expansion approximation - - constexpr real_t tauM = 12; - constexpr real_t tauM2= 144; - constexpr size_t maxN = 23; - constexpr std::array aN = { //Fourier coefficients - 2.758402332921771e-01, 2.245739552246158e-01, 1.594149382739117e-01, 9.866576641545419e-02, 5.324414078763941e-02 - , 2.505215000539365e-02, 1.027746567053954e-02, 3.676164332844847e-03, 1.146493641242233e-03, 3.117570150461976e-04 - , 7.391433429603010e-05, 1.527949342800837e-05, 2.753956608221073e-06, 4.327858781901246e-07, 5.930030408745884e-08 - , 7.084490307748205e-09, 7.379520635816785e-10, 6.702171606002010e-11, 5.307265163470805e-12, 3.664324113467642e-13 - , 2.205894944941035e-14, 1.157826862628556e-15, 5.298711429467307e-17}; - auto z1 = exp_i(tauM*z); - auto z2 = tauM2*sqr(z); - auto FE = sqrtpi/tauM*oneminus(z1)/z2; - for (size_t n = 1; n <= maxN; ++n) - FE += (aN[n-1]*dec(eve::sign_alternate(real_t(n))*z1)/(sqr(n*eve::pi(eve::as(tauM))) - z2)); - return iosqrtpi*tauM2*z*FE; - }; - - auto contfr = [iosqrtpi](auto z){ // the Laplace continued fraction approximation - constexpr std::array bN = {5.5000, 5.0000, 4.5000, 4.0000, 3.5000, - 3.0000, 2.5000, 2.0000, 1.5000, 1.0000, 0.5000}; - auto CF = bN[0]/z; - for (int k = 1; k <= 10; ++k) CF = bN[k]*rec(z - CF); - return iosqrtpi*rec(z - CF); - }; - - auto small_z= [sqrtpi](auto z){ - auto zP2=sqr(z); - auto zP4=sqr(zP2); - auto zP6=zP2*zP4; - - return (((6 - 6*zP2 + 3*zP4 - zP6)*(15*sqrtpi + complex(real_t(0), real_t(1))*z*(30 + 10*zP2 + 3*zP4)))/(90*sqrtpi)); - }; - - auto narr_band = [sqrtpi](auto z, auto aN, auto tauM){ // the narrow band - real_t tauM2 = sqr(tauM); - constexpr size_t maxN = 23; - auto z1 = cos(tauM*z); - auto z2 = tauM2*sqr(z); - Z NB{}; - for (size_t n = 1; n <= maxN; ++n) NB += (aN[n-1]*dec(eve::sign_alternate(real_t(n))*z1)/(sqr(n*eve::pi(eve::as(tauM))) - z2)); - return exp(-sqr(z)) - complex(real_t(0), real_t(1))*(dec(z1)/(tauM*z) - tauM2*z/sqrtpi*NB); - }; - - - auto smallim = [narr_band, small_z](auto z){ // approximation at small imag(z) - constexpr size_t maxN = 23; - constexpr std::array aN12 = { //Fourier coefficients - 2.758402332921771e-01, 2.245739552246158e-01, 1.594149382739117e-01, 9.866576641545419e-02, 5.324414078763941e-02 - , 2.505215000539365e-02, 1.027746567053954e-02, 3.676164332844847e-03, 1.146493641242233e-03, 3.117570150461976e-04 - , 7.391433429603010e-05, 1.527949342800837e-05, 2.753956608221073e-06, 4.327858781901246e-07, 5.930030408745884e-08 - , 7.084490307748205e-09, 7.379520635816785e-10, 6.702171606002010e-11, 5.307265163470805e-12, 3.664324113467642e-13 - , 2.205894944941035e-14, 1.157826862628556e-15, 5.298711429467307e-17}; - constexpr std::array aN12_1 = { //Fourier coefficients avoiding poles - 2.738693653257005e-01, 2.237253191645656e-01, 1.597109174949294e-01, 9.963269094255395e-02, 5.431463971403691e-02 - , 2.587496275171680e-02, 1.077185133049561e-02, 3.918760804860615e-03, 1.245818993506881e-03, 3.461058398397597e-04 - , 8.402542038611651e-05, 1.782626047239582e-05, 3.304894416742214e-06, 5.354300211460283e-07, 7.580461285556943e-08 - , 9.378563805151443e-09, 1.013969351293613e-09, 9.579902872635222e-11, 7.909429715194891e-12, 5.706594634502440e-13 - , 3.597962756907448e-14, 1.982367142666383e-15, 9.544634564831883e-17 }; - - auto x = real(z); - auto ind_0 = abs(x)<5e-3; - auto ind_poles = eve::false_(eve::as(x)); - for(int k = 1; k <= 23; ++k) ind_poles = ind_poles || (abs(x - k*eve::pi(eve::as(x))/12)<1e-4); - return if_else(!ind_0, if_else(ind_poles - , narr_band(z, aN12_1, real_t(12.1)) - , narr_band(z, aN12, real_t(12)) - ) - , small_z(z)); - }; - - auto indneg = eve::is_ltz(imag(z)); - z = if_else(indneg, conj(z), z); - auto r = complex(eve::nan(eve::as(real(z))), eve::nan(eve::as(real(z)))); // nan case treated here - r = if_else(eve::is_infinite(real(z)), Z{}, r); - auto notdone = is_finite(z); - if (eve::any(notdone)) - { - auto indin = sqr_abs(z) <= real_t(64); - auto indband = indin && (imag(z) < real_t(5.0e-3)); - notdone = next_interval(smallim, notdone, indband, r, z); - if( eve::any(notdone) ) - { - notdone = next_interval(fexp, notdone, indin, r, z); - if( eve::any(notdone) ) - { - notdone = next_interval(contfr, notdone, notdone, r, z); - } - } - } - return if_else(indneg, conj(r), r); ; - } -} diff --git a/include/kyosu/types/impl/complex/invtrig.hpp b/include/kyosu/types/impl/complex/invtrig.hpp deleted file mode 100644 index 5cee77b9..00000000 --- a/include/kyosu/types/impl/complex/invtrig.hpp +++ /dev/null @@ -1,583 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include - -namespace kyosu::_ -{ - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - // This implementation is a simd transcription and adaptation of the boost_math code - // which itself is a transcription of the pseudo-code in: - // - // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling." - // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang. - // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997. - // - auto [a0r, a0i] = a0; - using rtype = decltype(a0r); - const rtype a_crossover(1.5); - const rtype b_crossover(0.6417); - auto ltzra0 = eve::is_negative(a0r); - auto ltzia0 = eve::is_negative(a0i); - - // - // Begin by insuring real(a0) >= 0 and imag(a0) >= 0 : - // - - rtype x = eve::abs(a0r); - rtype y = eve::abs(a0i); - rtype proper_real = eve::asin(x); - auto lexone = (x <= eve::one(eve::as(a0r))); - auto is_proper_real = is_real(a0) && lexone; - - auto s_min = eve::sqrtsmallestposval(eve::as(x))*4; - auto s_max = eve::sqrtvalmax(eve::as(x))/8; - rtype xp1 = eve::inc(x); - rtype xm1 = eve::dec(x); - auto not_in_safe_zone = (((x > s_max) || (x < s_min)) || ((y > s_max) || (y < s_min))); - //compute for safe zone - rtype r, i; - rtype yy = eve::sqr(y); - rtype tr = eve::sqrt(eve::sqr(xp1) + yy);//hypot for pedantic ? - rtype ts = eve::sqrt(eve::sqr(xm1) + yy);//hypot for pedantic ? - rtype a = eve::average(tr, ts); - rtype b = x/a; - //compute r for b > b_crossover - rtype apx = a + x; - r = eve::if_else(lexone, - eve::atan(x/eve::sqrt(eve::half(eve::as(a0r))*apx*(yy/(tr+xp1)+(ts-xm1)))), - eve::atan(x/(y*eve::sqrt(eve::half(eve::as(a0r))*(apx/(tr+xp1)+apx/(ts+xm1))))) - ); - // r is computed - r = eve::if_else((b <= b_crossover), r, eve::asin(b)); - //compute am 1 temporary for i for a <= a_crossover - rtype tmp = yy/(tr+xp1); - rtype am1 = eve::if_else(lexone, - eve::average(tmp, yy/(ts-xm1)), - eve::average(tmp, (ts+xm1))); - i = eve::if_else((a <= a_crossover), - eve::log1p(am1 + eve::sqrt(am1*(eve::inc(a)))), - eve::log(a + eve::sqrt(eve::dec(eve::sqr(a)))) - ); - // i is computed - //compute for exception zone - if (eve::any(not_in_safe_zone)) - { - auto zone1 = (y <= eve::eps(eve::as(a0r))*eve::abs(xm1)); - if (eve::any(eve::logical_and(zone1, not_in_safe_zone))) - { - rtype rr = eve::if_else(lexone, eve::asin(x), eve::pio_2(eve::as(a0r))); - rtype ii = eve::if_else(lexone, y/eve::sqrt(xp1*xm1), - eve::if_else((eve::valmax(eve::as(a0r))/xp1 > xm1), - eve::log1p(xm1 + eve::sqrt(xp1*xm1)), - eve::log_2(eve::as(a0r)) + eve::log(x) - ) - ); - r = eve::if_else(zone1, rr, r); - i = eve::if_else(zone1, ii, i); - } - auto zone2 = (y <= s_min); - auto not_treated = eve::logical_notand(zone1, not_in_safe_zone); - if (eve::any(eve::logical_and(zone2, not_treated))) - { - r = eve::if_else(zone2, eve::pio_2(eve::as(a0r)) - eve::sqrt(y), r); - i = eve::if_else(zone2, eve::sqrt(y), i); - } - auto zone3 = (eve::dec(eve::eps(eve::as(a0r))*y) >= x); - not_treated = eve::logical_notand(zone2, not_treated); - if (eve::any(eve::logical_and(zone3, not_treated))) - { - r = eve::if_else(zone3, x/y, r); - i = eve::if_else(zone3, eve::log_2(eve::as(a0r)) + eve::log(y), i); - } - auto zone4 = (x > eve::one(eve::as(a0r))); - not_treated = eve::logical_notand(zone3, not_treated); - if (eve::any(eve::logical_and(zone4, not_treated))) - { - r = eve::if_else(zone4, eve::atan(x/y), r); - i = eve::if_else(zone4, eve::log_2(eve::as(a0r)) + eve::log(y) - + eve::half(eve::as(a0r))*eve::log1p(eve::sqr(x/y)), i); - } - not_treated = eve::logical_notand(zone4, not_treated); - if (eve::any(not_treated)) - { - rtype aa = eve::sqrt(eve::inc(eve::sqr(y))); - r = eve::if_else(not_treated, x/aa, r); - i = eve::if_else(not_treated,eve::half(eve::as(a0r))*eve::log1p(2*y*(y+aa)), i); - } - } - if (eve::any(is_not_finite(a0))) - { - auto nanx = eve::is_nan(x); - auto nany = eve::is_nan(y); - auto infx = (x == eve::inf(eve::as(a0r))) ; - auto infy = (y == eve::inf(eve::as(a0r))) ; - if (eve::any(nanx)) - { - r = eve::if_else(nanx, x, r); - r = eve::if_else(nanx && infy, x, r); - i = eve::if_else(nanx, x, i); - i = eve::if_else(nanx && infy, y, i); - } - if (eve::any(nany)) - { - auto isimag = is_imag(a0); - r = eve::if_else(isimag && nany, eve::zero, r); - r = eve::if_else(eve::logical_and(nany, infx),y, r); - i = eve::if_else(isimag && nany, eve::allbits, i); - i = eve::if_else(nany && infx, x, i); - } - auto test = eve::logical_notand(eve::logical_or(nanx, nany), infx); - if (eve::any(test)) - { - r = eve::if_else(infx && test, - eve::if_else(infy, eve::pio_4(eve::as(a0r)), eve::pio_2(eve::as(a0r))), - r); - } - test = eve::logical_notand(eve::is_nan(x) || eve::is_nan(y), - eve::logical_andnot(infy, infx)); - r = eve::if_else(test,eve::zero,r); - } - // use proper real results - - r = eve::if_else(is_proper_real, proper_real, r); - i = eve::if_else(is_proper_real, eve::zero, i); - // restore signs - r = eve::if_else(ltzra0, -r, r); - i = eve::if_else(ltzia0, -i, i); - return complex(r, i); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - // This implementation is a simd transcription and adaptation of the boost_math code - // which itself is a transcription of the pseudo-code in: - // - // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling." - // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang. - // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997. - // - auto [a0r, a0i] = a0; - using rtype = decltype(a0r); - const rtype a_crossover(1.5); - const rtype b_crossover(0.6417); - auto ltzra0 = eve::is_negative(a0r); - auto gtzia0 = eve::is_positive(a0i); - // - // Begin by insuring a0r >= 0 and imag(a0) >= 0 : - // - rtype x = eve::abs(a0r); - rtype y = eve::abs(a0i); - rtype proper_real = eve::acos(x); - auto lexone = (x <= eve::one(eve::as(x))); - auto is_proper_real = eve::logical_and(is_real(a0), lexone); - - auto s_min = eve::sqrtsmallestposval(eve::as(x))*4; - auto s_max = eve::sqrtvalmax(eve::as(x))/8; - rtype xp1 = eve::inc(x); - rtype xm1 = eve::dec(x); - auto not_in_safe_zone = (((x > s_max) || (x < s_min)) || ((y > s_max) || (y < s_min))); - //compute for safe zone - rtype r, i; - rtype yy = eve::sqr(y); - rtype tr = eve::sqrt(eve::sqr(xp1) + yy); //hypot for pedantic ? - rtype ts = eve::sqrt(eve::sqr(xm1) + yy); //hypot for pedantic ? - rtype a = eve::average(tr, ts); - rtype b = x/a; - //compute r for b > b_crossover - rtype apx = a + x; - r = eve::if_else(lexone, - eve::atan(eve::sqrt(eve::half(eve::as(x))*apx*(yy/(tr+xp1)+(ts-xm1)))/x), - eve::atan((y*eve::sqrt(eve::half(eve::as(x))*(apx/(tr+xp1)+apx/(ts+xm1))))/x) - ); - // r is computed - r = eve::if_else((b <= b_crossover), eve::acos(b), r); - //compute am1 temporary for i for a <= a_crossover - rtype tmp = yy/(tr+xp1); - rtype am1 = eve::if_else(lexone, - eve::average(tmp, yy/(ts-xm1)), - eve::average(tmp, (ts+xm1))); - i = eve::if_else((a <= a_crossover), - eve::log1p(am1 + eve::sqrt(am1*(eve::inc(a)))), - eve::log(a + eve::sqrt(eve::dec(eve::sqr(a)))) - ); - // i is computed - //compute for exception zone - if (eve::any(not_in_safe_zone)) - { - auto zone1 = (y <= eve::eps(eve::as(x))*eve::abs(xm1)); - if (eve::any(eve::logical_and(zone1, not_in_safe_zone))) - { - rtype rr = eve::if_else(lexone, proper_real, eve::zero); - rtype ii = eve::if_else(lexone, y/eve::sqrt(-xp1*xm1), - eve::if_else((eve::valmax(eve::as(x))/xp1 > xm1), - eve::log1p(xm1 + eve::sqrt(xp1*xm1)), - eve::log_2(eve::as(x)) + eve::log(x) - ) - ); - r = eve::if_else(zone1, rr, r); - i = eve::if_else(zone1, ii, i); - } - auto zone2 = (y <= s_min); - auto not_treated = eve::logical_notand(zone1, not_in_safe_zone); - if (eve::any(eve::logical_and(zone2, not_treated))) - { - rtype sqrty = eve::sqrt(y); - r = eve::if_else(zone2, sqrty, r); - i = eve::if_else(zone2, sqrty, i); - } - auto zone3 = (eve::dec(eve::eps(eve::as(x))*y) >= x); - not_treated = eve::logical_notand(zone2, not_treated); - if (eve::any(eve::logical_and(zone3, not_treated))) - { - r = eve::if_else(zone3, eve::pio_2(eve::as(x)), r); - i = eve::if_else(zone3, eve::log_2(eve::as(x)) + eve::log(y), i); - } - auto zone4 = (x > eve::one(eve::as(x))); - not_treated = eve::logical_notand(zone3, not_treated); - if (eve::any(eve::logical_and(zone4, not_treated))) - { - r = eve::if_else(zone4, eve::atan(y/x), r); - i = eve::if_else(zone4, eve::log_2(eve::as(x)) + eve::log(y) + eve::half(eve::as(x))*eve::log1p(eve::sqr(x/y)), i); - } - not_treated = eve::logical_notand(zone4, not_treated); - if (eve::any(not_treated)) - { - rtype aa = eve::sqrt(eve::inc(sqr(y))); - r = eve::if_else(not_treated, eve::pio_2(eve::as(x)), r); - i = eve::if_else(not_treated, eve::half(eve::as(x))*eve::log1p(2*y*(y+aa)), i); - } - } - if (eve::any(kyosu::is_not_finite(a0))) - { - auto nanx = eve::is_nan(x); - auto nany = eve::is_nan(y); - auto infx = (x == eve::inf(eve::as(x))) ; - auto infy = (y == eve::inf(eve::as(x))) ; - if (eve::any(infx)) - { - r = eve::if_else(infx, eve::zero, r); - i = eve::if_else(infx, eve::inf(eve::as(x)), i); - r = eve::if_else(eve::logical_and(infx, infy), eve::pio_4(eve::as(x)), r); - i = eve::if_else(eve::logical_and(infx, infy), eve::inf(eve::as(x)), i); - - r = eve::if_else(eve::logical_and(infx, nany), y, r); - i = eve::if_else(eve::logical_and(infx, nany), eve::minf(eve::as(x)), i); - } - if (eve::any(nanx)) - { - r = eve::if_else(nanx, x, r); - i = eve::if_else(nanx, x, i); - i = eve::if_else(eve::logical_and(nanx, infy), y, i); - } - auto test = eve::logical_notand(eve::logical_or(infx, nanx), infy); - if (eve::any(test)) - { - r = eve::if_else(eve::logical_and(infy, test), eve::pio_2(eve::as(x)), r); - i = eve::if_else(eve::logical_and(infy, test), y, i); - } - test = eve::logical_notand(eve::logical_or(infx, nanx), nany); - r = eve::if_else(test,eve::if_else(is_imag(a0), eve::pio_2(eve::as(x)), y), r); - i = eve::if_else(test,y,i); - } - // use proper real results - r = eve::if_else(is_proper_real, proper_real, r); - i = eve::if_else(is_proper_real, eve::zero, i); - // restore signs - r = eve::if_else(ltzra0, eve::pi(eve::as(x))-r, r); - i = eve::if_else(gtzia0, -i, i); - return complex(r, i); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - // acosh(a0) = +/-i acos(a0) - // Choosing the sign of multiplier to give real(acosh(a0)) >= 0 - // we have compatibility with C99. - auto [r, i] = kyosu::acos(a0); - auto lez = eve::is_negative(i);; - auto res = complex(-i, r); - res = eve::if_else(lez, res, -res); - auto nani = is_nan(i); - if (eve::any(nani)) - return eve::if_else(nani && eve::is_finite(r) - , complex(eve::nan(eve::as(r)), eve::nan(eve::as(r))) - , res); - else - return res; - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - auto [r, i] = a0; - auto [r1, i1] = kyosu::asin(complex(-i, r)); - return complex(i1, -r1); // -(eve::i*asin(eve::i*z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - // This implementation is a simd (i.e. no branch) transcription and adaptation of the - // boost_math code which itself is a transcription of the pseudo-code in: - // - // Eric W. Weisstein. "Inverse Hyperbolic Tangent." - // From MathWorld--A Wolfram Web Resource. - // http://mathworld.wolfram.com/InverseHyperbolicTangent.html - // - // Also: The Wolfram Functions Site, - // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/ - // - // Also "Abramowitz and Stegun. Handbook of Mathematical Functions." - // at : http://jove.prohosting.com/~skripty/toc.htm - // - auto [a0r, a0i] = a0; - auto realinf = eve::is_eqz(a0i) && eve::is_infinite(a0r); - using rtype = decltype(a0r); - const rtype alpha_crossover(0.3); - auto ltzra0 = eve::is_ltz(a0r); - auto ltzia0 = eve::is_ltz(a0i); - auto s_min = eve::sqrtsmallestposval(eve::as(a0r))*2; - auto s_max = eve::sqrtvalmax(eve::as(a0r))/2; - rtype const two = rtype(2); - rtype inf = eve::inf(eve::as(a0r)); - rtype x = eve::abs(a0r); - rtype y = eve::abs(a0i); - auto special = eve::is_eqz(y) && (x < eve::one(eve::as(a0r))); - auto sr = eve::atanh(a0r); - if (eve::all(special)) { - return complex(sr, eve::zero(eve::as(sr))); - } - - rtype r = eve::zero(eve::as(a0r)); - rtype i = eve::zero(eve::as(a0r)); - auto gtxmax = (x > s_max); - auto ltxmin = (x < s_min); - auto gtymax = (y > s_max); - auto ltymin = (y < s_min); - rtype xx = eve::sqr(x); - rtype yy = eve::sqr(y); - rtype sqrabs = xx + yy; - - auto not_in_safe_zone = ((gtxmax || ltxmin) || (gtymax || ltymin)); - if(eve::any(not_in_safe_zone)) - { - //treat underflow or overflow - // one or both of x and y are small, calculate divisor carefully: - rtype div = eve::one(eve::as(a0r)); - div += eve::if_else(ltxmin, xx, eve::zero); - div += eve::if_else(ltxmin, yy, eve::zero); - - rtype alpha = x/div; - alpha += alpha; - - auto test = gtymax; - // big y, medium x, divide through by y: - rtype tmp_alpha = (two*x/y) / (y + xx/y); - // small x and y, whatever alpha is, it's too small to calculate: - tmp_alpha = eve::if_else(x > eve::one(eve::as(a0r)), tmp_alpha, eve::zero); - alpha = eve::if_else(test && (x > eve::one(eve::as(a0r))), tmp_alpha, alpha); - - test = eve::logical_andnot(gtxmax, test); - - // big x small y, as above but neglect y^2/x: - tmp_alpha = two/x; - // big x: divide through by x: - tmp_alpha = eve::if_else((y > eve::one(eve::as(a0r))), two / (x + y*y/x), tmp_alpha); - // big x and y: divide alpha through by x*y: - tmp_alpha = eve::if_else(gtymax, (two/y) / (x/y + y/x), tmp_alpha); - // x or y are infinite: the result is 0 - tmp_alpha = eve::if_else((y == inf) || (x == inf), eve::zero, tmp_alpha); - - alpha = eve::if_else(test, tmp_alpha, alpha); - r = eve::if_else((alpha < alpha_crossover), - eve::log1p(alpha) - eve::log1p(-alpha), - eve::log(inc(two*x + xx)) - eve::log(sqr(dec(x))) - ); - test = (x == eve::one(eve::as(a0r))) && ltymin; - r = eve::if_else(test, -(two*(eve::log(y) - eve::log_2(eve::as(a0r)))), r); - r *= rtype(0.25); - //compute the imag part - // y^2 is negligible: - i = eve::atan2(two*y, eve::oneminus(xx)); - i = eve::if_else(gtymax || gtxmax, eve::pi(eve::as(a0r)), i); - rtype tmp_i = eve::if_else(ltymin, eve::atan2(two*y, eve::one(eve::as(a0r))), - eve::atan2(two*y, eve::oneminus(yy))); - i = eve::if_else(ltxmin, tmp_i, i); - } - auto test0 = (inf == x) && (inf == y); - if(eve::any(test0)) - { - //inf x, inf y - r = eve::if_else(test0, eve::zero, r); - i = eve::if_else(test0, eve::pi(eve::as(a0r)), r); - } - auto test = kyosu::is_nan(a0); - - if(eve::any(test)) - { - //nan x, inf y - r = eve::if_else(eve::is_nan(x) && (y == inf), eve::zero, r); - i = eve::if_else(eve::is_nan(x) && (y == inf), eve::pi(eve::as(a0r)), r); - - r = eve::if_else(is_nan(y) && (x == inf), eve::zero, r); - i = eve::if_else(is_nan(y) && (x == inf), y, i); - - r = eve::if_else(is_nan(y) && eve::is_eqz(x), eve::zero, r); - i = eve::if_else(is_nan(y) && is_eqz(x), eve::allbits, i); - } - //compute for safe zone::one - // the real part is given by: - // - // eve::real(atanh(z)) == eve::log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x)) - // - // however, when x is either large (x > 1/e) or very small - // (x < e) then this effectively simplifies - // to log(1), leading to wildly inaccurate results. - // by dividing the above (top and bottom) by (1 + x^2 + y^2) we get: - // - // eve::real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2)))) - // - // which is much more sensitive to the value of x, when x is not near 1 - // (remember we can compute log(1+x) for small x very accurately). - // - // the cross-over from eve::one method to the other has to be determined - // experimentally, the value used below appears correct to within a - // factor of 2 (and there are larger errors from other parts - // of the input domain anyway). - // - rtype alpha = x*two / (eve::inc(sqrabs)); - rtype sqrxm1 = eve::sqr(eve::dec(x)); - rtype tmp_r = eve::if_else((alpha < alpha_crossover), - eve::log1p(alpha) - eve::log1p(-alpha), - eve::log1p(x+x + sqrabs) - eve::log(sqrxm1 + yy) - )*rtype(0.25); - r = eve::if_else(not_in_safe_zone, r, tmp_r); - - // compute the imag part - i = eve::if_else(not_in_safe_zone, - i, - eve::atan2(y+y, (eve::oneminus(sqrabs))) - )*eve::half(eve::as(a0r)); - - r = eve::if_else( ltzra0,-r, r); - i = eve::if_else(eve::is_infinite(y), eve::pio_2(eve::as(a0r))*eve::sign(y), i); - i = eve::if_else( ltzia0,-i, i); - r = eve::if_else(realinf, eve::zero(eve::as(a0r)), r); - i = eve::if_else(realinf, -eve::sign(a0r)*eve::pio_2(eve::as(a0r)), i); - r = eve::if_else(special, sr, r); - i = eve::if_else(special, eve::zero, i); - return complex(r, i); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - - // C99 definition here; atan(z) = -i atanh(iz): - auto [r, i] = a0; - auto [r1, i1] = kyosu::atanh(complex(-i, r)); - return complex(i1, -r1); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return kyosu::acos(kyosu::rec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return kyosu::asin(kyosu::rec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return kyosu::acosh(kyosu::rec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return kyosu::asinh(kyosu::rec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::atan(kyosu::rec(a0)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::atanh(kyosu::rec(a0)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::radinpi(kyosu::acos(a0)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::radinpi(kyosu::acot(a0)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::radinpi(kyosu::acsc(a0)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::radinpi(kyosu::asec(a0)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::radinpi(kyosu::asin(a0)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& a0) noexcept - { - return kyosu::radinpi(kyosu::atan(a0)); - } - -} diff --git a/include/kyosu/types/impl/complex/predicates.hpp b/include/kyosu/types/impl/complex/predicates.hpp deleted file mode 100644 index 9398fbcc..00000000 --- a/include/kyosu/types/impl/complex/predicates.hpp +++ /dev/null @@ -1,21 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include - -namespace kyosu::_ -{ - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& c) noexcept - { - return eve::is_eqz(kyosu::real(c)); - } - -} diff --git a/include/kyosu/types/impl/complex/special.hpp b/include/kyosu/types/impl/complex/special.hpp deleted file mode 100644 index ab51c7cc..00000000 --- a/include/kyosu/types/impl/complex/special.hpp +++ /dev/null @@ -1,423 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include -#include - -namespace kyosu::_ -{ - //===------------------------------------------------------------------------------------------- - // Unary functions : tgamma - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, Z const& a0) noexcept - { - // 15 sig. digits for 0<=real(z)<=171 - // coeffs should sum to about g*g/2+23/24 - // - using r_t = eve::element_type_t>; - auto g=r_t(607)/r_t(128); - // best results when 4<=g<=5 - constexpr int N = 15; - std::array c = - { 0.99999999999999709182, - 57.156235665862923517, - -59.597960355475491248, - 14.136097974741747174, - -0.49191381609762019978, - .33994649984811888699e-4, - .46523628927048575665e-4, - -.98374475304879564677e-4, - .15808870322491248884e-3, - -.21026444172410488319e-3, - .21743961811521264320e-3, - -.16431810653676389022e-3, - .84418223983852743293e-4, - -.26190838401581408670e-4, - .36899182659531622704e-5 - }; - - //Num Recipes used g=5 with 7 terms - //for a less effective approximation - - auto negra0 = eve::is_negative(real(a0)); - auto z = if_else(negra0, -a0, a0); - z = dec(z); - auto zh = z+eve::half(eve::as()); - auto zgh=zh+g; - //trick for avoiding FP overflow above z=141 - auto zp=pow(zgh,(zh*eve::half(eve::as()))); - auto ss = Z{}; - for(int pp = N-1; pp >= 1; --pp){ - ss+= c[pp]/(z+pp); - } - auto sq2pi = r_t(2.5066282746310005024157652848110); - auto f=(sq2pi*(c[0]+ss))*((zp*exp(-zgh))*zp); - auto o = eve::one(eve::as()); - f = if_else(is_eqz(z) || z == o, o, f); - //adjust for negative real parts - auto reala0 = is_real(a0); - if(eve::any(negra0)) - { - f = if_else(negra0, rec(-eve::inv_pi(eve::as(real(a0)))*a0*f*sinpi(a0)), eve::zero); - f = if_else (negra0 && reala0 && eve::is_flint(real(a0)), complex(eve::nan(eve::as(sq2pi)), eve::inf(eve::as(sq2pi))), f); - } - f = if_else (reala0, complex(eve::tgamma(real(a0))), f); - f = if_else (eve::is_nan(real(f)), complex(eve::nan(eve::as(sq2pi)), eve::inf(eve::as(sq2pi))), f); - f = if_else (is_eqz(a0), complex(eve::inf(eve::as(g))*eve::pedantic(eve::signnz)(real(a0))), f); - return f; - } - - //===------------------------------------------------------------------------------------------- - // Unary functions : log_gamma - //===------------------------------------------------------------------------------------------- - template - EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& a0) noexcept - { - return log_abs(tgamma(a0)); - } - - template - auto dispatch(eve::tag_of, Z const& a0) noexcept - { - // 15 sig. digits for 0<=real(z)<=171 - // coeffs should sum to about g*g/2+23/24 - // - using r_t = eve::element_type_t>; - auto g=r_t(607)/r_t(128); - // best results when 4<=g<=5 - constexpr int N = 15; - std::array c = - { 0.99999999999999709182, - 57.156235665862923517, - -59.597960355475491248, - 14.136097974741747174, - -0.49191381609762019978, - .33994649984811888699e-4, - .46523628927048575665e-4, - -.98374475304879564677e-4, - .15808870322491248884e-3, - -.21026444172410488319e-3, - .21743961811521264320e-3, - -.16431810653676389022e-3, - .84418223983852743293e-4, - -.26190838401581408670e-4, - .36899182659531622704e-5 - }; - - //Num Recipes used g=5 with 7 terms - //for a less effective approximation - - auto negra0 = eve::is_negative(real(a0)); - auto z = if_else(negra0, -a0, a0); - Z ss{}; - for(int pp = N-1; pp >= 1; --pp){ - ss += c[pp]*rec(z+eve::dec(pp)); - } - auto zg = z+g-eve::half(eve::as(g)); - auto lsq2pi = r_t(0.9189385332046727417803297); - auto f=(lsq2pi + log(c[0]+ss)) - zg + (z-eve::half(eve::as()))*log(zg); - auto zer = eve::zero(eve::as(g)); - auto o = eve::one(eve::as(g)); - auto t = o+o; - f = if_else(z == t|| z == o, zer, f); - //adjust for negative real parts - if(eve::any(negra0)) - { - auto lpi = r_t(1.14472988584940017414342735); - auto reala0 = is_real(a0); - f = kyosu::if_else(negra0, lpi-log(a0*sinpi(-a0))-f, f); - f = kyosu::if_else (negra0 && reala0 && eve::is_flint(real(a0)) - , complex(eve::nan(eve::as(g)), eve::inf(eve::as(g))) - , f); - } - return f; - } - - //===------------------------------------------------------------------------------------------- - // Unary functions : digamma - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, Z const& a0) noexcept - { - // 15 sig. digits for 0<=real(z)<=171 - // coeffs should sum to about g*g/2+23/24 - // - using r_t = eve::element_type_t>; - auto g=r_t(607)/r_t(128); - // best results when 4<=g<=5 - constexpr int N = 15; - std::array c = - { 0.99999999999999709182, - 57.156235665862923517, - -59.597960355475491248, - 14.136097974741747174, - -0.49191381609762019978, - .33994649984811888699e-4, - .46523628927048575665e-4, - -.98374475304879564677e-4, - .15808870322491248884e-3, - -.21026444172410488319e-3, - .21743961811521264320e-3, - -.16431810653676389022e-3, - .84418223983852743293e-4, - -.26190838401581408670e-4, - .36899182659531622704e-5 - }; - - //Num Recipes used g=5 with 7 terms - //for a less effective approximation - - auto reflection = real(a0) < eve::half(eve::as(real(a0))); - auto z = if_else(reflection, oneminus(a0), a0); - - Z d{}; - auto n = d; - for(int pp = N-1; pp >= 1; --pp){ - auto dz = rec(z+eve::dec(pp)); - auto dd = c[pp]*dz; - d += dd; - n -= dd*dz; - } - d+= c[0]; - auto zg = z+g-eve::half(eve::as(g)); - auto f = log(zg) + (n/d - g/zg); - - if(eve::any(reflection)) - { - f = if_else(reflection, f-eve::pi(eve::as(g))*cotpi(a0), f); - f = if_else (reflection && is_real(a0) && eve::is_flint(real(a0)) - , complex(eve::nan(eve::as(g)), eve::inf(eve::as(g))), f); - } - return f; - } - - //===------------------------------------------------------------------------------------------- - // Binary functions : rising_factorial, lrising_factorial, lbeta, beta - //===------------------------------------------------------------------------------------------- - // WAITING FOR DECORATORS - template - auto dispatch(eve::tag_of - , Z1 const& a0, Z2 const& a1) noexcept - requires (kyosu::concepts::complex || kyosu::concepts::complex) - { - using r_t = as_cayley_dickson_t; - return kyosu::if_else(is_eqz(a1), r_t(1), tgamma(a0+a1)/tgamma(a0)); //pedantic(div)(tgamma(a0+a1),tgamma(a0))); - } - - template - auto dispatch(eve::tag_of - , Z1 const& a0, Z2 const& a1) noexcept - requires (kyosu::concepts::complex || kyosu::concepts::complex) - { - using r_t = as_cayley_dickson_t; - return kyosu::if_else(is_eqz(a1), r_t{}, log(tgamma(a0+a1)/tgamma(a0))); //pedantic(divide)(tgamma(a0+a1),tgamma(a0)))); - } - - template - auto dispatch(eve::tag_of, Z1 const& a0, Z2 const& a1) noexcept - requires (kyosu::concepts::complex || kyosu::concepts::complex) - { - return log(beta(a0, a1)); - } - - template - auto dispatch(eve::tag_of, Z1 const& a0, Z2 const& a1) noexcept - requires (kyosu::concepts::complex || kyosu::concepts::complex) - { - auto y = a0 + a1; - return tgamma(a0)*tgamma(a1)/tgamma(y); - } - - template - auto dispatch(eve::tag_of - , K const & kk - , Z z) noexcept - { - // 15 sig. digits for 0<=real(z)<=171 - // coeffs should sum to about g*g/2+23/24 - // - using v_t = kyosu::as_real_t; - using real_t = eve::element_type_t; - auto k = real_t(kk); - auto [rz, iz] = z; - auto iszero = is_eqz(z); - auto isreal = eve::is_eqz(iz); - auto isneg = eve::is_ltz(rz); - auto isnegreal = isreal && isneg; - auto isodd = isnegreal && eve::is_odd(rz); - auto iseven = isnegreal && eve::is_even(rz); - - auto r=eve::half(eve::as(rz)); // reflection point - auto reflect = rz < r; - z = if_else(reflect, oneminus(z), z); - - std::array cm = { - .27105054312137610850e-19 - -.17889335846010823161e-17, - .58167446553847312884e-16, - -.12421162189080181548e-14, - .19593322190397666205e-13, - -.24347803504257137241e-12, - .24823251635643084345e-11, - -.21352608103961806529e-10, - .15816215942184366772e-9, - -.10246226511017621219e-8, - .58768014045093054654e-8, - -.30137695171547022183e-7, - .13931171712321674741e-6, - -.58440580661848562719e-6, - .22376124247437700378e-5, - -.78585149263697370338e-5, - .25423835243950883896e-4, - -.76053287924037718971e-4, - .21106516173760261250e-3, - -.54504190222378945440e-3, - .13131884053420191908e-2, - -.29592166263096543401e-2, - .62512730682449246388e-2, - -.12405987285776082154e-1, - .23176737166455607805e-1, - -.40840766970770029873e-1, - .68016197438946063823e-1, - -.10726959700408922397, - .16054206784249779846, - -.22851039270529494523, - .31007238254065152134, - -.40215850009669926857, - .50000000000000000000, - -.59784149990330073143, - .68992761745934847866, - -.77148960729470505477, - .83945793215750220154, - -.89273040299591077603, - .93198380256105393618, - -.95915923302922997013, - .97682326283354439220, - -.98759401271422391785, - .99374872693175507536, - -.99704078337369034566, - .99868681159465798081, - -.99945495809777621055, - .99978893483826239739, - -.99992394671207596228, - .99997457616475604912, - -.99999214148507363026, - .99999776238757525623, - -.99999941559419338151, - .99999986068828287678, - -.99999996986230482845, - .99999999412319859549, - -.99999999897537734890, - .99999999984183784058, - -.99999999997864739190, - .99999999999751767484, - -.99999999999975652196, - .99999999999998040668, - -.99999999999999875788, - .99999999999999994183, - -.99999999999999999821, - .99999999999999999997 - }; - - Z f{}; - auto j = inc(63*k); - for(size_t i=0; i < cm.size(); ++i, j -= k) - { - f += cm[i]*pow(j, -z); - } - if (eve::none(reflect)) return f; - auto reflection = [&](auto f){ - if (k == 1) - { - auto u = dec(exp2(z)); - auto t = -2*u/dec(u)/pow(eve::pi(eve::as(real(z))), z);; - auto g = t*cospi(eve::half(eve::as(real(z)))*z)*tgamma(z)*f; - g = if_else(iszero, complex(eve::half(eve::as(real(f)))), g); - return if_else(iseven, complex(eve::zero(eve::as(real(f)))), g); - } - else if (k == 2){ - auto g= pow(eve::two_o_pi(eve::as(real(f))), z)*sinpi(z)*tgamma(z)*f; - return if_else(isodd, eve::zero(eve::as(real(f))), g); - } - else - { - return complex(eve::allbits(eve::as(real(f)))); - } - }; - return if_else(reflect, reflection(f), f); - } - - //===------------------------------------------------------------------------------------------- - // Unary functions : eta - //===------------------------------------------------------------------------------------------- - template - EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept - { - return deta(1u, z); - } - - //===------------------------------------------------------------------------------------------- - // Unary functions : zeta - //===------------------------------------------------------------------------------------------- - template - EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept - { - auto zz=exp2(z); - auto k = zz/(zz-2); - auto g = if_else(z == Z(1), complex(eve::nan(eve::as(real(z)))), k*eta(z)); - return if_else(real(z) == eve::inf(eve::as(real(z))), complex(eve::one(eve::as(real(z)))), g); - } - - //===------------------------------------------------------------------------------------------- - // Unary functions : lambda - //===------------------------------------------------------------------------------------------- - template - EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept - { - auto zz=exp2(z); - auto k = (z-1)/(z-2); - auto r = if_else(z == complex(eve::one(eve::as(real(z)))), complex(eve::inf(eve::as(real(z)))), k*deta(1u, zz)); - imag(r) = eve::if_else(is_real(z), eve::zero, imag(r)); - return r; - } -} - -#include -#include - -namespace kyosu::_ -{ - template - EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept - { - auto realz = is_real(z); - if (eve::all(realz)) - return complex(eve::erfcx(real(z))); - else if (eve::none(realz)) - return faddeeva(complex(-imag(z), real(z))); - else - return if_else(realz, complex(eve::erfcx(real(z))), faddeeva(complex(-imag(z), real(z)))); - } - - template - EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept - { - auto realz = is_real(z); - if (eve::all(realz)) - { - return kyosu::erfi(real(z)); - } - else - { - auto [rz, iz] = z; - auto tmp = erf(complex(-iz, rz)); - return complex(imag(tmp), -real(tmp)); - } - } -} diff --git a/include/kyosu/types/impl/erf.hpp b/include/kyosu/types/impl/erf.hpp new file mode 100644 index 00000000..d3b755db --- /dev/null +++ b/include/kyosu/types/impl/erf.hpp @@ -0,0 +1,298 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include +#include +#include + + +namespace kyosu::_ +{ + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + auto numeric_sqr = [](auto z){ // UNTIL DECORATORS ARE AT HAND + auto [zr, zi] = z; + return complex((zr-zi)*(zi+zr), 2 * zr * zi); + }; + auto x = real(z); + auto y = imag(z); + using real_t = decltype(x); + using r_t = eve::element_type_t; + auto mz2 = -numeric_sqr(z);// -z^2, being careful of overflow + using A9 = kumi::result::generate_t<9, r_t>; + static constexpr std::array< A9, 97> coefs = + { + A9{0.28351593328822191546e-2,0.28494783221378400759e-2,0.14427470563276734183e-4,0.10939723080231588129e-6,0.92474307943275042045e-9,0.89128907666450075245e-11,0.92974121935111111110e-13,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.85927161243940350562e-2,0.29085312941641339862e-2,0.15106783707725582090e-4,0.11716709978531327367e-6,0.10197387816021040024e-8,0.10122678863073360769e-10,0.10917479678400000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.14471159831187703054e-1,0.29703978970263836210e-2,0.15835096760173030976e-4,0.12574803383199211596e-6,0.11278672159518415848e-8,0.11547462300333495797e-10,0.12894535335111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.20476320420324610618e-1,0.30352843012898665856e-2,0.16617609387003727409e-4,0.13525429711163116103e-6,0.12515095552507169013e-8,0.13235687543603382345e-10,0.15326595042666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.26614461952489004566e-1,0.31034189276234947088e-2,0.17460268109986214274e-4,0.14582130824485709573e-6,0.13935959083809746345e-8,0.15249438072998932900e-10,0.18344741882133333333e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.32892330248093586215e-1,0.31750557067975068584e-2,0.18369907582308672632e-4,0.15761063702089457882e-6,0.15577638230480894382e-8,0.17663868462699097951e-10,0.22126732680711111111e-12,0.30273474177737853668e-14,0.00000000000000000000e+00}, + A9{0.39317207681134336024e-1,0.32504779701937539333e-2,0.19354426046513400534e-4,0.17081646971321290539e-6,0.17485733959327106250e-8,0.20593687304921961410e-10,0.26917401949155555556e-12,0.38562123837725712270e-14,0.00000000000000000000e+00}, + A9{0.45896976511367738235e-1,0.33300031273110976165e-2,0.20423005398039037313e-4,0.18567412470376467303e-6,0.19718038363586588213e-8,0.24175006536781219807e-10,0.33059982791466666666e-12,0.49756574284439426165e-14,0.00000000000000000000e+00}, + A9{0.52640192524848962855e-1,0.34139883358846720806e-2,0.21586390240603337337e-4,0.20247136501568904646e-6,0.22348696948197102935e-8,0.28597516301950162548e-10,0.41045502119111111110e-12,0.65151614515238361946e-14,0.00000000000000000000e+00}, + A9{0.59556171228656770456e-1,0.35028374386648914444e-2,0.22857246150998562824e-4,0.22156372146525190679e-6,0.25474171590893813583e-8,0.34122390890697400584e-10,0.51593189879111111110e-12,0.86775076853908006938e-14,0.00000000000000000000e+00}, + A9{0.66655089485108212551e-1,0.35970095381271285568e-2,0.24250626164318672928e-4,0.24339561521785040536e-6,0.29221990406518411415e-8,0.41117013527967776467e-10,0.65786450716444444445e-12,0.11791885745450623331e-13,0.00000000000000000000e+00}, + A9{0.73948106345519174661e-1,0.36970297216569341748e-2,0.25784588137312868792e-4,0.26853012002366752770e-6,0.33763958861206729592e-8,0.50111549981376976397e-10,0.85313857496888888890e-12,0.16417079927706899860e-13,0.00000000000000000000e+00}, + A9{0.81447508065002963203e-1,0.38035026606492705117e-2,0.27481027572231851896e-4,0.29769200731832331364e-6,0.39336816287457655076e-8,0.61895471132038157624e-10,0.11292303213511111111e-11,0.23558532213703884304e-13,0.00000000000000000000e+00}, + A9{0.89166884027582716628e-1,0.39171301322438946014e-2,0.29366827260422311668e-4,0.33183204390350724895e-6,0.46276006281647330524e-8,0.77692631378169813324e-10,0.15335153258844444444e-11,0.35183103415916026911e-13,0.00000000000000000000e+00}, + A9{0.97121342888032322019e-1,0.40387340353207909514e-2,0.31475490395950776930e-4,0.37222714227125135042e-6,0.55074373178613809996e-8,0.99509175283990337944e-10,0.21552645758222222222e-11,0.55728651431872687605e-13,0.00000000000000000000e+00}, + A9{0.10532778218603311137e00,0.41692873614065380607e-2,0.33849549774889456984e-4,0.42064596193692630143e-6,0.66494579697622432987e-8,0.13094103581931802337e-09,0.31896187409777777778e-11,0.97271974184476560742e-13,0.00000000000000000000e+00}, + A9{0.11380523107427108222e00,0.43099572287871821013e-2,0.36544324341565929930e-4,0.47965044028581857764e-6,0.81819034238463698796e-8,0.17934133239549647357e-09,0.50956666166186293627e-11,0.18850487318190638010e-12,0.79697813173519853340e-14}, + A9{0.12257529703447467345e00,0.44621675710026986366e-2,0.39634304721292440285e-4,0.55321553769873381819e-6,0.10343619428848520870e-7,0.26033830170470368088e-09,0.87743837749108025357e-11,0.34427092430230063401e-12,0.10205506615709843189e-13}, + A9{0.13166276955656699478e00,0.46276970481783001803e-2,0.43225026380496399310e-4,0.64799164020016902656e-6,0.13580082794704641782e-7,0.39839800853954313927e-09,0.14431142411840000000e-10,0.42193457308830027541e-12,0.00000000000000000000e+00}, + A9{0.14109647869803356475e00,0.48088424418545347758e-2,0.47474504753352150205e-4,0.77509866468724360352e-6,0.18536851570794291724e-7,0.60146623257887570439e-09,0.18533978397305276318e-10,0.41033845938901048380e-13,-0.46160680279304825485e-13}, + A9{0.15091057940548936603e00,0.50086864672004685703e-2,0.52622482832192230762e-4,0.95034664722040355212e-6,0.25614261331144718769e-7,0.80183196716888606252e-09,0.12282524750534352272e-10,-0.10531774117332273617e-11,-0.86157181395039646412e-13}, + A9{0.16114648116017010770e00,0.52314661581655369795e-2,0.59005534545908331315e-4,0.11885518333915387760e-5,0.33975801443239949256e-7,0.82111547144080388610e-09,-0.12357674017312854138e-10,-0.24355112256914479176e-11,-0.75155506863572930844e-13}, + A9{0.17185551279680451144e00,0.54829002967599420860e-2,0.67013226658738082118e-4,0.14897400671425088807e-5,0.40690283917126153701e-7,0.44060872913473778318e-09,-0.52641873433280000000e-10,-0.30940587864543343124e-11,0.00000000000000000000e+00}, + A9{0.18310194559815257381e00,0.57701559375966953174e-2,0.76948789401735193483e-4,0.18227569842290822512e-5,0.41092208344387212276e-7,-0.44009499965694442143e-09,-0.92195414685628803451e-10,-0.22657389705721753299e-11,0.10004784908106839254e-12}, + A9{0.19496527191546630345e00,0.61010853144364724856e-2,0.88812881056342004864e-4,0.21180686746360261031e-5,0.30652145555130049203e-7,-0.16841328574105890409e-08,-0.11008129460612823934e-09,-0.12180794204544515779e-12,0.15703325634590334097e-12}, + A9{0.20754006813966575720e00,0.64825787724922073908e-2,0.10209599627522311893e-3,0.22785233392557600468e-5,0.73495224449907568402e-8,-0.29442705974150112783e-08,-0.94082603434315016546e-10,0.23609990400179321267e-11,0.14141908654269023788e-12}, + A9{0.22093185554845172146e00,0.69182878150187964499e-2,0.11568723331156335712e-3,0.22060577946323627739e-5,-0.26929730679360840096e-7,-0.38176506152362058013e-08,-0.47399503861054459243e-10,0.40953700187172127264e-11,0.69157730376118511127e-13}, + A9{0.23524827304057813918e00,0.74063350762008734520e-2,0.12796333874615790348e-3,0.18327267316171054273e-5,-0.66742910737957100098e-7,-0.40204740975496797870e-08,0.14515984139495745330e-10,0.44921608954536047975e-11,-0.18583341338983776219e-13}, + A9{0.25058626331812744775e00,0.79377285151602061328e-2,0.13704268650417478346e-3,0.11427511739544695861e-5,-0.10485442447768377485e-6,-0.34850364756499369763e-08,0.72656453829502179208e-10,0.36195460197779299406e-11,-0.84882136022200714710e-13}, + A9{0.26701724900280689785e00,0.84959936119625864274e-2,0.14112359443938883232e-3,0.17800427288596909634e-6,-0.13443492107643109071e-6,-0.23512456315677680293e-08,0.11245846264695936769e-09,0.19850501334649565404e-11,-0.11284666134635050832e-12}, + A9{0.28457293586253654144e00,0.90581563892650431899e-2,0.13880520331140646738e-3,-0.97262302362522896157e-6,-0.15077100040254187366e-6,-0.88574317464577116689e-09,0.12760311125637474581e-09,0.20155151018282695055e-12,-0.10514169375181734921e-12}, + A9{0.30323425595617385705e00,0.95968346790597422934e-2,0.12931067776725883939e-3,-0.21938741702795543986e-5,-0.15202888584907373963e-6,0.61788350541116331411e-09,0.11957835742791248256e-09,-0.12598179834007710908e-11,-0.75151817129574614194e-13}, + A9{0.32292521181517384379e00,0.10082957727001199408e-1,0.11257589426154962226e-3,-0.33670890319327881129e-5,-0.13910529040004008158e-6,0.19170714373047512945e-08,0.94840222377720494290e-10,-0.21650018351795353201e-11,-0.37875211678024922689e-13}, + A9{0.34351233557911753862e00,0.10488575435572745309e-1,0.89209444197248726614e-4,-0.43893459576483345364e-5,-0.11488595830450424419e-6,0.28599494117122464806e-08,0.61537542799857777779e-10,-0.24935749227658002212e-11,0.00000000000000000000e+00}, + A9{0.36480946642143669093e00,0.10789304203431861366e-1,0.60357993745283076834e-4,-0.51855862174130669389e-5,-0.83291664087289801313e-7,0.33898011178582671546e-08,0.27082948188277716482e-10,-0.23603379397408694974e-11,0.19328087692252869842e-13}, + A9{0.38658679935694939199e00,0.10966119158288804999e-1,0.27521612041849561426e-4,-0.57132774537670953638e-5,-0.48404772799207914899e-7,0.35268354132474570493e-08,-0.32383477652514618094e-11,-0.19334202915190442501e-11,0.32333189861286460270e-13}, + A9{0.40858275583808707870e00,0.11006378016848466550e-1,-0.76396376685213286033e-5,-0.59609835484245791439e-5,-0.13834610033859313213e-7,0.33406952974861448790e-08,-0.26474915974296612559e-10,-0.13750229270354351983e-11,0.36169366979417390637e-13}, + A9{0.43051714914006682977e00,0.10904106549500816155e-1,-0.43477527256787216909e-4,-0.59429739547798343948e-5,0.17639200194091885949e-7,0.29235991689639918688e-08,-0.41718791216277812879e-10,-0.81023337739508049606e-12,0.33618915934461994428e-13}, + A9{0.45210428135559607406e00,0.10659670756384400554e-1,-0.78488639913256978087e-4,-0.56919860886214735936e-5,0.44181850467477733407e-7,0.23694306174312688151e-08,-0.49492621596685443247e-10,-0.31827275712126287222e-12,0.27494438742721623654e-13}, + A9{0.47306491195005224077e00,0.10279006119745977570e-1,-0.11140268171830478306e-3,-0.52518035247451432069e-5,0.64846898158889479518e-7,0.17603624837787337662e-08,-0.51129481592926104316e-10,0.62674584974141049511e-13,0.20055478560829935356e-13}, + A9{0.49313638965719857647e00,0.97725799114772017662e-2,-0.14122854267291533334e-3,-0.46707252568834951907e-5,0.79421347979319449524e-7,0.11603027184324708643e-08,-0.48269605844397175946e-10,0.32477251431748571219e-12,0.12831052634143527985e-13}, + A9{0.51208057433416004042e00,0.91542422354009224951e-2,-0.16726530230228647275e-3,-0.39964621752527649409e-5,0.88232252903213171454e-7,0.61343113364949928501e-09,-0.42516755603130443051e-10,0.47910437172240209262e-12,0.66784341874437478953e-14}, + A9{0.52968945458607484524e00,0.84400880445116786088e-2,-0.18908729783854258774e-3,-0.32725905467782951931e-5,0.91956190588652090659e-7,0.14593989152420122909e-09,-0.35239490687644444445e-10,0.54613829888448694898e-12,0.00000000000000000000e+00}, + A9{0.54578857454330070965e00,0.76474155195880295311e-2,-0.20651230590808213884e-3,-0.25364339140543131706e-5,0.91455367999510681979e-7,-0.23061359005297528898e-09,-0.27512928625244444444e-10,0.54895806008493285579e-12,0.00000000000000000000e+00}, + A9{0.56023851910298493910e00,0.67938321739997196804e-2,-0.21956066613331411760e-3,-0.18181127670443266395e-5,0.87650335075416845987e-7,-0.51548062050366615977e-09,-0.20068462174044444444e-10,0.50912654909758187264e-12,0.00000000000000000000e+00}, + A9{0.57293478057455721150e00,0.58965321010394044087e-2,-0.22841145229276575597e-3,-0.11404605562013443659e-5,0.81430290992322326296e-7,-0.71512447242755357629e-09,-0.13372664928000000000e-10,0.44461498336689298148e-12,0.00000000000000000000e+00}, + A9{0.58380635448407827360e00,0.49717469530842831182e-2,-0.23336001540009645365e-3,-0.51952064448608850822e-6,0.73596577815411080511e-7,-0.84020916763091566035e-09,-0.76700972702222222221e-11,0.36914462807972467044e-12,0.00000000000000000000e+00}, + A9{0.59281340237769489597e00,0.40343592069379730568e-2,-0.23477963738658326185e-3,0.34615944987790224234e-7,0.64832803248395814574e-7,-0.90329163587627007971e-09,-0.30421940400000000000e-11,0.29237386653743536669e-12,0.00000000000000000000e+00}, + A9{0.59994428743114271918e00,0.30976579788271744329e-2,-0.23308875765700082835e-3,0.51681681023846925160e-6,0.55694594264948268169e-7,-0.91719117313243464652e-09,0.53982743680000000000e-12,0.22050829296187771142e-12,0.00000000000000000000e+00}, + A9{0.60521224471819875444e00,0.21732138012345456060e-2,-0.22872428969625997456e-3,0.92588959922653404233e-6,0.46612665806531930684e-7,-0.89393722514414153351e-09,0.31718550353777777778e-11,0.15705458816080549117e-12,0.00000000000000000000e+00}, + A9{0.60865189969791123620e00,0.12708480848877451719e-2,-0.22212090111534847166e-3,0.12636236031532793467e-5,0.37904037100232937574e-7,-0.84417089968101223519e-09,0.49843180828444444445e-11,0.10355439441049048273e-12,0.00000000000000000000e+00}, + A9{0.61031580103499200191e00,0.39867436055861038223e-3,-0.21369573439579869291e-3,0.15339402129026183670e-5,0.29787479206646594442e-7,-0.77687792914228632974e-09,0.61192452741333333334e-11,0.60216691829459295780e-13,0.00000000000000000000e+00}, + A9{0.61027109047879835868e00,-0.43680904508059878254e-3,-0.20383783788303894442e-3,0.17421743090883439959e-5,0.22400425572175715576e-7,-0.69934719320045128997e-09,0.67152759655111111110e-11,0.26419960042578359995e-13,0.00000000000000000000e+00}, + A9{0.60859639489217430521e00,-0.12305921390962936873e-2,-0.19290150253894682629e-3,0.18944904654478310128e-5,0.15815530398618149110e-7,-0.61726850580964876070e-09,0.68987888999111111110e-11,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.60537899426486075181e00,-0.19790062241395705751e-2,-0.18120271393047062253e-3,0.19974264162313241405e-5,0.10055795094298172492e-7,-0.53491997919318263593e-09,0.67794550295111111110e-11,-0.17059208095741511603e-13,0.00000000000000000000e+00}, + A9{0.60071229457904110537e00,-0.26795676776166354354e-2,-0.16901799553627508781e-3,0.20575498324332621581e-5,0.51077165074461745053e-8,-0.45536079828057221858e-09,0.64488005516444444445e-11,-0.29311677573152766338e-13,0.00000000000000000000e+00}, + A9{0.59469361520112714738e00,-0.33308208190600993470e-2,-0.15658501295912405679e-3,0.20812116912895417272e-5,0.93227468760614182021e-9,-0.38066673740116080415e-09,0.59806790359111111110e-11,-0.36887077278950440597e-13,0.00000000000000000000e+00}, + A9{0.58742228631775388268e00,-0.39321858196059227251e-2,-0.14410441141450122535e-3,0.20743790018404020716e-5,-0.25261903811221913762e-8,-0.31212416519526924318e-09,0.54328422462222222221e-11,-0.40864152484979815972e-13,0.00000000000000000000e+00}, + A9{0.57899804200033018447e00,-0.44838157005618913447e-2,-0.13174245966501437965e-3,0.20425306888294362674e-5,-0.53330296023875447782e-8,-0.25041289435539821014e-09,0.48490437205333333334e-11,-0.42162206939169045177e-13,0.00000000000000000000e+00}, + A9{0.56951968796931245974e00,-0.49864649488074868952e-2,-0.11963416583477567125e-3,0.19906021780991036425e-5,-0.75580140299436494248e-8,-0.19576060961919820491e-09,0.42613011928888888890e-11,-0.41539443304115604377e-13,0.00000000000000000000e+00}, + A9{0.55908401930063918964e00,-0.54413711036826877753e-2,-0.10788661102511914628e-3,0.19229663322982839331e-5,-0.92714731195118129616e-8,-0.14807038677197394186e-09,0.36920870298666666666e-11,-0.39603726688419162617e-13,0.00000000000000000000e+00}, + A9{0.54778496152925675315e00,-0.58501497933213396670e-2,-0.96582314317855227421e-4,0.18434405235069270228e-5,-0.10541580254317078711e-7,-0.10702303407788943498e-09,0.31563175582222222222e-11,-0.36829748079110481422e-13,0.00000000000000000000e+00}, + A9{0.53571290831682823999e00,-0.62147030670760791791e-2,-0.85782497917111760790e-4,0.17553116363443470478e-5,-0.11432547349815541084e-7,-0.72157091369041330520e-10,0.26630811607111111111e-11,-0.33578660425893164084e-13,0.00000000000000000000e+00}, + A9{0.52295422962048434978e00,-0.65371404367776320720e-2,-0.75530164941473343780e-4,0.16613725797181276790e-5,-0.12003521296598910761e-7,-0.42929753689181106171e-10,0.22170894940444444444e-11,-0.30117697501065110505e-13,0.00000000000000000000e+00}, + A9{0.50959092577577886140e00,-0.68197117603118591766e-2,-0.65852936198953623307e-4,0.15639654113906716939e-5,-0.12308007991056524902e-7,-0.18761997536910939570e-10,0.18198628922666666667e-11,-0.26638355362285200932e-13,0.00000000000000000000e+00}, + A9{0.49570040481823167970e00,-0.70647509397614398066e-2,-0.56765617728962588218e-4,0.14650274449141448497e-5,-0.12393681471984051132e-7,0.92904351801168955424e-12,0.14706755960177777778e-11,-0.23272455351266325318e-13,0.00000000000000000000e+00}, + A9{0.48135536250935238066e00,-0.72746293327402359783e-2,-0.48272489495730030780e-4,0.13661377309113939689e-5,-0.12302464447599382189e-7,0.16707760028737074907e-10,0.11672928324444444444e-11,-0.20105801424709924499e-13,0.00000000000000000000e+00}, + A9{0.46662374675511439448e00,-0.74517177649528487002e-2,-0.40369318744279128718e-4,0.12685621118898535407e-5,-0.12070791463315156250e-7,0.29105507892605823871e-10,0.90653314645333333334e-12,-0.17189503312102982646e-13,0.00000000000000000000e+00}, + A9{0.45156879030168268778e00,-0.75983560650033817497e-2,-0.33045110380705139759e-4,0.11732956732035040896e-5,-0.11729986947158201869e-7,0.38611905704166441308e-10,0.68468768305777777779e-12,-0.14549134330396754575e-13,0.00000000000000000000e+00}, + A9{0.43624909769330896904e00,-0.77168291040309554679e-2,-0.26283612321339907756e-4,0.10811018836893550820e-5,-0.11306707563739851552e-7,0.45670446788529607380e-10,0.49782492549333333334e-12,-0.12191983967561779442e-13,0.00000000000000000000e+00}, + A9{0.42071877443548481181e00,-0.78093484015052730097e-2,-0.20064596897224934705e-4,0.99254806680671890766e-6,-0.10823412088884741451e-7,0.50677203326904716247e-10,0.34200547594666666666e-12,-0.10112698698356194618e-13,0.00000000000000000000e+00}, + A9{0.40502758809710844280e00,-0.78780384460872937555e-2,-0.14364940764532853112e-4,0.90803709228265217384e-6,-0.10298832847014466907e-7,0.53981671221969478551e-10,0.21342751381333333333e-12,-0.82975901848387729274e-14,0.00000000000000000000e+00}, + A9{0.38922115269731446690e00,-0.79249269708242064120e-2,-0.91595258799106970453e-5,0.82783535102217576495e-6,-0.97484311059617744437e-8,0.55889029041660225629e-10,0.10851981336888888889e-12,-0.67278553237853459757e-14,0.00000000000000000000e+00}, + A9{0.37334112915460307335e00,-0.79519385109223148791e-2,-0.44219833548840469752e-5,0.75209719038240314732e-6,-0.91848251458553190451e-8,0.56663266668051433844e-10,0.23995894257777777778e-13,-0.53819475285389344313e-14,0.00000000000000000000e+00}, + A9{0.35742543583374223085e00,-0.79608906571527956177e-2,-0.12530071050975781198e-6,0.68088605744900552505e-6,-0.86181844090844164075e-8,0.56530784203816176153e-10,-0.43120012248888888890e-13,-0.42372603392496813810e-14,0.00000000000000000000e+00}, + A9{0.34150846431979618536e00,-0.79534924968773806029e-2,0.37576885610891515813e-5,0.61419263633090524326e-6,-0.80565865409945960125e-8,0.55684175248749269411e-10,-0.95486860764444444445e-13,-0.32712946432984510595e-14,0.00000000000000000000e+00}, + A9{0.32562129649136346824e00,-0.79313448067948884309e-2,0.72539159933545300034e-5,0.55195028297415503083e-6,-0.75063365335570475258e-8,0.54281686749699595941e-10,-0.13545424295111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.30979191977078391864e00,-0.78959416264207333695e-2,0.10389774377677210794e-4,0.49404804463196316464e-6,-0.69722488229411164685e-8,0.52469254655951393842e-10,-0.16507860650666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.29404543811214459904e00,-0.78486728990364155356e-2,0.13190885683106990459e-4,0.44034158861387909694e-6,-0.64578942561562616481e-8,0.50354306498006928984e-10,-0.18614473550222222222e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.27840427686253660515e00,-0.77908279176252742013e-2,0.15681928798708548349e-4,0.39066226205099807573e-6,-0.59658144820660420814e-8,0.48030086420373141763e-10,-0.20018995173333333333e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.26288838011163800908e00,-0.77235993576119469018e-2,0.17886516796198660969e-4,0.34482457073472497720e-6,-0.54977066551955420066e-8,0.45572749379147269213e-10,-0.20852924954666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.24751539954181029717e00,-0.76480877165290370975e-2,0.19827114835033977049e-4,0.30263228619976332110e-6,-0.50545814570120129947e-8,0.43043879374212005966e-10,-0.21228012028444444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.23230087411688914593e00,-0.75653060136384041587e-2,0.21524991113020016415e-4,0.26388338542539382413e-6,-0.46368974069671446622e-8,0.40492715758206515307e-10,-0.21238627815111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.21725840021297341931e00,-0.74761846305979730439e-2,0.23000194404129495243e-4,0.22837400135642906796e-6,-0.42446743058417541277e-8,0.37958104071765923728e-10,-0.20963978568888888889e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.20239979200788191491e00,-0.73815761980493466516e-2,0.24271552727631854013e-4,0.19590154043390012843e-6,-0.38775884642456551753e-8,0.35470192372162901168e-10,-0.20470131678222222222e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.18773523211558098962e00,-0.72822604530339834448e-2,0.25356688567841293697e-4,0.16626710297744290016e-6,-0.35350521468015310830e-8,0.33051896213898864306e-10,-0.19811844544000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.17327341258479649442e00,-0.71789490089142761950e-2,0.26272046822383820476e-4,0.13927732375657362345e-6,-0.32162794266956859603e-8,0.30720156036105652035e-10,-0.19034196304000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.15902166648328672043e00,-0.70722899934245504034e-2,0.27032932310132226025e-4,0.11474573347816568279e-6,-0.29203404091754665063e-8,0.28487010262547971859e-10,-0.18174029063111111111e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.14498609036610283865e00,-0.69628725220045029273e-2,0.27653554229160596221e-4,0.92493727167393036470e-7,-0.26462055548683583849e-8,0.26360506250989943739e-10,-0.17261211260444444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.13117165798208050667e00,-0.68512309830281084723e-2,0.28147075431133863774e-4,0.72351212437979583441e-7,-0.23927816200314358570e-8,0.24345469651209833155e-10,-0.16319736960000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.11758232561160626306e00,-0.67378491192463392927e-2,0.28525664781722907847e-4,0.54156999310046790024e-7,-0.21589405340123827823e-8,0.22444150951727334619e-10,-0.15368675584000000000e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.10422112945361673560e00,-0.66231638959845581564e-2,0.28800551216363918088e-4,0.37758983397952149613e-7,-0.19435423557038933431e-8,0.20656766125421362458e-10,-0.14422990012444444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.91090275493541084785e-1,-0.65075691516115160062e-2,0.28982078385527224867e-4,0.23014165807643012781e-7,-0.17454532910249875958e-8,0.18981946442680092373e-10,-0.13494234691555555556e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.78191222288771379358e-1,-0.63914190297303976434e-2,0.29079759021299682675e-4,0.97885458059415717014e-8,-0.15635596116134296819e-8,0.17417110744051331974e-10,-0.12591151763555555556e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.65524757106147402224e-1,-0.62750311956082444159e-2,0.29102328354323449795e-4,-0.20430838882727954582e-8,-0.13967781903855367270e-8,0.15958771833747057569e-10,-0.11720175765333333333e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.53091065838453612773e-1,-0.61586898417077043662e-2,0.29057796072960100710e-4,-0.12597414620517987536e-7,-0.12440642607426861943e-8,0.14602787128447932137e-10,-0.10885859114666666667e-12,0.00000000000000000000e+00,0.00000000000000000000e+00}, + A9{0.40889797115352738582e-1,-0.60426484889413678200e-2,0.28953496450191694606e-4,-0.21982952021823718400e-7,-0.11044169117553026211e-8,0.13344562332430552171e-10,-0.10091231402844444444e-12,0.00000000000000000000e+00,0.00000000000000000000e+00} + }; + + auto w_im_y100 = [](T y100, T x) //purely scalar + { + if (eve::is_nan(x)) return x; + auto t = eve::dec(2*eve::frac(y100)); + int index = ((int) y100); + if (index < 97) return eve::reverse_horner(t, coefs[index]); + else + { + // use Taylor expansion for small x (|x| <= 0.0309...) + // (2/sqrt(pi)) * (x -2/3 x^3 , 4/15 x^5 -8/105 x^7, 16/945 x^9) + auto mx2(-eve::sqr(x)); + kumi::result::generate_t<5, r_t> a{1.1283791670955125739, 0.75225277806367504925, 0.30090111122547001970 + , 0.085971746064420005629, 0.016931216931216931217}; + return x*eve::reverse_horner(mx2, a); + } + }; + + auto taylor = [mz2, z](){ + using r_t = eve::element_type_t; + kumi::result::generate_t<5, r_t> a{1.1283791670955125739, 0.37612638903183752464, 0.11283791670955125739 + , 0.026866170645131251760, 0.0052239776254421878422}; + return eve::reverse_horner(mz2, a)*z; + }; + + if constexpr(eve::scalar_value) + { + auto w_im = [w_im_y100](real_t xx){ + const r_t ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) + auto signxx = sign(xx); + xx = eve::abs(xx); + if (xx > 45) { // continued-fraction expansion is faster + if (xx > real_t(5e7)) // 1-term expansion, important to avoid overflow + return ispi / (xx*signxx); + /* 5-term expansion (rely on compiler for CSE), simplified from: + ispi / (xx-0.5/(xx-1/(xx-1.5/(xx-2/xx)))) */ + return ispi*((xx*xx) * (xx*xx-real_t(4.5)) + 2) / (xx * ((xx*xx) * (xx*xx-5) + real_t(3.75)))*signxx; + } + return w_im_y100(100/(1+xx), xx)*signxx; + }; + + if(eve::is_eqz(y)) + { + return complex(eve::erf(x), y); //call to real implementation + } + else if(eve::is_eqz(x)) + { + return complex(x, (eve::sqr(y) > real_t(720) ? eve::inf(eve::as(y))*eve::sign(y) : eve::expx2(y) * w_im(y))); + } + if (real(mz2) < -750) + { + return complex(eve::signnz(x)); // underflow + } + /* Handle positive and negative x via different formulas, + using the mirror symmetries of w, to avoid overflow/underflow + problems from multiplying exponentially large and small quantities. */ + auto signx = eve::sign(x); + auto ax = eve::abs(x); + auto smallx = ax < 8e-2; + auto smally = eve::abs(y) < 1e-2; + if (smallx && smally) return taylor(); + /* don't use complex exp function, since that will produce spurious NaN + values when multiplying w in an overflow situation. */ + auto [mRe_z2, mIm_z2] = mz2; + auto [s, c] = eve::sincos(mIm_z2); + return oneminus(eve::exp(mRe_z2)*(complex(c, s)*faddeeva(complex(-y, x)*signx)))*signx; + } + else //simd + { + auto signx = eve::signnz(x); + auto w_im = [signx, w_im_y100](real_t x){ + const real_t ispi = real_t(0.56418958354775628694807945156); // 1 / sqrt(pi) + auto greater5e7 = [ispi, signx](auto xx){ + return ispi / (xx*signx); + }; + auto greater45 = [ispi, signx](auto xx){ + auto xx2 = xx*xx; + return signx*ispi*((xx2) * (xx2-real_t(4.5)) + 2) / (xx * ((xx2) * (xx2-5) + real_t(3.75))); + }; + auto remain = [w_im_y100, signx](auto xx){ + + if constexpr(eve::scalar_value) + { + return w_im_y100(100/(1+xx), xx)*signx; + } + else + { + return eve::detail::map(w_im_y100, 100/(1+xx), xx)*signx; + } + }; + + auto ax = eve::abs(x); + auto notdone = eve::true_(eve::as(ax)); + auto r = eve::nan(eve::as(x)); + if( eve::any(notdone) ) + { + notdone = next_interval(greater5e7, notdone, ax > real_t(5e7), r, ax); + if( eve::any(notdone) ) + { + notdone = next_interval(greater45, notdone, ax > real_t(45), r, ax); + if( eve::any(notdone) ) + { + if (eve::any(notdone)) r = remain(ax); + notdone = eve::false_(eve::as(r)); + } + } + } + return r; + }; + + auto signy = eve::signnz(y); + auto no_underflow = real(mz2) >= -750; + auto nfin = eve::is_not_finite(y); //|| is_nan(z); + auto r = if_else(no_underflow || is_nan( real(mz2)) + , complex(eve::nan(eve::as(real_t(0))), eve::nan(eve::as(real_t(0)))) + , complex(signx, real_t(0))); // treat underflow and nan + auto notdone = (no_underflow && !nfin) || eve::is_eqz(y); // no underflow + r = if_else(nfin, complex(eve::nan(eve::as(y)), eve::nan(eve::as(y))), r); + r = if_else(eve::is_eqz(x) && nfin, complex(real_t(0), y), r); + auto ax = eve::abs(x); + auto xsmall = ax < 8e-2; + auto ysmall = eve::abs(y) < 1e-2; + auto nully = [](auto x, auto y){ + return complex(eve::erf(x), y); + }; + auto nullx = [signy, w_im](auto x, auto y){ + return complex(x, eve::if_else(eve::sqr(y) > real_t(720), eve::inf(eve::as(y)), eve::expx2(y) * w_im(y)))*signy; + }; + + auto smallxy = [taylor](){ + return taylor(); + }; + + auto remain = [mz2, signx](auto x, auto y){ + auto [mRe_z2, mIm_z2] = mz2; + auto [s, c] = eve::sincos(mIm_z2); + auto zz = oneminus(eve::exp(mRe_z2)*(complex(c, s)*faddeeva(complex(-y, x)*signx)))*signx; + return zz; + }; + + if( eve::any(notdone) ) + { + notdone = next_interval(nully, notdone, eve::is_eqz(y), r, x, y); + if( eve::any(notdone) ) + { + notdone = next_interval(nullx, notdone, eve::is_eqz(x), r, x, y); + if( eve::any(notdone) ) + { + notdone = next_interval(smallxy, notdone, xsmall && ysmall, r); + if( eve::any(notdone) ) + { + notdone = last_interval(remain, notdone, r, x, y); + } + } + } + } + return r; + } + } + else + { + return cayley_extend(erf, z); + } + } +} diff --git a/include/kyosu/types/impl/faddeeva.hpp b/include/kyosu/types/impl/faddeeva.hpp new file mode 100644 index 00000000..16860fcc --- /dev/null +++ b/include/kyosu/types/impl/faddeeva.hpp @@ -0,0 +1,123 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include + +namespace kyosu::_ +{ + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using v_t = as_real_t; + using real_t = eve::element_type_t; + auto const sqrtpi = eve::sqrt_pi(eve::as()); + auto const iosqrtpi = complex(real_t(0), eve::rec(sqrtpi)); + + auto fexp = [iosqrtpi, sqrtpi](auto z){//Fourier expansion approximation + + constexpr real_t tauM = 12; + constexpr real_t tauM2= 144; + constexpr size_t maxN = 23; + constexpr std::array aN = { //Fourier coefficients + 2.758402332921771e-01, 2.245739552246158e-01, 1.594149382739117e-01, 9.866576641545419e-02, 5.324414078763941e-02 + , 2.505215000539365e-02, 1.027746567053954e-02, 3.676164332844847e-03, 1.146493641242233e-03, 3.117570150461976e-04 + , 7.391433429603010e-05, 1.527949342800837e-05, 2.753956608221073e-06, 4.327858781901246e-07, 5.930030408745884e-08 + , 7.084490307748205e-09, 7.379520635816785e-10, 6.702171606002010e-11, 5.307265163470805e-12, 3.664324113467642e-13 + , 2.205894944941035e-14, 1.157826862628556e-15, 5.298711429467307e-17}; + auto z1 = exp_i(tauM*z); + auto z2 = tauM2*sqr(z); + auto FE = sqrtpi/tauM*oneminus(z1)/z2; + for (size_t n = 1; n <= maxN; ++n) + FE += (aN[n-1]*dec(eve::sign_alternate(real_t(n))*z1)/(sqr(n*eve::pi(eve::as(tauM))) - z2)); + return iosqrtpi*tauM2*z*FE; + }; + + auto contfr = [iosqrtpi](auto z){ // the Laplace continued fraction approximation + constexpr std::array bN = {5.5000, 5.0000, 4.5000, 4.0000, 3.5000, + 3.0000, 2.5000, 2.0000, 1.5000, 1.0000, 0.5000}; + auto CF = bN[0]/z; + for (int k = 1; k <= 10; ++k) CF = bN[k]*rec(z - CF); + return iosqrtpi*rec(z - CF); + }; + + auto small_z= [sqrtpi](auto z){ + auto zP2=sqr(z); + auto zP4=sqr(zP2); + auto zP6=zP2*zP4; + + return (((6 - 6*zP2 + 3*zP4 - zP6)*(15*sqrtpi + complex(real_t(0), real_t(1))*z*(30 + 10*zP2 + 3*zP4)))/(90*sqrtpi)); + }; + + auto narr_band = [sqrtpi](auto z, auto aN, auto tauM){ // the narrow band + real_t tauM2 = sqr(tauM); + constexpr size_t maxN = 23; + auto z1 = cos(tauM*z); + auto z2 = tauM2*sqr(z); + Z NB{}; + for (size_t n = 1; n <= maxN; ++n) NB += (aN[n-1]*dec(eve::sign_alternate(real_t(n))*z1)/(sqr(n*eve::pi(eve::as(tauM))) - z2)); + return exp(-sqr(z)) - complex(real_t(0), real_t(1))*(dec(z1)/(tauM*z) - tauM2*z/sqrtpi*NB); + }; + + + auto smallim = [narr_band, small_z](auto z){ // approximation at small imag(z) + constexpr size_t maxN = 23; + constexpr std::array aN12 = { //Fourier coefficients + 2.758402332921771e-01, 2.245739552246158e-01, 1.594149382739117e-01, 9.866576641545419e-02, 5.324414078763941e-02 + , 2.505215000539365e-02, 1.027746567053954e-02, 3.676164332844847e-03, 1.146493641242233e-03, 3.117570150461976e-04 + , 7.391433429603010e-05, 1.527949342800837e-05, 2.753956608221073e-06, 4.327858781901246e-07, 5.930030408745884e-08 + , 7.084490307748205e-09, 7.379520635816785e-10, 6.702171606002010e-11, 5.307265163470805e-12, 3.664324113467642e-13 + , 2.205894944941035e-14, 1.157826862628556e-15, 5.298711429467307e-17}; + constexpr std::array aN12_1 = { //Fourier coefficients avoiding poles + 2.738693653257005e-01, 2.237253191645656e-01, 1.597109174949294e-01, 9.963269094255395e-02, 5.431463971403691e-02 + , 2.587496275171680e-02, 1.077185133049561e-02, 3.918760804860615e-03, 1.245818993506881e-03, 3.461058398397597e-04 + , 8.402542038611651e-05, 1.782626047239582e-05, 3.304894416742214e-06, 5.354300211460283e-07, 7.580461285556943e-08 + , 9.378563805151443e-09, 1.013969351293613e-09, 9.579902872635222e-11, 7.909429715194891e-12, 5.706594634502440e-13 + , 3.597962756907448e-14, 1.982367142666383e-15, 9.544634564831883e-17 }; + + auto x = real(z); + auto ind_0 = abs(x)<5e-3; + auto ind_poles = eve::false_(eve::as(x)); + for(int k = 1; k <= 23; ++k) ind_poles = ind_poles || (abs(x - k*eve::pi(eve::as(x))/12)<1e-4); + return if_else(!ind_0, if_else(ind_poles + , narr_band(z, aN12_1, real_t(12.1)) + , narr_band(z, aN12, real_t(12)) + ) + , small_z(z)); + }; + + auto indneg = eve::is_ltz(imag(z)); + z = if_else(indneg, conj(z), z); + auto r = complex(eve::nan(eve::as(real(z))), eve::nan(eve::as(real(z)))); // nan case treated here + r = if_else(eve::is_infinite(real(z)), Z{}, r); + auto notdone = is_finite(z); + if (eve::any(notdone)) + { + auto indin = sqr_abs(z) <= real_t(64); + auto indband = indin && (imag(z) < real_t(5.0e-3)); + notdone = next_interval(smallim, notdone, indband, r, z); + if( eve::any(notdone) ) + { + notdone = next_interval(fexp, notdone, indin, r, z); + if( eve::any(notdone) ) + { + notdone = next_interval(contfr, notdone, notdone, r, z); + } + } + } + return if_else(indneg, conj(r), r); ; + } + else + { + return cayley_extend(faddeeva, z); + } + } +} diff --git a/include/kyosu/types/impl/invtrig.hpp b/include/kyosu/types/impl/invtrig.hpp new file mode 100644 index 00000000..e846d2be --- /dev/null +++ b/include/kyosu/types/impl/invtrig.hpp @@ -0,0 +1,625 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::_ +{ + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // This implementation is a simd transcription and adaptation of the boost_math code + // which itself is a transcription of the pseudo-code in: + // + // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling." + // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang. + // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997. + // + auto [a0r, a0i] = a0; + using rtype = decltype(a0r); + const rtype a_crossover(1.5); + const rtype b_crossover(0.6417); + auto ltzra0 = eve::is_negative(a0r); + auto ltzia0 = eve::is_negative(a0i); + + // + // Begin by insuring real(a0) >= 0 and imag(a0) >= 0 : + // + + rtype x = eve::abs(a0r); + rtype y = eve::abs(a0i); + rtype proper_real = eve::asin(x); + auto lexone = (x <= eve::one(eve::as(a0r))); + auto is_proper_real = is_real(a0) && lexone; + + auto s_min = eve::sqrtsmallestposval(eve::as(x))*4; + auto s_max = eve::sqrtvalmax(eve::as(x))/8; + rtype xp1 = eve::inc(x); + rtype xm1 = eve::dec(x); + auto not_in_safe_zone = (((x > s_max) || (x < s_min)) || ((y > s_max) || (y < s_min))); + //compute for safe zone + rtype r, i; + rtype yy = eve::sqr(y); + rtype tr = eve::sqrt(eve::sqr(xp1) + yy);//hypot for pedantic ? + rtype ts = eve::sqrt(eve::sqr(xm1) + yy);//hypot for pedantic ? + rtype a = eve::average(tr, ts); + rtype b = x/a; + //compute r for b > b_crossover + rtype apx = a + x; + r = eve::if_else(lexone, + eve::atan(x/eve::sqrt(eve::half(eve::as(a0r))*apx*(yy/(tr+xp1)+(ts-xm1)))), + eve::atan(x/(y*eve::sqrt(eve::half(eve::as(a0r))*(apx/(tr+xp1)+apx/(ts+xm1))))) + ); + // r is computed + r = eve::if_else((b <= b_crossover), r, eve::asin(b)); + //compute am 1 temporary for i for a <= a_crossover + rtype tmp = yy/(tr+xp1); + rtype am1 = eve::if_else(lexone, + eve::average(tmp, yy/(ts-xm1)), + eve::average(tmp, (ts+xm1))); + i = eve::if_else((a <= a_crossover), + eve::log1p(am1 + eve::sqrt(am1*(eve::inc(a)))), + eve::log(a + eve::sqrt(eve::dec(eve::sqr(a)))) + ); + // i is computed + //compute for exception zone + if (eve::any(not_in_safe_zone)) + { + auto zone1 = (y <= eve::eps(eve::as(a0r))*eve::abs(xm1)); + if (eve::any(eve::logical_and(zone1, not_in_safe_zone))) + { + rtype rr = eve::if_else(lexone, eve::asin(x), eve::pio_2(eve::as(a0r))); + rtype ii = eve::if_else(lexone, y/eve::sqrt(xp1*xm1), + eve::if_else((eve::valmax(eve::as(a0r))/xp1 > xm1), + eve::log1p(xm1 + eve::sqrt(xp1*xm1)), + eve::log_2(eve::as(a0r)) + eve::log(x) + ) + ); + r = eve::if_else(zone1, rr, r); + i = eve::if_else(zone1, ii, i); + } + auto zone2 = (y <= s_min); + auto not_treated = eve::logical_notand(zone1, not_in_safe_zone); + if (eve::any(eve::logical_and(zone2, not_treated))) + { + r = eve::if_else(zone2, eve::pio_2(eve::as(a0r)) - eve::sqrt(y), r); + i = eve::if_else(zone2, eve::sqrt(y), i); + } + auto zone3 = (eve::dec(eve::eps(eve::as(a0r))*y) >= x); + not_treated = eve::logical_notand(zone2, not_treated); + if (eve::any(eve::logical_and(zone3, not_treated))) + { + r = eve::if_else(zone3, x/y, r); + i = eve::if_else(zone3, eve::log_2(eve::as(a0r)) + eve::log(y), i); + } + auto zone4 = (x > eve::one(eve::as(a0r))); + not_treated = eve::logical_notand(zone3, not_treated); + if (eve::any(eve::logical_and(zone4, not_treated))) + { + r = eve::if_else(zone4, eve::atan(x/y), r); + i = eve::if_else(zone4, eve::log_2(eve::as(a0r)) + eve::log(y) + + eve::half(eve::as(a0r))*eve::log1p(eve::sqr(x/y)), i); + } + not_treated = eve::logical_notand(zone4, not_treated); + if (eve::any(not_treated)) + { + rtype aa = eve::sqrt(eve::inc(eve::sqr(y))); + r = eve::if_else(not_treated, x/aa, r); + i = eve::if_else(not_treated,eve::half(eve::as(a0r))*eve::log1p(2*y*(y+aa)), i); + } + } + if (eve::any(is_not_finite(a0))) + { + auto nanx = eve::is_nan(x); + auto nany = eve::is_nan(y); + auto infx = (x == eve::inf(eve::as(a0r))) ; + auto infy = (y == eve::inf(eve::as(a0r))) ; + if (eve::any(nanx)) + { + r = eve::if_else(nanx, x, r); + r = eve::if_else(nanx && infy, x, r); + i = eve::if_else(nanx, x, i); + i = eve::if_else(nanx && infy, y, i); + } + if (eve::any(nany)) + { + auto isimag = is_imag(a0); + r = eve::if_else(isimag && nany, eve::zero, r); + r = eve::if_else(eve::logical_and(nany, infx),y, r); + i = eve::if_else(isimag && nany, eve::allbits, i); + i = eve::if_else(nany && infx, x, i); + } + auto test = eve::logical_notand(eve::logical_or(nanx, nany), infx); + if (eve::any(test)) + { + r = eve::if_else(infx && test, + eve::if_else(infy, eve::pio_4(eve::as(a0r)), eve::pio_2(eve::as(a0r))), + r); + } + test = eve::logical_notand(eve::is_nan(x) || eve::is_nan(y), + eve::logical_andnot(infy, infx)); + r = eve::if_else(test,eve::zero,r); + } + // use proper real results + + r = eve::if_else(is_proper_real, proper_real, r); + i = eve::if_else(is_proper_real, eve::zero, i); + // restore signs + r = eve::if_else(ltzra0, -r, r); + i = eve::if_else(ltzia0, -i, i); + return complex(r, i); + } + else + { + return cayley_extend(asin, a0); + } + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // This implementation is a simd transcription and adaptation of the boost_math code + // which itself is a transcription of the pseudo-code in: + // + // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling." + // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang. + // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997. + // + auto [a0r, a0i] = a0; + using rtype = decltype(a0r); + const rtype a_crossover(1.5); + const rtype b_crossover(0.6417); + auto ltzra0 = eve::is_negative(a0r); + auto gtzia0 = eve::is_positive(a0i); + // + // Begin by insuring a0r >= 0 and imag(a0) >= 0 : + // + rtype x = eve::abs(a0r); + rtype y = eve::abs(a0i); + rtype proper_real = eve::acos(x); + auto lexone = (x <= eve::one(eve::as(x))); + auto is_proper_real = eve::logical_and(is_real(a0), lexone); + + auto s_min = eve::sqrtsmallestposval(eve::as(x))*4; + auto s_max = eve::sqrtvalmax(eve::as(x))/8; + rtype xp1 = eve::inc(x); + rtype xm1 = eve::dec(x); + auto not_in_safe_zone = (((x > s_max) || (x < s_min)) || ((y > s_max) || (y < s_min))); + //compute for safe zone + rtype r, i; + rtype yy = eve::sqr(y); + rtype tr = eve::sqrt(eve::sqr(xp1) + yy); //hypot for pedantic ? + rtype ts = eve::sqrt(eve::sqr(xm1) + yy); //hypot for pedantic ? + rtype a = eve::average(tr, ts); + rtype b = x/a; + //compute r for b > b_crossover + rtype apx = a + x; + r = eve::if_else(lexone, + eve::atan(eve::sqrt(eve::half(eve::as(x))*apx*(yy/(tr+xp1)+(ts-xm1)))/x), + eve::atan((y*eve::sqrt(eve::half(eve::as(x))*(apx/(tr+xp1)+apx/(ts+xm1))))/x) + ); + // r is computed + r = eve::if_else((b <= b_crossover), eve::acos(b), r); + //compute am1 temporary for i for a <= a_crossover + rtype tmp = yy/(tr+xp1); + rtype am1 = eve::if_else(lexone, + eve::average(tmp, yy/(ts-xm1)), + eve::average(tmp, (ts+xm1))); + i = eve::if_else((a <= a_crossover), + eve::log1p(am1 + eve::sqrt(am1*(eve::inc(a)))), + eve::log(a + eve::sqrt(eve::dec(eve::sqr(a)))) + ); + // i is computed + //compute for exception zone + if (eve::any(not_in_safe_zone)) + { + auto zone1 = (y <= eve::eps(eve::as(x))*eve::abs(xm1)); + if (eve::any(eve::logical_and(zone1, not_in_safe_zone))) + { + rtype rr = eve::if_else(lexone, proper_real, eve::zero); + rtype ii = eve::if_else(lexone, y/eve::sqrt(-xp1*xm1), + eve::if_else((eve::valmax(eve::as(x))/xp1 > xm1), + eve::log1p(xm1 + eve::sqrt(xp1*xm1)), + eve::log_2(eve::as(x)) + eve::log(x) + ) + ); + r = eve::if_else(zone1, rr, r); + i = eve::if_else(zone1, ii, i); + } + auto zone2 = (y <= s_min); + auto not_treated = eve::logical_notand(zone1, not_in_safe_zone); + if (eve::any(eve::logical_and(zone2, not_treated))) + { + rtype sqrty = eve::sqrt(y); + r = eve::if_else(zone2, sqrty, r); + i = eve::if_else(zone2, sqrty, i); + } + auto zone3 = (eve::dec(eve::eps(eve::as(x))*y) >= x); + not_treated = eve::logical_notand(zone2, not_treated); + if (eve::any(eve::logical_and(zone3, not_treated))) + { + r = eve::if_else(zone3, eve::pio_2(eve::as(x)), r); + i = eve::if_else(zone3, eve::log_2(eve::as(x)) + eve::log(y), i); + } + auto zone4 = (x > eve::one(eve::as(x))); + not_treated = eve::logical_notand(zone3, not_treated); + if (eve::any(eve::logical_and(zone4, not_treated))) + { + r = eve::if_else(zone4, eve::atan(y/x), r); + i = eve::if_else(zone4, eve::log_2(eve::as(x)) + eve::log(y) + eve::half(eve::as(x))*eve::log1p(eve::sqr(x/y)), i); + } + not_treated = eve::logical_notand(zone4, not_treated); + if (eve::any(not_treated)) + { + rtype aa = eve::sqrt(eve::inc(sqr(y))); + r = eve::if_else(not_treated, eve::pio_2(eve::as(x)), r); + i = eve::if_else(not_treated, eve::half(eve::as(x))*eve::log1p(2*y*(y+aa)), i); + } + } + if (eve::any(kyosu::is_not_finite(a0))) + { + auto nanx = eve::is_nan(x); + auto nany = eve::is_nan(y); + auto infx = (x == eve::inf(eve::as(x))) ; + auto infy = (y == eve::inf(eve::as(x))) ; + if (eve::any(infx)) + { + r = eve::if_else(infx, eve::zero, r); + i = eve::if_else(infx, eve::inf(eve::as(x)), i); + r = eve::if_else(eve::logical_and(infx, infy), eve::pio_4(eve::as(x)), r); + i = eve::if_else(eve::logical_and(infx, infy), eve::inf(eve::as(x)), i); + + r = eve::if_else(eve::logical_and(infx, nany), y, r); + i = eve::if_else(eve::logical_and(infx, nany), eve::minf(eve::as(x)), i); + } + if (eve::any(nanx)) + { + r = eve::if_else(nanx, x, r); + i = eve::if_else(nanx, x, i); + i = eve::if_else(eve::logical_and(nanx, infy), y, i); + } + auto test = eve::logical_notand(eve::logical_or(infx, nanx), infy); + if (eve::any(test)) + { + r = eve::if_else(eve::logical_and(infy, test), eve::pio_2(eve::as(x)), r); + i = eve::if_else(eve::logical_and(infy, test), y, i); + } + test = eve::logical_notand(eve::logical_or(infx, nanx), nany); + r = eve::if_else(test,eve::if_else(is_imag(a0), eve::pio_2(eve::as(x)), y), r); + i = eve::if_else(test,y,i); + } + // use proper real results + r = eve::if_else(is_proper_real, proper_real, r); + i = eve::if_else(is_proper_real, eve::zero, i); + // restore signs + r = eve::if_else(ltzra0, eve::pi(eve::as(x))-r, r); + i = eve::if_else(gtzia0, -i, i); + return complex(r, i); + } + else + { + return cayley_extend(acos, a0); + } + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // acosh(a0) = +/-i acos(a0) + // Choosing the sign of multiplier to give real(acosh(a0)) >= 0 + // we have compatibility with C99. + auto [r, i] = kyosu::acos(a0); + auto lez = eve::is_negative(i);; + auto res = complex(-i, r); + res = eve::if_else(lez, res, -res); + auto nani = is_nan(i); + if (eve::any(nani)) + return eve::if_else(nani && eve::is_finite(r) + , complex(eve::nan(eve::as(r)), eve::nan(eve::as(r))) + , res); + else + return res; + } + else + { + return cayley_extend(acosh, a0); + } + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + if constexpr(concepts::complex ) + { + auto [r, i] = a0; + auto [r1, i1] = kyosu::asin(complex(-i, r)); + return complex(i1, -r1); // -(eve::i*asin(eve::i*a0)); + } + else + { + return cayley_extend(asinh, a0); + } + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // This implementation is a simd (i.e. no branch) transcription and adaptation of the + // boost_math code which itself is a transcription of the pseudo-code in: + // + // Eric W. Weisstein. "Inverse Hyperbolic Tangent." + // From MathWorld--A Wolfram Web Resource. + // http://mathworld.wolfram.com/InverseHyperbolicTangent.html + // + // Also: The Wolfram Functions Site, + // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/ + // + // Also "Abramowitz and Stegun. Handbook of Mathematical Functions." + // at : http://jove.prohosting.com/~skripty/toc.htm + // + auto [a0r, a0i] = a0; + auto realinf = eve::is_eqz(a0i) && eve::is_infinite(a0r); + using rtype = decltype(a0r); + const rtype alpha_crossover(0.3); + auto ltzra0 = eve::is_ltz(a0r); + auto ltzia0 = eve::is_ltz(a0i); + auto s_min = eve::sqrtsmallestposval(eve::as(a0r))*2; + auto s_max = eve::sqrtvalmax(eve::as(a0r))/2; + rtype const two = rtype(2); + rtype inf = eve::inf(eve::as(a0r)); + rtype x = eve::abs(a0r); + rtype y = eve::abs(a0i); + auto special = eve::is_eqz(y) && (x < eve::one(eve::as(a0r))); + auto sr = eve::atanh(a0r); + if (eve::all(special)) { + return complex(sr, eve::zero(eve::as(sr))); + } + + rtype r = eve::zero(eve::as(a0r)); + rtype i = eve::zero(eve::as(a0r)); + auto gtxmax = (x > s_max); + auto ltxmin = (x < s_min); + auto gtymax = (y > s_max); + auto ltymin = (y < s_min); + rtype xx = eve::sqr(x); + rtype yy = eve::sqr(y); + rtype sqrabs = xx + yy; + + auto not_in_safe_zone = ((gtxmax || ltxmin) || (gtymax || ltymin)); + if(eve::any(not_in_safe_zone)) + { + //treat underflow or overflow + // one or both of x and y are small, calculate divisor carefully: + rtype div = eve::one(eve::as(a0r)); + div += eve::if_else(ltxmin, xx, eve::zero); + div += eve::if_else(ltxmin, yy, eve::zero); + + rtype alpha = x/div; + alpha += alpha; + + auto test = gtymax; + // big y, medium x, divide through by y: + rtype tmp_alpha = (two*x/y) / (y + xx/y); + // small x and y, whatever alpha is, it's too small to calculate: + tmp_alpha = eve::if_else(x > eve::one(eve::as(a0r)), tmp_alpha, eve::zero); + alpha = eve::if_else(test && (x > eve::one(eve::as(a0r))), tmp_alpha, alpha); + + test = eve::logical_andnot(gtxmax, test); + + // big x small y, as above but neglect y^2/x: + tmp_alpha = two/x; + // big x: divide through by x: + tmp_alpha = eve::if_else((y > eve::one(eve::as(a0r))), two / (x + y*y/x), tmp_alpha); + // big x and y: divide alpha through by x*y: + tmp_alpha = eve::if_else(gtymax, (two/y) / (x/y + y/x), tmp_alpha); + // x or y are infinite: the result is 0 + tmp_alpha = eve::if_else((y == inf) || (x == inf), eve::zero, tmp_alpha); + + alpha = eve::if_else(test, tmp_alpha, alpha); + r = eve::if_else((alpha < alpha_crossover), + eve::log1p(alpha) - eve::log1p(-alpha), + eve::log(inc(two*x + xx)) - eve::log(sqr(dec(x))) + ); + test = (x == eve::one(eve::as(a0r))) && ltymin; + r = eve::if_else(test, -(two*(eve::log(y) - eve::log_2(eve::as(a0r)))), r); + r *= rtype(0.25); + //compute the imag part + // y^2 is negligible: + i = eve::atan2(two*y, eve::oneminus(xx)); + i = eve::if_else(gtymax || gtxmax, eve::pi(eve::as(a0r)), i); + rtype tmp_i = eve::if_else(ltymin, eve::atan2(two*y, eve::one(eve::as(a0r))), + eve::atan2(two*y, eve::oneminus(yy))); + i = eve::if_else(ltxmin, tmp_i, i); + } + auto test0 = (inf == x) && (inf == y); + if(eve::any(test0)) + { + //inf x, inf y + r = eve::if_else(test0, eve::zero, r); + i = eve::if_else(test0, eve::pi(eve::as(a0r)), r); + } + auto test = kyosu::is_nan(a0); + + if(eve::any(test)) + { + //nan x, inf y + r = eve::if_else(eve::is_nan(x) && (y == inf), eve::zero, r); + i = eve::if_else(eve::is_nan(x) && (y == inf), eve::pi(eve::as(a0r)), r); + + r = eve::if_else(is_nan(y) && (x == inf), eve::zero, r); + i = eve::if_else(is_nan(y) && (x == inf), y, i); + + r = eve::if_else(is_nan(y) && eve::is_eqz(x), eve::zero, r); + i = eve::if_else(is_nan(y) && is_eqz(x), eve::allbits, i); + } + //compute for safe zone::one + // the real part is given by: + // + // eve::real(atanh(z)) == eve::log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x)) + // + // however, when x is either large (x > 1/e) or very small + // (x < e) then this effectively simplifies + // to log(1), leading to wildly inaccurate results. + // by dividing the above (top and bottom) by (1 + x^2 + y^2) we get: + // + // eve::real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2)))) + // + // which is much more sensitive to the value of x, when x is not near 1 + // (remember we can compute log(1+x) for small x very accurately). + // + // the cross-over from eve::one method to the other has to be determined + // experimentally, the value used below appears correct to within a + // factor of 2 (and there are larger errors from other parts + // of the input domain anyway). + // + rtype alpha = x*two / (eve::inc(sqrabs)); + rtype sqrxm1 = eve::sqr(eve::dec(x)); + rtype tmp_r = eve::if_else((alpha < alpha_crossover), + eve::log1p(alpha) - eve::log1p(-alpha), + eve::log1p(x+x + sqrabs) - eve::log(sqrxm1 + yy) + )*rtype(0.25); + r = eve::if_else(not_in_safe_zone, r, tmp_r); + + // compute the imag part + i = eve::if_else(not_in_safe_zone, + i, + eve::atan2(y+y, (eve::oneminus(sqrabs))) + )*eve::half(eve::as(a0r)); + + r = eve::if_else( ltzra0,-r, r); + i = eve::if_else(eve::is_infinite(y), eve::pio_2(eve::as(a0r))*eve::sign(y), i); + i = eve::if_else( ltzia0,-i, i); + r = eve::if_else(realinf, eve::zero(eve::as(a0r)), r); + i = eve::if_else(realinf, -eve::sign(a0r)*eve::pio_2(eve::as(a0r)), i); + r = eve::if_else(special, sr, r); + i = eve::if_else(special, eve::zero, i); + return complex(r, i); + } + else + { + return cayley_extend(atanh, a0); + } + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // C99 definition here; atan(a0) = -i atanh(ia0): + auto [r, i] = a0; + auto [r1, i1] = kyosu::atanh(complex(-i, r)); + return complex(i1, -r1); + } + else + { + return cayley_extend(atan, a0); + } + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& z) noexcept + { + return kyosu::acos(kyosu::rec(z)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& z) noexcept + { + return kyosu::asin(kyosu::rec(z)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& z) noexcept + { + return kyosu::acosh(kyosu::rec(z)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& z) noexcept + { + return kyosu::asinh(kyosu::rec(z)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::atan(kyosu::rec(a0)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::atanh(kyosu::rec(a0)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::radinpi(kyosu::acos(a0)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::radinpi(kyosu::acot(a0)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::radinpi(kyosu::acsc(a0)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::radinpi(kyosu::asec(a0)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::radinpi(kyosu::asin(a0)); + } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& a0) noexcept + { + return kyosu::radinpi(kyosu::atan(a0)); + } + +} diff --git a/include/kyosu/types/impl/math.hpp b/include/kyosu/types/impl/math.hpp index 5523f4ea..42a7490b 100644 --- a/include/kyosu/types/impl/math.hpp +++ b/include/kyosu/types/impl/math.hpp @@ -10,6 +10,7 @@ #include #include #include +#include namespace kyosu::_ { @@ -29,11 +30,7 @@ namespace kyosu::_ } else { - auto p = pure(z); - auto az = abs(p); - auto r = exp(real(z)); - auto w = r*eve::sinc(az); - return r*eve::cos(az) + w*p; + return cayley_extend(exp, z); } } @@ -60,10 +57,7 @@ namespace kyosu::_ } else { - auto p = pure(z); - auto az = abs(p); - auto c = expm1(complex(real(z), az)); - return real(c) + ipart(c)*sign(p); + return cayley_extend(expm1, z); } } @@ -160,18 +154,19 @@ namespace kyosu::_ } else { - auto az = kyosu::abs(z); - auto v = kyosu::pure(z); - auto s = kyosu::real(z); - auto z1 = (eve::acos(s/az)/abs(v))*v+eve::log(az); - auto tmp = kyosu::if_else( kyosu::is_real(z) - , kyosu::log(kyosu::real(z)) - , z1 - ); - return kyosu::if_else( kyosu::is_eqz(z) - , eve::minf(eve::as(az)) - , tmp - ); + return cayley_extend(log, z); +// auto az = kyosu::abs(z); +// auto v = kyosu::pure(z); +// auto s = kyosu::real(z); +// auto z1 = (eve::acos(s/az)/abs(v))*v+eve::log(az); +// auto tmp = kyosu::if_else( kyosu::is_real(z) +// , kyosu::log(kyosu::real(z)) +// , z1 +// ); +// return kyosu::if_else( kyosu::is_eqz(z) +// , eve::minf(eve::as(az)) +// , tmp +// ); } } @@ -198,8 +193,9 @@ namespace kyosu::_ } else { - using e_t = eve::underlying_type_t; - return log(z)*eve::invlog_10(eve::as()); + return cayley_extend(log10, z); +// using e_t = eve::underlying_type_t; +// return log(z)*eve::invlog_10(eve::as()); } } @@ -225,9 +221,10 @@ namespace kyosu::_ } else { - using e_t = eve::underlying_type_t; - return log(z)*eve::invlog_2(eve::as()); - } + return cayley_extend(log2, z); +// using e_t = eve::underlying_type_t; +// return log(z)*eve::invlog_2(eve::as()); + } } template @@ -247,20 +244,21 @@ namespace kyosu::_ } else { - auto incz = inc(z); - auto az = kyosu::abs(incz); - auto az2 = kyosu::sqr_abs(z) + 2*real(z); - auto v = kyosu::pure(z); - auto s = kyosu::real(incz); - auto z1 = (eve::acos(s/az)/abs(v))*v+ eve::half(eve::as())*eve::log1p(az2); - auto tmp = kyosu::if_else( kyosu::is_real(z) - , kyosu::log(kyosu::real(z)) - , z1 - ); - return kyosu::if_else( kyosu::is_eqz(z) - , eve::minf(eve::as(az)) - , tmp - ); + return cayley_extend(log1p, z); +// auto incz = inc(z); +// auto az = kyosu::abs(incz); +// auto az2 = kyosu::sqr_abs(z) + 2*real(z); +// auto v = kyosu::pure(z); +// auto s = kyosu::real(incz); +// auto z1 = (eve::acos(s/az)/abs(v))*v+ eve::half(eve::as())*eve::log1p(az2); +// auto tmp = kyosu::if_else( kyosu::is_real(z) +// , kyosu::log(kyosu::real(z)) +// , z1 +// ); +// return kyosu::if_else( kyosu::is_eqz(z) +// , eve::minf(eve::as(az)) +// , tmp +// ); } } @@ -314,13 +312,14 @@ namespace kyosu::_ } else { - auto r = kyosu::abs(z); - auto theta = eve::acos(real(z)/r); - auto u = kyosu::sign(kyosu::pure(z)); - auto [s, c] = eve::sincos(theta*eve::half(eve::as(theta))); - auto res = u*s; - kyosu::real(res) = c; - return kyosu::if_else(eve::is_eqz(r), eve::zero(eve::as(z)), res*eve::sqrt(r)); + return cayley_extend(sqrt, z); +// auto r = kyosu::abs(z); +// auto theta = eve::acos(real(z)/r); +// auto u = kyosu::sign(kyosu::pure(z)); +// auto [s, c] = eve::sincos(theta*eve::half(eve::as(theta))); +// auto res = u*s; +// kyosu::real(res) = c; +// return kyosu::if_else(eve::is_eqz(r), eve::zero(eve::as(z)), res*eve::sqrt(r)); } } @@ -361,18 +360,25 @@ namespace kyosu::_ { if constexpr( eve::unsigned_value ) { - C0 base = c0; - C1 expo = c1; - auto const o = eve::one(eve::as()); - r_t result(o); - while(true) + if constexpr(kyosu::concepts::complex) { - if (eve::all(eve::is_eqz(expo))) break; - result = kyosu::if_else(eve::is_odd(expo), result*base, o); - expo = (expo >> 1); - base = kyosu::sqr(base); + C0 base = c0; + C1 expo = c1; + auto const o = eve::one(eve::as()); + r_t result(o); + while(true) + { + if (eve::all(eve::is_eqz(expo))) break; + result = kyosu::if_else(eve::is_odd(expo), result*base, o); + expo = (expo >> 1); + base = kyosu::sqr(base); + } + return result; + } + else + { + return cayley_extend(pow, c0, c1); } - return result; } else { @@ -414,8 +420,6 @@ namespace kyosu::_ } else if constexpr( kyosu::concepts::complex)// c0 and c1 are complex { -// auto rc1 = real(c1); -// auto ic1 = imag(c1); auto [rc1, ic1] = c1; auto lc0 = kyosu::log_abs(c0); auto argc0 = kyosu::arg(c0); @@ -434,16 +438,27 @@ namespace kyosu::_ } else { - auto cc0 = kyosu::convert(c0, eve::as()); - auto cc1 = kyosu::convert(c1, eve::as()); + if constexpr(eve::floating_value) //c0 cayley c1 real + { + return cayley_extend(pow, c0, c1); + } + else if constexpr(eve::floating_value)//c1 cayley c0 real + { + return cayley_extend_rev(pow, c0, c1); + } + else + { + auto cc0 = kyosu::convert(c0, eve::as()); + auto cc1 = kyosu::convert(c1, eve::as()); - auto r = kyosu::exp(kyosu::log(cc0)*cc1); - return kyosu::if_else (kyosu::is_eqz(cc0) - , eve::if_else(kyosu::is_eqz(cc1) + auto r = kyosu::exp(kyosu::log(cc0)*cc1); + return kyosu::if_else (kyosu::is_eqz(cc0) + , eve::if_else(kyosu::is_eqz(cc1) , eve::one(eve::as()) , eve::zero(eve::as())) - , r - ); + , r + ); + } } } diff --git a/include/kyosu/types/impl/predicates.hpp b/include/kyosu/types/impl/predicates.hpp index 7995325f..d8cc2aee 100644 --- a/include/kyosu/types/impl/predicates.hpp +++ b/include/kyosu/types/impl/predicates.hpp @@ -31,8 +31,15 @@ namespace kyosu::_ KYOSU_FORCEINLINE constexpr auto dispatch(eve::tag_of const&, C c) noexcept { - get<0>(c) = eve::zero(eve::as(get<0>(c))); + if constexpr(kyosu::concepts::complex) + { + return eve::is_eqz(ipart(c)); + } + else + { + get<0>(c) = eve::zero(eve::as(get<0>(c))); return kumi::all_of(c, [](auto const& e) { return eve::is_eqz(e); }); + } } template diff --git a/include/kyosu/types/impl/special.hpp b/include/kyosu/types/impl/special.hpp new file mode 100644 index 00000000..415c2cb7 --- /dev/null +++ b/include/kyosu/types/impl/special.hpp @@ -0,0 +1,477 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include + +namespace kyosu::_ +{ + //===------------------------------------------------------------------------------------------- + // Unary functions : tgamma + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // 15 sig. digits for 0<=real(z)<=171 + // coeffs should sum to about g*g/2+23/24 + // + using r_t = eve::element_type_t>; + auto g=r_t(607)/r_t(128); + // best results when 4<=g<=5 + constexpr int N = 15; + std::array c = + { 0.99999999999999709182, + 57.156235665862923517, + -59.597960355475491248, + 14.136097974741747174, + -0.49191381609762019978, + .33994649984811888699e-4, + .46523628927048575665e-4, + -.98374475304879564677e-4, + .15808870322491248884e-3, + -.21026444172410488319e-3, + .21743961811521264320e-3, + -.16431810653676389022e-3, + .84418223983852743293e-4, + -.26190838401581408670e-4, + .36899182659531622704e-5 + }; + + //Num Recipes used g=5 with 7 terms + //for a less effective approximation + + auto negra0 = eve::is_negative(real(a0)); + auto z = if_else(negra0, -a0, a0); + z = dec(z); + auto zh = z+eve::half(eve::as()); + auto zgh=zh+g; + //trick for avoiding FP overflow above z=141 + auto zp=pow(zgh,(zh*eve::half(eve::as()))); + auto ss = Z{}; + for(int pp = N-1; pp >= 1; --pp){ + ss+= c[pp]/(z+pp); + } + auto sq2pi = r_t(2.5066282746310005024157652848110); + auto f=(sq2pi*(c[0]+ss))*((zp*exp(-zgh))*zp); + auto o = eve::one(eve::as()); + f = if_else(is_eqz(z) || z == o, o, f); + //adjust for negative real parts + auto reala0 = is_real(a0); + if(eve::any(negra0)) + { + f = if_else(negra0, rec(-eve::inv_pi(eve::as(real(a0)))*a0*f*sinpi(a0)), eve::zero); + f = if_else (negra0 && reala0 && eve::is_flint(real(a0)), complex(eve::nan(eve::as(sq2pi)), eve::inf(eve::as(sq2pi))), f); + } + f = if_else (reala0, complex(eve::tgamma(real(a0))), f); + f = if_else (eve::is_nan(real(f)), complex(eve::nan(eve::as(sq2pi)), eve::inf(eve::as(sq2pi))), f); + f = if_else (is_eqz(a0), complex(eve::inf(eve::as(g))*eve::pedantic(eve::signnz)(real(a0))), f); + return f; + } + else + { + return cayley_extend(tgamma, a0); + } + } + + //===------------------------------------------------------------------------------------------- + // Unary functions : log_gamma + //===------------------------------------------------------------------------------------------- + template + EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& a0) noexcept + { + return log_abs(tgamma(a0)); + } + + template + auto dispatch(eve::tag_of, Z const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // 15 sig. digits for 0<=real(z)<=171 + // coeffs should sum to about g*g/2+23/24 + // + using r_t = eve::element_type_t>; + auto g=r_t(607)/r_t(128); + // best results when 4<=g<=5 + constexpr int N = 15; + std::array c = + { 0.99999999999999709182, + 57.156235665862923517, + -59.597960355475491248, + 14.136097974741747174, + -0.49191381609762019978, + .33994649984811888699e-4, + .46523628927048575665e-4, + -.98374475304879564677e-4, + .15808870322491248884e-3, + -.21026444172410488319e-3, + .21743961811521264320e-3, + -.16431810653676389022e-3, + .84418223983852743293e-4, + -.26190838401581408670e-4, + .36899182659531622704e-5 + }; + + //Num Recipes used g=5 with 7 terms + //for a less effective approximation + + auto negra0 = eve::is_negative(real(a0)); + auto z = if_else(negra0, -a0, a0); + Z ss{}; + for(int pp = N-1; pp >= 1; --pp){ + ss += c[pp]*rec(z+eve::dec(pp)); + } + auto zg = z+g-eve::half(eve::as(g)); + auto lsq2pi = r_t(0.9189385332046727417803297); + auto f=(lsq2pi + log(c[0]+ss)) - zg + (z-eve::half(eve::as()))*log(zg); + auto zer = eve::zero(eve::as(g)); + auto o = eve::one(eve::as(g)); + auto t = o+o; + f = if_else(z == t|| z == o, zer, f); + //adjust for negative real parts + if(eve::any(negra0)) + { + auto lpi = r_t(1.14472988584940017414342735); + auto reala0 = is_real(a0); + f = kyosu::if_else(negra0, lpi-log(a0*sinpi(-a0))-f, f); + f = kyosu::if_else (negra0 && reala0 && eve::is_flint(real(a0)) + , complex(eve::nan(eve::as(g)), eve::inf(eve::as(g))) + , f); + } + return f; + } + else + { + return cayley_extend(log_gamma, a0); + } + } + + //===------------------------------------------------------------------------------------------- + // Unary functions : digamma + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z const& a0) noexcept + { + if constexpr(concepts::complex ) + { + // 15 sig. digits for 0<=real(z)<=171 + // coeffs should sum to about g*g/2+23/24 + // + using r_t = eve::element_type_t>; + auto g=r_t(607)/r_t(128); + // best results when 4<=g<=5 + constexpr int N = 15; + std::array c = + { 0.99999999999999709182, + 57.156235665862923517, + -59.597960355475491248, + 14.136097974741747174, + -0.49191381609762019978, + .33994649984811888699e-4, + .46523628927048575665e-4, + -.98374475304879564677e-4, + .15808870322491248884e-3, + -.21026444172410488319e-3, + .21743961811521264320e-3, + -.16431810653676389022e-3, + .84418223983852743293e-4, + -.26190838401581408670e-4, + .36899182659531622704e-5 + }; + + //Num Recipes used g=5 with 7 terms + //for a less effective approximation + + auto reflection = real(a0) < eve::half(eve::as(real(a0))); + auto z = if_else(reflection, oneminus(a0), a0); + + Z d{}; + auto n = d; + for(int pp = N-1; pp >= 1; --pp){ + auto dz = rec(z+eve::dec(pp)); + auto dd = c[pp]*dz; + d += dd; + n -= dd*dz; + } + d+= c[0]; + auto zg = z+g-eve::half(eve::as(g)); + auto f = log(zg) + (n/d - g/zg); + + if(eve::any(reflection)) + { + f = if_else(reflection, f-eve::pi(eve::as(g))*cotpi(a0), f); + f = if_else (reflection && is_real(a0) && eve::is_flint(real(a0)) + , complex(eve::nan(eve::as(g)), eve::inf(eve::as(g))), f); + } + return f; + } + else + { + return cayley_extend(digamma, a0); + } + } + + //===------------------------------------------------------------------------------------------- + // Binary functions : rising_factorial, lrising_factorial, lbeta, beta + //===------------------------------------------------------------------------------------------- + // WAITING FOR DECORATORS + template + auto dispatch(eve::tag_of + , Z1 const& a0, Z2 const& a1) noexcept + requires (kyosu::concepts::complex || kyosu::concepts::complex) + { + using r_t = as_cayley_dickson_t; + return kyosu::if_else(is_eqz(a1), r_t(1), tgamma(a0+a1)/tgamma(a0)); //pedantic(div)(tgamma(a0+a1),tgamma(a0))); + } + + template + auto dispatch(eve::tag_of + , Z1 const& a0, Z2 const& a1) noexcept + { + using r_t = as_cayley_dickson_t; + return kyosu::if_else(is_eqz(a1), r_t{}, log(tgamma(a0+a1)/tgamma(a0))); //pedantic(divide)(tgamma(a0+a1),tgamma(a0)))); + } + + template + auto dispatch(eve::tag_of, Z1 const& a0, Z2 const& a1) noexcept + { + return log(beta(a0, a1)); + } + + template + auto dispatch(eve::tag_of, Z1 const& a0, Z2 const& a1) noexcept + { + auto y = a0 + a1; + return tgamma(a0)*tgamma(a1)/tgamma(y); + } + + template + auto dispatch(eve::tag_of + , K const & kk + , Z z) noexcept + { + if constexpr(concepts::complex ) + { + // 15 sig. digits for 0<=real(z)<=171 + // coeffs should sum to about g*g/2+23/24 + // + using v_t = kyosu::as_real_t; + using real_t = eve::element_type_t; + auto k = real_t(kk); + auto [rz, iz] = z; + auto iszero = is_eqz(z); + auto isreal = eve::is_eqz(iz); + auto isneg = eve::is_ltz(rz); + auto isnegreal = isreal && isneg; + auto isodd = isnegreal && eve::is_odd(rz); + auto iseven = isnegreal && eve::is_even(rz); + + auto r=eve::half(eve::as(rz)); // reflection point + auto reflect = rz < r; + z = if_else(reflect, oneminus(z), z); + + std::array cm = { + .27105054312137610850e-19 + -.17889335846010823161e-17, + .58167446553847312884e-16, + -.12421162189080181548e-14, + .19593322190397666205e-13, + -.24347803504257137241e-12, + .24823251635643084345e-11, + -.21352608103961806529e-10, + .15816215942184366772e-9, + -.10246226511017621219e-8, + .58768014045093054654e-8, + -.30137695171547022183e-7, + .13931171712321674741e-6, + -.58440580661848562719e-6, + .22376124247437700378e-5, + -.78585149263697370338e-5, + .25423835243950883896e-4, + -.76053287924037718971e-4, + .21106516173760261250e-3, + -.54504190222378945440e-3, + .13131884053420191908e-2, + -.29592166263096543401e-2, + .62512730682449246388e-2, + -.12405987285776082154e-1, + .23176737166455607805e-1, + -.40840766970770029873e-1, + .68016197438946063823e-1, + -.10726959700408922397, + .16054206784249779846, + -.22851039270529494523, + .31007238254065152134, + -.40215850009669926857, + .50000000000000000000, + -.59784149990330073143, + .68992761745934847866, + -.77148960729470505477, + .83945793215750220154, + -.89273040299591077603, + .93198380256105393618, + -.95915923302922997013, + .97682326283354439220, + -.98759401271422391785, + .99374872693175507536, + -.99704078337369034566, + .99868681159465798081, + -.99945495809777621055, + .99978893483826239739, + -.99992394671207596228, + .99997457616475604912, + -.99999214148507363026, + .99999776238757525623, + -.99999941559419338151, + .99999986068828287678, + -.99999996986230482845, + .99999999412319859549, + -.99999999897537734890, + .99999999984183784058, + -.99999999997864739190, + .99999999999751767484, + -.99999999999975652196, + .99999999999998040668, + -.99999999999999875788, + .99999999999999994183, + -.99999999999999999821, + .99999999999999999997 + }; + + Z f{}; + auto j = inc(63*k); + for(size_t i=0; i < cm.size(); ++i, j -= k) + { + f += cm[i]*pow(j, -z); + } + if (eve::none(reflect)) return f; + auto reflection = [&](auto f){ + if (k == 1) + { + auto u = dec(exp2(z)); + auto t = -2*u/dec(u)/pow(eve::pi(eve::as(real(z))), z);; + auto g = t*cospi(eve::half(eve::as(real(z)))*z)*tgamma(z)*f; + g = if_else(iszero, complex(eve::half(eve::as(real(f)))), g); + return if_else(iseven, complex(eve::zero(eve::as(real(f)))), g); + } + else if (k == 2){ + auto g= pow(eve::two_o_pi(eve::as(real(f))), z)*sinpi(z)*tgamma(z)*f; + return if_else(isodd, eve::zero(eve::as(real(f))), g); + } + else + { + return complex(eve::allbits(eve::as(real(f)))); + } + }; + return if_else(reflect, reflection(f), f); + } + else + { + return cayley_extend_rev(deta, kk, z); + } + } + + //===------------------------------------------------------------------------------------------- + // Unary functions : eta + //===------------------------------------------------------------------------------------------- + template + EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept + { + return deta(1u, z); + } + + //===------------------------------------------------------------------------------------------- + // Unary functions : zeta + //===------------------------------------------------------------------------------------------- + template + EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept + { + if constexpr(concepts::complex ) + { + auto zz=exp2(z); + auto k = zz/(zz-2); + auto g = if_else(z == Z(1), complex(eve::nan(eve::as(real(z)))), k*eta(z)); + return if_else(real(z) == eve::inf(eve::as(real(z))), complex(eve::one(eve::as(real(z)))), g); + } + else + { + return cayley_extend(zeta, z); + } + } + + //===------------------------------------------------------------------------------------------- + // Unary functions : lambda + //===------------------------------------------------------------------------------------------- + template + EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept + { + if constexpr(concepts::complex ) + { + auto zz=exp2(z); + auto k = (z-1)/(z-2); + auto r = if_else(z == complex(eve::one(eve::as(real(z)))), complex(eve::inf(eve::as(real(z)))), k*deta(1u, zz)); + imag(r) = eve::if_else(is_real(z), eve::zero, imag(r)); + return r; + } + else + { + return cayley_extend(cosh, z); + } + } +} + +#include +#include + +namespace kyosu::_ +{ + template + EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept + { + if constexpr(concepts::complex ) + { + auto realz = is_real(z); + if (eve::all(realz)) + return complex(eve::erfcx(real(z))); + else if (eve::none(realz)) + return faddeeva(complex(-imag(z), real(z))); + else + return if_else(realz, complex(eve::erfcx(real(z))), faddeeva(complex(-imag(z), real(z)))); + } + else + { + return cayley_extend(erfcx, z); + } + } + + template + EVE_FORCEINLINE auto dispatch(eve::tag_of, Z const& z) noexcept + { + if constexpr(concepts::complex ) + { + auto realz = is_real(z); + if (eve::all(realz)) + { + return kyosu::erfi(real(z)); + } + else + { + auto [rz, iz] = z; + auto tmp = erf(complex(-iz, rz)); + return complex(imag(tmp), -real(tmp)); + } + } + else + { + return cayley_extend(erfi, z); + } + } +} diff --git a/include/kyosu/types/impl/trigo.hpp b/include/kyosu/types/impl/trigo.hpp index 2fba01c5..dd566c56 100644 --- a/include/kyosu/types/impl/trigo.hpp +++ b/include/kyosu/types/impl/trigo.hpp @@ -9,6 +9,7 @@ #include #include +#include namespace kyosu::_ { @@ -34,8 +35,7 @@ namespace kyosu::_ } else { - auto e = kyosu::exp(z); - return kyosu::average(e, rec(e)); + return cayley_extend(cosh, z); } } @@ -73,8 +73,7 @@ namespace kyosu::_ } else { - auto e = kyosu::exp(z); - return eve::average(e, -kyosu::rec(e)); + return cayley_extend(sinh, z); } } @@ -122,9 +121,7 @@ namespace kyosu::_ } else { - auto e = kyosu::exp(z); - auto ie = kyosu::rec(e); - return kumi::tuple{ eve::average(e, -ie), eve::average(e, ie)}; + return cayley_extend2(sinhcosh, z); } } @@ -145,8 +142,7 @@ namespace kyosu::_ } else { - auto e = kyosu::expm1(z+z); - return e/(e+2); + return cayley_extend(tanh, z); } } @@ -160,8 +156,7 @@ namespace kyosu::_ } else { - auto e = kyosu::expm1(z+z); - return (e+2)/e; + return cayley_extend(coth, z); } } @@ -175,11 +170,7 @@ namespace kyosu::_ } else { - auto p = kyosu::pure(z); - auto az = kyosu::abs(p); - auto [s, c] = eve::sincos(real(z)); - auto w = -s*eve::sinhc(az); - return c*cosh(az)+w*pure(z); + return cayley_extend(cos, z); } } @@ -201,11 +192,7 @@ namespace kyosu::_ } else { - auto p = kyosu::pure(z); - auto az = kyosu::abs(p); - auto [s, c] = eve::sincos(real(z)); - auto w = c*eve::sinhc(az); - return s*cosh(az)+w*pure(z); + return cayley_extend(sin, z); } } @@ -228,16 +215,7 @@ namespace kyosu::_ } else { - auto p = kyosu::pure(z); - auto az = kyosu::abs(p); - auto [s, c] = eve::sincos(kyosu::real(z)); - auto shc = eve::sinhc(az); - auto ch = eve::cosh(az); - auto wc = c*shc; - auto ws =-s*shc; - auto sq = s*ch + wc*p; - auto cq = c*ch + ws*p; - return kumi::tuple{sq, cq}; + return cayley_extend2(sincos, z); } } @@ -252,8 +230,7 @@ namespace kyosu::_ } else { - auto [s, c] = sincos(z); - return s/c; + return cayley_extend(tan, z); } } @@ -268,8 +245,7 @@ namespace kyosu::_ } else { - auto [s, c] = sincos(z); - return c/s; + return cayley_extend(cot, z); } } @@ -277,96 +253,15 @@ namespace kyosu::_ KYOSU_FORCEINLINE constexpr auto dispatch(eve::tag_of const&, C const& z) noexcept { - auto s = kyosu::sin(z); - using u_t = eve::underlying_type_t; - return kyosu::if_else(kyosu::abs(z) < eve::eps(eve::as(u_t())), eve::one(eve::as(u_t())), s/z); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return eve::half(eve::as(real(z)))*(log(inc(z))-log(oneminus(z))); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return eve::half(eve::as(real(z)))*(log(inc(z))-log(dec(z))); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return log(z+sqrt(inc(z))*sqrt(dec(z))); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return log(z+sqrt(inc(sqr(z)))); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return acosh(rec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return sinh(rec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - auto s = sign(pure(z)); - return -s*asinh(z*s); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return sign(pure(z))*acosh(z); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - auto s = sign(pure(z)); - return -s*atanh(z*s); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - auto s = sign(pure(z)); - return s*acoth(z*s); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return acos(rec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return asin(rec(z)); + if constexpr(concepts::complex ) + { + auto s = kyosu::sin(z); + using u_t = eve::underlying_type_t; + return kyosu::if_else(kyosu::abs(z) < eve::eps(eve::as(u_t())), eve::one(eve::as(u_t())), s/z); + } + else + { + return cayley_extend(sinc, z); + } } - } diff --git a/include/kyosu/types/impl/trigo_pi.hpp b/include/kyosu/types/impl/trigo_pi.hpp index 473845f8..5d9d3798 100644 --- a/include/kyosu/types/impl/trigo_pi.hpp +++ b/include/kyosu/types/impl/trigo_pi.hpp @@ -36,13 +36,7 @@ namespace kyosu::_ } else { - using u_t = eve::underlying_type_t; - auto pi = eve::pi(eve::as()); - auto p = kyosu::pure(z); - auto az = kyosu::abs(p)*eve::pi(eve::as(kyosu::abs(p))); - auto [s, c] = eve::sinpicospi(real(z)); - auto w = -s*eve::sinhc(az)*pi; - return c*cosh(az)+w*pure(z); + return cayley_extend(cospi, z); } } @@ -79,13 +73,7 @@ namespace kyosu::_ } else { - using u_t = eve::underlying_type_t; - auto pi = eve::pi(eve::as()); - auto p = kyosu::pure(z); - auto az = kyosu::abs(p)*pi; - auto [s, c] = eve::sinpicospi(kyosu::real(z)); - auto w = c*eve::sinhc(az)*pi; - return s*cosh(az)+w*p; + return cayley_extend(sinpi, z); } } @@ -137,18 +125,7 @@ namespace kyosu::_ } else { - using u_t = eve::underlying_type_t; - auto pi = eve::pi(eve::as()); - auto p = kyosu::pure(z); - auto az = kyosu::abs(p)*pi; - auto [s, c] = eve::sinpicospi(kyosu::real(z)); - auto shc = eve::sinhc(az)*pi; - auto ch = eve::cosh(az); - auto wc = c*shc; - auto ws =-s*shc; - auto sq = s*ch + wc*p; - auto cq = c*ch + ws*p; - return kumi::tuple{sq, cq}; + return cayley_extend2(sinpicospi, z); } } @@ -172,8 +149,7 @@ namespace kyosu::_ } else { - auto [s, c] = kyosu::sinpicospi(z); - return s/c; + return cayley_extend(tanpi, z); } } @@ -190,51 +166,7 @@ namespace kyosu::_ } else { - auto [s, c] = kyosu::sinpicospi(z); - return c/s; + return cayley_extend(cotpi, z); } } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return radinpi(asin(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return radinpi(acos(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return radinpi(atan(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return radinpi(acot(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return radinpi(asec(z)); - } - - template - KYOSU_FORCEINLINE constexpr - auto dispatch(eve::tag_of const&, C const& z) noexcept - { - return radinpi(acsc(z)); - } - }