diff --git a/CMakeLists.txt b/CMakeLists.txt index da75d79..1080005 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,23 +9,23 @@ target_compile_features(twofloat INTERFACE cxx_std_17) add_subdirectory(test) -option(BUILD_DOC "Build documentation" ON) -find_package(Doxygen) -if (DOXYGEN_FOUND) - # set input and output files - set(DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/docs/Doxyfile.in) - set(DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile) - - # request to configure the file - configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY) - message("Doxygen build started") - - # note the option ALL which allows to build the docs together with the application - add_custom_target( doc_doxygen ALL - COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT} - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMENT "Generating API documentation with Doxygen" - VERBATIM ) -else (DOXYGEN_FOUND) - message("Doxygen need to be installed to generate the doxygen documentation") -endif (DOXYGEN_FOUND) \ No newline at end of file +#option(BUILD_DOC "Build documentation" ON) +#find_package(Doxygen) +#if (DOXYGEN_FOUND) +# # set input and output files +# set(DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/docs/Doxyfile.in) +# set(DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile) +# +# # request to configure the file +# configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY) +# message("Doxygen build started") +# +# # note the option ALL which allows to build the docs together with the application +# add_custom_target( doc_doxygen ALL +# COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT} +# WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +# COMMENT "Generating API documentation with Doxygen" +# VERBATIM ) +#else (DOXYGEN_FOUND) +# message("Doxygen need to be installed to generate the doxygen documentation") +#endif (DOXYGEN_FOUND) \ No newline at end of file diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index e42202b..4aa4e62 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -233,5 +233,423 @@ inline two div(const two &x, const two &y) { } else static_assert(sizeof(T) == 0, "Unsupported mode"); } + +/// \brief QD constants +/// Reference: QD / dd_const.cpp / inline.h / dd_real.cpp +template +static const two _2pi = {6.283185307179586232e+00, 2.449293598294706414e-16}; + +template +static const two _pi2 = {1.570796326794896558e+00, 6.123233995736766036e-17}; + +template +static const two _pi16 = {1.963495408493620697e-01, 7.654042494670957545e-18}; + +template +static const constexpr T _nan = std::numeric_limits::quiet_NaN(); + +template +static const constexpr T _eps = 4.93038065763132e-32; // 2^-104 + +static const int n_inv_fact = 15; + +/// \brief Some recurrently-used contants +template +static const constexpr T pointfive = 0.5; + +template +static const constexpr T twopointzero = 2.0; + +template +static const constexpr T zeropointzero = 0.0; + +template +static const constexpr T onepointzero = 1.0; + +/// \brief Tables containing constant data for trigonometric functions +template +struct constants_trig_tables { + + static const constexpr T inv_fact[n_inv_fact][2] = { + { 1.66666666666666657e-01, 9.25185853854297066e-18}, + { 4.16666666666666644e-02, 2.31296463463574266e-18}, + { 8.33333333333333322e-03, 1.15648231731787138e-19}, + { 1.38888888888888894e-03, -5.30054395437357706e-20}, + { 1.98412698412698413e-04, 1.72095582934207053e-22}, + { 2.48015873015873016e-05, 2.15119478667758816e-23}, + { 2.75573192239858925e-06, -1.85839327404647208e-22}, + { 2.75573192239858883e-07, 2.37677146222502973e-23}, + { 2.50521083854417202e-08, -1.44881407093591197e-24}, + { 2.08767569878681002e-09, -1.20734505911325997e-25}, + { 1.60590438368216133e-10, 1.25852945887520981e-26}, + { 1.14707455977297245e-11, 2.06555127528307454e-28}, + { 7.64716373181981641e-13, 7.03872877733453001e-30}, + { 4.77947733238738525e-14, 4.39920548583408126e-31}, + { 2.81145725434552060e-15, 1.65088427308614326e-31} + }; + + // Reference: QD / dd_real.cpp + // Table of sin(k * pi/16) and cos(k * pi/16). + static const constexpr T cos_table [4][2] = { + {9.807852804032304306e-01, 1.854693999782500573e-17}, + {9.238795325112867385e-01, 1.764504708433667706e-17}, + {8.314696123025452357e-01, 1.407385698472802389e-18}, + {7.071067811865475727e-01, -4.833646656726456726e-17} + }; + static const constexpr T sin_table [4][2] = { + {1.950903220161282758e-01, -7.991079068461731263e-18}, + {3.826834323650897818e-01, -1.005077269646158761e-17}, + {5.555702330196021776e-01, 4.709410940561676821e-17}, + {7.071067811865475727e-01, -4.833646656726456726e-17} + }; +}; + +/// \brief Computes fl(input * input) and err(input * input) +/// Reference: QD / inline.h (within qd namespace) +template +inline T two_sqr(T input, T &err) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + T hi, lo; + T q = input * input; + two temp = twofloat::algorithms::Split(input); + hi = temp.h; + lo = temp.l; + err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; + return q; +} + +/// Reference QD / inline.h (within qd namespace) +/// TODO: make sure this has not been already implemented in twofloat +template +inline T two_sum(T a, T b, T &err) { + T s = a + b; + T bb = s - a; + err = (a - (s - bb)) - (b - bb); + return s; +} + +/// Reference: QD / dd_inline.h (within dd_real namespace) +/// TODO: make sure this has not been already implemented in twofloat +template +inline two add(T a, T b) { + T s, e; + s = two_sum(a, b, e); + return two (s, e); +} + +/// \brief Computes the nearest integer to input. +/// Reference: QD / inline.h (within qd namespace) +template +inline T nint(T input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + if(input = std::floor(input)) { + return input; + } + return std::floor(input + pointfive); +} + +/// \brief Round to nearest integer. +/// Reference: QD / dd_inline.h +template +inline two nint(const two &input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + T hi = nint(input.h); + T lo; + if(hi == input.h) { + // High word is an integer already. Round the low word + lo = nint(input.l); + + // Renormalize. This is needed if h = some integer, l = 1/2 + //hi = qd::quick_two_sum(hi, lo, lo); + two temp = twofloat::algorithms::FastTwoSum(hi, lo); + hi = temp.h; + lo = temp.l; + } + else { + // High word is not an integer + lo = zeropointzero; + if(std::abs(hi - input.h) == pointfive && input.l < zeropointzero) { + // There is a tie in the high word, consult the low word to break the tie + hi -= onepointzero; /* NOTE: This does not cause INEXACT. */ + } + } + + return two{hi, lo}; +} + +/// Reference: QD / dd_inline.h (within dd_real namespace) +template +inline two sqr(T input) { + T p1, p2; + p1 = two_sqr(input, p2); + return two{p1, p2}; +} + +/// Reference: QD / dd_inline.h +template +inline two sqr(const two &input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + T p1, p2; + T s1, s2; + p1 = two_sqr(input.h, p2); + p2 += twopointzero * input.h * input.l; + p2 += input.l * input.l; + two temp = twofloat::algorithms::FastTwoSum(p1, p2); + s1 = temp.h; + s2 = temp.l; + return two{s1, s2}; +} + +/// \brief Computes sin(a) using Taylor series. +/// Assumes |a| <= pi/32 +/// Reference: QD / dd_real.cpp +template +static two sin_taylor(const two &input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + + const T* local_ptr_inv_fact = &constants_trig_tables::inv_fact[0][0]; + T local_eps = _eps; + const T thresh = pointfive * std::abs(input.eval()) * local_eps; + two r, s, t, x; + + if (input.eval() == zeropointzero) { + return two{zeropointzero, zeropointzero}; + } + + int i = 0; + two temp = sqr(input); + x = two{-temp.h, -temp.l}; + s = input; + r = input; + + do { + r = mul(r, x); + t = mul(r, two{(local_ptr_inv_fact+i)[0], (local_ptr_inv_fact+i)[1]}); + s = add

(s, t); + i += 2; + } while (i < n_inv_fact && std::abs(t.eval()) > thresh); + + return s; +} + +/// \brief double-double * double, where double is a power of 2. +/// Reference: QD / dd_inline.h +template +inline two mul_pwr2(const two &input, T b) { + return two(input.h * b, input.l * b); +} + +/// Reference: QD / dd_real.cpp +template +static two cos_taylor(const two &input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + + const T* local_ptr_inv_fact = &constants_trig_tables::inv_fact[0][0]; + T local_eps = _eps; + const T thresh = pointfive * local_eps; + two r, s, t, x; + + if (input.eval() == zeropointzero) { + return two{onepointzero, zeropointzero}; + } + + two temp = sqr(input); + x = two{-temp.h, -temp.l}; + r = x; + s = add(onepointzero, mul_pwr2(r, pointfive)); + int i = 1; + do { + r = mul(r, x); + t = mul(r, two{(local_ptr_inv_fact+i)[0], (local_ptr_inv_fact+i)[1]}); + s = add

(s, t); + i += 2; + } while (i < n_inv_fact && std::abs(t.eval()) > thresh); + + return s; +} + +/// \brief Computes the square root of the double-double number dd. +/// NOTE: dd must be a non-negative number. +/// +/// Strategy: Use Karp's trick: if x is an approximation +/// to sqrt(a), then +/// +/// sqrt(a) = a*x + [a - (a*x)^2] * x / 2 (approx) +/// +/// The approximation is accurate to twice the accuracy of x. +/// Also, the multiplication (a*x) and [-]*x can be done with +/// only half the precision. +/// +/// Reference: QD / dd_real.cpp +template +two sqrt(const two &input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + T local_nan = _nan; + + if (input.eval() == zeropointzero) { + return two{zeropointzero, zeropointzero}; + } + + // ERROR: Negative argument + // if (input.is_negative()) + if (input.h < zeropointzero) { + return two{local_nan, local_nan}; + } + + T x = onepointzero / std::sqrt(input.h); + T ax = input.h * x; + two temp = sub

(input, sqr(ax)); + return add(ax, temp.h * x * pointfive); +} + +/// Reference: QD / dd_real.cpp +template +static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + + if (input.eval() == zeropointzero) { + sin_a.h = zeropointzero; + sin_a.l = zeropointzero; + cos_a.h = onepointzero; + cos_a.l = zeropointzero; + return; + } + + sin_a = sin_taylor(input); + cos_a = sqrt

(sub(onepointzero, sqr(sin_a))); +} + +/// Reference: QD / dd_real.cpp +template +inline two sin(const two &input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + + two local_2pi = _2pi; + two local_pi2 = _pi2; + two local_pi16 = _pi16; + T local_nan = _nan; + + const T* local_ptr_cos_table = &constants_trig_tables::cos_table[0][0]; + const T* local_ptr_sin_table = &constants_trig_tables::sin_table[0][0]; + + if (input.eval() == zeropointzero) { + return two{zeropointzero, zeropointzero}; + } + + // Approximately reducing modulo 2*pi + two z = nint(div(input, local_2pi)); + two r = sub

(input, mul(local_2pi, z)); + + /* + std::cout << "LSV starts ..." << std::endl; + std::cout << "input = " << input.eval() << std::endl; + std::cout << "z = " << z.eval() << std::endl; + std::cout << "r = " << r.eval() << std::endl; + std::cout << "... LSV ends" << std::endl; + */ + + // Approximately reducing modulo pi/2 and then modulo pi/16 + T q = std::floor(r.h / local_pi2.h + pointfive); + two t = sub

(r, mul(local_pi2, q)); + int j = static_cast(q); + q = std::floor(t.h / local_pi16.h + pointfive); + t = sub

(t, mul(local_pi16, q)); + int k = static_cast(q); + int abs_k = std::abs(k); + + /* + std::cout << "LSV starts ..." << std::endl; + std::cout << "t = " << t.eval() << std::endl; + std::cout << "q = " << q << std::endl; + std::cout << "j = " << j << std::endl; + std::cout << "k = " << k << std::endl; + std::cout << "abs_k = " << abs_k << std::endl; + std::cout << "... LSV ends" << std::endl; + */ + + // ERROR: cannot reduce modulo pi/2\n + if (j < -2 || j > 2) { + return two{local_nan, local_nan}; + } + + // ERROR: cannot reduce modulo pi/16 + if (abs_k > 4) { + return two{local_nan, local_nan}; + } + + if (k == 0) { + switch(j) { + case 0: + return sin_taylor(t); + case 1: + return cos_taylor(t); + case -1: + return two{-cos_taylor(t).h, -cos_taylor(t).l}; // Reference: QD / dd_inline.h + default: + return two{-sin_taylor(t).h, -sin_taylor(t).l}; + } + } + + // IMPORTANT: notice the factor 2 when accessing the cos() and sin() look-up tables + // This factor is required to access data correctly as done in QD + two u{(local_ptr_cos_table + 2 * (abs_k - 1))[0], (local_ptr_cos_table + 2 * (abs_k - 1))[1]}; + two v{(local_ptr_sin_table + 2 * (abs_k - 1))[0], (local_ptr_sin_table + 2 * (abs_k - 1))[1]}; + + two sin_t, cos_t; + + sincos_taylor(t, sin_t, cos_t); + + /* + std::cout << "LSV starts ..." << std::endl; + std::cout << "sin_t = " << sin_t.eval() << std::endl; + std::cout << "cos_t = " << cos_t.eval() << std::endl; + std::cout << "... LSV ends" << std::endl; + */ + + if (j == 0) { + if (k > 0) { + r = add

(mul(u, sin_t), mul(v, cos_t)); + } else { + r = sub

(mul(u, sin_t), mul(v, cos_t)); + } + } else if (j == 1) { + if (k > 0) { + r = sub

(mul(u, cos_t), mul(v, sin_t)); + } else { + r = add

(mul(u, cos_t), mul(v, sin_t)); + /* + std::cout << "LSV starts ..." << std::endl; + std::cout << "u = " << u.eval() << std::endl; + std::cout << "cos_t = " << cos_t.eval() << std::endl; + std::cout << "v = " << v.eval() << std::endl; + std::cout << "sin_t = " << sin_t.eval() << std::endl; + std::cout << "r = " << r.eval() << std::endl; + std::cout << "... LSV ends" << std::endl; + */ + } + } else if (j == -1) { + if (k > 0) { + r = sub

(mul(v, sin_t), mul(u, cos_t)); + } else if (k < 0) { + two temp = mul(u, cos_t); + two neg_temp = two{-temp.h, -temp.l}; + r = sub

(neg_temp, mul(v, sin_t)); + } + } else { + if (k > 0) { + two temp = mul(u, sin_t); + two neg_temp = two{-temp.h, -temp.l}; + r = sub

(neg_temp, mul(v, cos_t)); + } else { + r = sub

(mul(v, cos_t), mul(u, sin_t)); + } + } + /* + std::cout << "LSV starts ..." << std::endl; + std::cout << "r = " << r.eval() << std::endl; + std::cout << "... LSV ends" << std::endl; + */ + return r; +} } // namespace doubleword } // namespace twofloat \ No newline at end of file diff --git a/include/libtwofloat/twofloat.hpp b/include/libtwofloat/twofloat.hpp index 9134072..18e7403 100644 --- a/include/libtwofloat/twofloat.hpp +++ b/include/libtwofloat/twofloat.hpp @@ -6,6 +6,7 @@ #include #include #include +#include /// \brief A namespace for algorithms that operate on pairs of floating point namespace twofloat { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 522bb9d..0ed77ec 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -7,7 +7,8 @@ FetchContent_Declare( ) FetchContent_MakeAvailable(googletest) -add_executable(twofloat_test pair-arithmetic.test.cpp double-word-arithmetic.test.cpp) +#add_executable(twofloat_test pair-arithmetic.test.cpp double-word-arithmetic.test.cpp) +add_executable(twofloat_test double-word-arithmetic.test.cpp) target_link_libraries(twofloat_test twofloat GTest::gtest_main) include(GoogleTest) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index d85614b..ed0b690 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -255,6 +255,30 @@ TEST(DoubleWordArithmetic, DivFPTest) { divFPTest(); } TEST(DoubleWordArithmetic, DivFPFMATest) { divFPTest(); } +TEST(DoubleWordArithmetic, SinTest) { + float x = 1; + float sinx = std::sin(x); + + two two_x = {x, 0.0f}; + two sin_two_x = doubleword::sin(two_x); + + // Check the expected and actual values within /build/Testing/Temporary/LastTest.log + EXPECT_FLOAT_EQ(sinx, sin_two_x.eval()); + //EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); +} + +TEST(DoubleWordArithmetic, SinTest2) { + double x = 1; + double sinx = std::sin(x); + + two two_x = {x, 0.0}; + two sin_two_x = doubleword::sin(two_x); + + // Check the expected and actual values within /build/Testing/Temporary/LastTest.log + EXPECT_DOUBLE_EQ(sinx, sin_two_x.eval()); + //EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); +} + } // namespace test } // namespace doubleword } // namespace twofloat \ No newline at end of file