Skip to content

Commit

Permalink
Implement and test asin()
Browse files Browse the repository at this point in the history
  • Loading branch information
@ckormanyos
ckormanyos committed Feb 5, 2024
1 parent fe6afb1 commit a429de8
Show file tree
Hide file tree
Showing 2 changed files with 166 additions and 8 deletions.
110 changes: 103 additions & 7 deletions math/softfloat/soft_double.h
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
///////////////////////////////////////////////////////////////////
// Copyright Christopher Kormanyos 2012 - 2023. //
// Copyright Christopher Kormanyos 2012 - 2024. //
// Distributed under the Boost Software License, //
// Version 1.0. (See accompanying file LICENSE_1_0.txt //
// or copy at http://www.boost.org/LICENSE_1_0.txt) //
Expand Down Expand Up @@ -592,6 +592,9 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
constexpr auto sin (soft_double x) -> soft_double;
constexpr auto cos (soft_double x) -> soft_double;
constexpr auto tan (soft_double x) -> soft_double;
constexpr auto asin (soft_double x) -> soft_double;
constexpr auto acos (soft_double x) -> soft_double;
constexpr auto atan (soft_double x) -> soft_double;
constexpr auto sinh (soft_double x) -> soft_double;
constexpr auto cosh (soft_double x) -> soft_double;
constexpr auto tanh (soft_double x) -> soft_double;
Expand All @@ -614,7 +617,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

using representation_type = std::uint64_t;

constexpr soft_double() noexcept = default;
constexpr soft_double() noexcept : my_value(static_cast<std::uint64_t>(UINT8_C(0))) { }

template<typename UnsignedIntegralType,
typename std::enable_if<( std::is_integral<UnsignedIntegralType>::value
Expand Down Expand Up @@ -2409,6 +2412,57 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
return soft_double::my_value_one() + ((x2 * top) / (bot * static_cast<int>(INT8_C(24))));
}

constexpr auto asin_pade(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param)
{
// Usa a Pade approximation of (asin(x) / x).
// Then simplify, extract the coefficients of numerator and denominator.
// Scale both numeroator and denominator by a factor of 10^29.
// Subsequently extract the integer parts of the coefficients.

// PadeApproximant[ArcSin[x]/x, {x, 0, {12, 12}}]
// FullSimplify[%]
// IntegerPart[N[CoefficientList[Numerator[Out[2]]/10^29, x^2], 24]]
// IntegerPart[N[CoefficientList[Denominator[Out[2]]/10^29, x^2], 24]]

// Numerator:
// 140095773199074572
// -388690947968573359
// 402271633106633823
// -190922593544635779
// 40839234741969911
// -3218334247465525
// 45053024618672

// Denominator:
// 140095773199074572
// -412040243501752455
// 460437824033661973
// -243013488210192408
// 60946917703998365
// -6321063389564933
// 174545740275468

const auto x2 = x * x;

const soft_double top = ((((((( + INT64_C(45053024618672))
* x2 - INT64_C(3218334247465525))
* x2 + INT64_C(40839234741969911))
* x2 - INT64_C(190922593544635779))
* x2 + INT64_C(402271633106633823))
* x2 - INT64_C(388690947968573359))
* x2 + INT64_C(140095773199074572));

const soft_double bot = ((((((( + INT64_C(174545740275468))
* x2 - INT64_C(6321063389564933))
* x2 + INT64_C(60946917703998365))
* x2 - INT64_C(243013488210192408))
* x2 + INT64_C(460437824033661973))
* x2 - INT64_C(412040243501752455))
* x2 + INT64_C(140095773199074572));

return (x * top) / bot;
}

} // namespace detail

constexpr auto sin(soft_double x) -> soft_double // NOLINT(misc-no-recursion)
Expand Down Expand Up @@ -2607,6 +2661,53 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
return c;
}

constexpr auto tan(soft_double x) -> soft_double
{
return sin(x) / cos(x);
}

constexpr auto asin(soft_double x) -> soft_double // NOLINT(misc-no-recursion)
{
soft_double result { };

if(x < static_cast<int>(INT8_C(0)))
{
result = -asin(-x);
}
else if(x > static_cast<int>(INT8_C(0)))
{
if(x > static_cast<int>(INT8_C(1)))
{
result = soft_double::my_value_quiet_NaN();
}
else if(x < static_cast<int>(INT8_C(1)))
{
if(x < soft_double::my_value_epsilon())
{
result = x;
}
else if(x < soft_double::my_value_half())
{
result = detail::asin_pade(x);
}
else
{
result = soft_double::my_value_pi_half() - (detail::asin_pade(sqrt((soft_double::my_value_one() - x) / 2)) * 2);
}
}
else
{
result = soft_double::my_value_pi_half();
}
}
else
{
// x == 0 and result = 0.
}

return result;
}

constexpr auto floor(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param)
{
auto result = soft_double { };
Expand Down Expand Up @@ -2781,11 +2882,6 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
return exp(a * log(x)); // NOLINT(performance-unnecessary-value-param)
}

constexpr auto tan(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param)
{
return sin(x) / cos(x);
}

constexpr auto sinh(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param)
{
const soft_double ep = exp(x); // NOLINT(performance-unnecessary-value-param)
Expand Down
64 changes: 63 additions & 1 deletion test/test_soft_double.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
///////////////////////////////////////////////////////////////////
// Copyright Christopher Kormanyos 2012 - 2023. //
// Copyright Christopher Kormanyos 2012 - 2024. //
// Distributed under the Boost Software License, //
// Version 1.0. (See accompanying file LICENSE_1_0.txt //
// or copy at http://www.boost.org/LICENSE_1_0.txt) //
Expand Down Expand Up @@ -454,6 +454,66 @@ auto test_cos() -> bool
return result_is_ok;
}

auto test_asin() -> bool
{
bool result_is_ok = true;

for(int i = 0; i < 50; ++i) // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
{
const math::softfloat::float64_t x = math::softfloat::float64_t(i) / 100; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
const auto d = static_cast<double>(x);

using std::asin;

const math::softfloat::float64_t as_x = asin(x);
const double as_d = asin(d);

if(i == 0)
{
const auto result_asin_zero_is_ok = (as_x.crepresentation() == 0U);

result_is_ok = (result_asin_zero_is_ok && result_is_ok);
}
else
{
const double closeness = std::fabs(1.0 - (static_cast<double>(as_x) / as_d));

const auto result_asin_is_ok = (closeness < std::numeric_limits<double>::epsilon() * 32.0); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)

result_is_ok = (result_asin_is_ok && result_is_ok);
}
}

{
const math::softfloat::float64_t x = math::softfloat::float64_t(-33) / 100;
const auto d = static_cast<double>(x);

using std::asin;

const math::softfloat::float64_t as_x = asin(x);
const double as_d = asin(d);

const double closeness = std::fabs(1.0 - (static_cast<double>(as_x) / as_d));

const auto result_asin_is_ok = ((as_x < 0) && (closeness < std::numeric_limits<double>::epsilon() * 32.0)); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)

result_is_ok = (result_asin_is_ok && result_is_ok);
}

for(int i = 121; i < 124; ++i) // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
{
const math::softfloat::float64_t x = math::softfloat::float64_t(i) / 100; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)

const math::softfloat::float64_t as_x = asin(x);

const auto result_asin_is_ok = isnan(as_x);

result_is_ok = (result_asin_is_ok && result_is_ok);
}

return result_is_ok;
}

auto test_floor(const std::uint32_t n) -> bool
{
bool result_is_ok = true;
Expand Down Expand Up @@ -751,6 +811,7 @@ auto test_soft_double() -> bool
std::cout << "testing log____... "; const auto result_log____is_ok = test_log ( ); std::cout << std::boolalpha << result_log____is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp)
std::cout << "testing sin____... "; const auto result_sin____is_ok = test_sin ( ); std::cout << std::boolalpha << result_sin____is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp)
std::cout << "testing cos____... "; const auto result_cos____is_ok = test_cos ( ); std::cout << std::boolalpha << result_cos____is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp)
std::cout << "testing asin___... "; const auto result_asin___is_ok = test_asin ( ); std::cout << std::boolalpha << result_asin___is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp)

local::eng32.seed(::util::util_pseudorandom_time_point_seed::value<typename local::eng32_type::result_type>());
local::eng64.seed(::util::util_pseudorandom_time_point_seed::value<typename local::eng64_type::result_type>());
Expand Down Expand Up @@ -798,6 +859,7 @@ auto test_soft_double() -> bool
&& result_log____is_ok
&& result_sin____is_ok
&& result_cos____is_ok
&& result_asin___is_ok
&& result_floor__is_ok
&& result_ceil___is_ok
&& result_add____is_ok
Expand Down

0 comments on commit a429de8

Please sign in to comment.