From decdb4bd281421973cc4b57ee742dc5c53b39a1d Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 09:44:45 -0500 Subject: [PATCH 001/137] Commenting out Doxygen generation in CMakeLists.txt. --- CMakeLists.txt | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) 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 From 5cee9647e5668fc3a6970a2c078c6b58b5280414 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 09:46:16 -0500 Subject: [PATCH 002/137] Starting to implement two sin(). --- .../arithmetics/double-word-arithmetic.hpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index e42202b..1e2db12 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -233,5 +233,24 @@ inline two div(const two &x, const two &y) { } else static_assert(sizeof(T) == 0, "Unsupported mode"); } + +// TODO: transform this to FP32 (this is FP64) +const two fp32_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); + +// TODO: check it was correctly copied from QD +const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); + +template +inline two sin(const two &input) { + + // TODO: make sure this is already supported in twofloat + if (input.is_zero()) { + return 0.0; + } + + // Approximately reducing modulo 2*pi + //two z = + +} } // namespace doubleword } // namespace twofloat \ No newline at end of file From 08aa7499bbd7caac35621f462650597174c6cf6f Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 09:56:13 -0500 Subject: [PATCH 003/137] For single-precision: reduce the number of precised decimal digits down to 7. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 1e2db12..7fe97d6 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -235,7 +235,7 @@ inline two div(const two &x, const two &y) { } // TODO: transform this to FP32 (this is FP64) -const two fp32_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); +const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); // TODO: check it was correctly copied from QD const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); From 56f722f2c01e9cfee4bf21c9da116eb9e05230f3 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 10:15:10 -0500 Subject: [PATCH 004/137] Adding implementation of nint(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7fe97d6..198e5a2 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -234,12 +234,22 @@ inline two div(const two &x, const two &y) { static_assert(sizeof(T) == 0, "Unsupported mode"); } +// Reference: QD / dd_const.cpp // TODO: transform this to FP32 (this is FP64) const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); // TODO: check it was correctly copied from QD const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); +// Reference: QD / inline.h +/* Computes the nearest integer to input. */ +inline float nint(float input) { + if(input = std::floor(input)) { + return input; + } + return std::floor(input + 0.5f); +} + template inline two sin(const two &input) { From 495a23aa94973d1ddd739c5df4a5ba85e59a6bb4 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 10:15:59 -0500 Subject: [PATCH 005/137] Adding reference implementation. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 198e5a2..7fdb2d8 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -250,6 +250,7 @@ inline float nint(float input) { return std::floor(input + 0.5f); } +// Reference: QD / dd_real.cpp template inline two sin(const two &input) { From a976da80e993c2c2763a06e47916c44b165e2e18 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 10:27:52 -0500 Subject: [PATCH 006/137] Parameterizing nint()'s implementation + extra statements. --- .../arithmetics/double-word-arithmetic.hpp | 27 ++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7fdb2d8..3340369 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -243,11 +243,22 @@ const two fp64_2pi = two(6.283185307179586232e+00, 2.44929359829 // Reference: QD / inline.h /* Computes the nearest integer to input. */ -inline float nint(float input) { +template +inline T nint(T input) { + T pointfive; + + if constexpr (std::is_same_v) { + pointfive = 0.5f; + } else if constexpr (std::is_same_v) { + pointfive = 0.5; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + if(input = std::floor(input)) { return input; } - return std::floor(input + 0.5f); + return std::floor(input + pointfive); } // Reference: QD / dd_real.cpp @@ -259,8 +270,18 @@ inline two sin(const two &input) { return 0.0; } + T local_2pi; + if constexpr (std::is_same_v) { + local_2pi = fp32_2pi; + } else if constexpr (std::is_same_v) { + local_2pi = fp64_2pi; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + // Approximately reducing modulo 2*pi - //two z = + two z = nint(input / local_2pi); + two r = input - mul(local_2pi, z); } } // namespace doubleword From 4ebf10499e8512db37691074d72b42c24c7e31ac Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 10:50:07 -0500 Subject: [PATCH 007/137] Adding definition of . --- .../arithmetics/double-word-arithmetic.hpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 3340369..f98fa46 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -237,9 +237,11 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp // TODO: transform this to FP32 (this is FP64) const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); +const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); // TODO: check it was correctly copied from QD const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); +const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); // Reference: QD / inline.h /* Computes the nearest integer to input. */ @@ -270,11 +272,15 @@ inline two sin(const two &input) { return 0.0; } - T local_2pi; + T local_2pi, local_pi2, pointfive; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; + local_pi2 = fp32_pi2; + pointfive = 0.5f; } else if constexpr (std::is_same_v) { local_2pi = fp64_2pi; + local_pi2 = fp64_pi2; + pointfive = 0.5; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } @@ -283,6 +289,12 @@ inline two sin(const two &input) { two z = nint(input / local_2pi); two r = input - mul(local_2pi, z); + // Approximately reducing modulo pi/2 and then modulo pi/16 + two t; + + //TODO: original type is double, here it is templated + T q = std::floor(r.h / local_pi2.h + pointfive); + } } // namespace doubleword } // namespace twofloat \ No newline at end of file From 97ca2b90c7fb0f5528ea70fc5f85924da767fc51 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 10:55:44 -0500 Subject: [PATCH 008/137] Adding definition of . --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index f98fa46..312ee40 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -290,10 +290,10 @@ inline two sin(const two &input) { two r = input - mul(local_2pi, z); // Approximately reducing modulo pi/2 and then modulo pi/16 - two t; //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); + two t = r - mul(local_pi2, q) } } // namespace doubleword From 718f62d890313528ee63684cb6ad821dbce5aafd Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 11:03:37 -0500 Subject: [PATCH 009/137] Adding missing ; --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 312ee40..69d3440 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -293,7 +293,7 @@ inline two sin(const two &input) { //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); - two t = r - mul(local_pi2, q) + two t = r - mul(local_pi2, q); } } // namespace doubleword From ec988141ca5977a20f8b148a8cf0ccc7a90117c1 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 11:12:13 -0500 Subject: [PATCH 010/137] Recalculation of . --- .../arithmetics/double-word-arithmetic.hpp | 28 +++++++++++-------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 69d3440..d883b20 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -236,12 +236,14 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp // TODO: transform this to FP32 (this is FP64) -const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); -const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); +const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); +const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); +const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); // TODO: check it was correctly copied from QD -const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); -const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); +const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); +const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); +const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); // Reference: QD / inline.h /* Computes the nearest integer to input. */ @@ -272,15 +274,17 @@ inline two sin(const two &input) { return 0.0; } - T local_2pi, local_pi2, pointfive; + T local_2pi, local_pi2, pointfive, local_pi16; if constexpr (std::is_same_v) { - local_2pi = fp32_2pi; - local_pi2 = fp32_pi2; - pointfive = 0.5f; + local_2pi = fp32_2pi; + local_pi2 = fp32_pi2; + pointfive = 0.5f; + local_pi16 = fp32_pi16: } else if constexpr (std::is_same_v) { - local_2pi = fp64_2pi; - local_pi2 = fp64_pi2; - pointfive = 0.5; + local_2pi = fp64_2pi; + local_pi2 = fp64_pi2; + pointfive = 0.5; + local_pi16 = fp64_pi16: } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } @@ -294,6 +298,8 @@ inline two sin(const two &input) { //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); two t = r - mul(local_pi2, q); + int j = static_cast(q); + q = std::floor(t.h / local_pi16.h + pointfive); } } // namespace doubleword From fc9fd58cbccf53388252a3a995a91422c9c55346 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 11:15:00 -0500 Subject: [PATCH 011/137] Definitions of and . --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index d883b20..de609b0 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -300,6 +300,9 @@ inline two sin(const two &input) { two t = r - mul(local_pi2, q); int j = static_cast(q); q = std::floor(t.h / local_pi16.h + pointfive); + t -= mul(local_pi16, q); + int k = static_cast(q); + int abs_k = std:abs(k); } } // namespace doubleword From 06b7948ae590238a7864c5b23411ff031343cd88 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 11:36:41 -0500 Subject: [PATCH 012/137] Adding conditional possibly returning nan. --- .../arithmetics/double-word-arithmetic.hpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index de609b0..f9d87d3 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -234,16 +234,18 @@ inline two div(const two &x, const two &y) { static_assert(sizeof(T) == 0, "Unsupported mode"); } -// Reference: QD / dd_const.cpp +// Reference: QD / dd_const.cpp / inline.h // TODO: transform this to FP32 (this is FP64) const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); +static const float fp32_nan = std::numeric_limits::quiet_NaN(); // TODO: check it was correctly copied from QD const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); +static const double fp64_nan = std::numeric_limits::quiet_NaN(); // Reference: QD / inline.h /* Computes the nearest integer to input. */ @@ -274,17 +276,19 @@ inline two sin(const two &input) { return 0.0; } - T local_2pi, local_pi2, pointfive, local_pi16; + T local_2pi, local_pi2, pointfive, local_pi16, local_nan; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; local_pi2 = fp32_pi2; pointfive = 0.5f; local_pi16 = fp32_pi16: + local_nan = fp32_nan; } else if constexpr (std::is_same_v) { local_2pi = fp64_2pi; local_pi2 = fp64_pi2; pointfive = 0.5; - local_pi16 = fp64_pi16: + local_pi16 = fp64_pi16; + local_nan = fp64_nan; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } @@ -304,6 +308,11 @@ inline two sin(const two &input) { int k = static_cast(q); int abs_k = std:abs(k); + if (j < -2 || j > 2) { + std::error("LSV: cannot reduce modulo pi/2"); + return local_nan; + } + } } // namespace doubleword } // namespace twofloat \ No newline at end of file From 6f0061b5087738517f5928fe4db9e4fe5ba67d5b Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 11:39:55 -0500 Subject: [PATCH 013/137] Adding conditional (related to pi/16) possibly returning nan. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index f9d87d3..81a2e23 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -313,6 +313,11 @@ inline two sin(const two &input) { return local_nan; } + if (abs_k > 4) { + std::error("LSV: cannot reduce modulo pi/16"); + return local_nan; + } + } } // namespace doubleword } // namespace twofloat \ No newline at end of file From 5f7543d017af14f03aa5ecd52faabc7d9369dc41 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 12:23:25 -0500 Subject: [PATCH 014/137] Starting incorporating sin_taylor(). --- .../arithmetics/double-word-arithmetic.hpp | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 81a2e23..5b028fa 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -240,12 +240,14 @@ const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); static const float fp32_nan = std::numeric_limits::quiet_NaN(); +const double fp32_eps = 4.9303806e-32; // 2^-104 // TODO: check it was correctly copied from QD const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); static const double fp64_nan = std::numeric_limits::quiet_NaN(); +const double fp64_eps = 4.93038065763132e-32; // 2^-104 // Reference: QD / inline.h /* Computes the nearest integer to input. */ @@ -267,6 +269,34 @@ inline T nint(T input) { return std::floor(input + pointfive); } +// Reference: QD / dd_inline.h +/* Cast to double. */ +template +inline T to_double(const two &input) { + return input.h; +} + +// Reference: QD / dd_real.cpp +/* Computes sin(a) using Taylor series. + Assumes |a| <= pi/32 */ +template +static two sin_taylor(const two &input) { + + T pointfive, local_eps; + if constexpr (std::is_same_v) { + pointfive = 0.5f; + local_eps = fp32_eps; + } else if constexpr (std::is_same_v) { + pointfive = 0.5; + local_eps = fp64_eps; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + + const T thresh = pointfive * std::abs(to_double(input)) * local_eps; + +} + // Reference: QD / dd_real.cpp template inline two sin(const two &input) { From f3d34f4a188918842ce7ab80a0360252cf84c203 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 12:32:40 -0500 Subject: [PATCH 015/137] Adding prototype for sqr(). --- .../arithmetics/double-word-arithmetic.hpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 5b028fa..08cf20e 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -276,6 +276,12 @@ inline T to_double(const two &input) { return input.h; } +// Reference: QD / dd_inline.h +template +inline two sqr(const two &input) { + +} + // Reference: QD / dd_real.cpp /* Computes sin(a) using Taylor series. Assumes |a| <= pi/32 */ @@ -295,6 +301,16 @@ static two sin_taylor(const two &input) { const T thresh = pointfive * std::abs(to_double(input)) * local_eps; + two r, s, t, x; + + // TODO: make sure this is already supported in twofloat + if (input.is_zero()) { + return 0.0; + } + + int i = 0; + x = -sqr(input); + } // Reference: QD / dd_real.cpp From b347a47a9d7979f4718e5a7f559a45bdad26b904 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 14:27:26 -0500 Subject: [PATCH 016/137] Adding implementation of two_sqr(). Missing implementation of split(). --- .../arithmetics/double-word-arithmetic.hpp | 25 ++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 08cf20e..a50c6d2 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -276,10 +276,33 @@ inline T to_double(const two &input) { return input.h; } +// Reference: QD / inline.h +/* Computes fl(input * input) and err(input * input) */ +template +inline T two_sqr(T input, T &err) { + + T twopointzero; + if constexpr (std::is_same_v) { + twopointzero = 2.0f; + } else if constexpr (std::is_same_v) { + twopointzero = 2.0; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + + T hi, lo; + T q = input * input; + split(input, hi, lo); + err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; + return q; +} + // Reference: QD / dd_inline.h template inline two sqr(const two &input) { - + T p1, p2; + T s1, s2; + p1 = two_sqr(input.h, p2); } // Reference: QD / dd_real.cpp From 959e99a1bc7f99ab208411f45854feedad2d9dc9 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 14:31:03 -0500 Subject: [PATCH 017/137] Adding split() function signature. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index a50c6d2..8d87aee 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -276,6 +276,13 @@ inline T to_double(const two &input) { return input.h; } +// Reference: QD / inline.h +/* Computes high word and low word of input */ +template +inline void split(T input, T &hi, T &lo) { + +} + // Reference: QD / inline.h /* Computes fl(input * input) and err(input * input) */ template From 3c8c66dc413cf6c6a3dc81bb3b0fee7a23c9ab1f Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 15:06:07 -0500 Subject: [PATCH 018/137] Adding split() implementation. --- .../arithmetics/double-word-arithmetic.hpp | 33 ++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 8d87aee..c7c94e9 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -236,18 +236,22 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp / inline.h // TODO: transform this to FP32 (this is FP64) +// TODO: check ranges (i.e., exponents) const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); static const float fp32_nan = std::numeric_limits::quiet_NaN(); -const double fp32_eps = 4.9303806e-32; // 2^-104 +const float fp32_eps = 4.9303806e-32; // 2^-104 +const double fp32_qd_split_thresh = 6.6969287e+299 // = 2^996 // TODO: check it was correctly copied from QD +// TODO: check ranges (i.e., exponents) const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); static const double fp64_nan = std::numeric_limits::quiet_NaN(); const double fp64_eps = 4.93038065763132e-32; // 2^-104 +const double fp64_qd_split_thresh = 6.69692879491417e+299 // = 2^996 // Reference: QD / inline.h /* Computes the nearest integer to input. */ @@ -280,7 +284,34 @@ inline T to_double(const two &input) { /* Computes high word and low word of input */ template inline void split(T input, T &hi, T &lo) { + T temp; + T local_qd_split_thresh; + T local_split_factor; + if constexpr (std::is_same_v) { + local_qd_split_thresh = fp32_qd_split_thresh; + local_split_factor = 3.7252902e-09; // 2^-28 + } else if constexpr (std::is_same_v) { + local_qd_split_thresh = fp64_qd_split_thresh; + local_split_factor = 3.7252902984619140625e-09; // 2^-28 + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + + const T qd_splitter = 134217729.0; // = 2^27 + 1 + + if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { + input *= local_split_factor; // 2^-28 + temp = qd_splitter * input; + hi = temp - (temp - input); + lo = input - hi; + hi *= 268435456.0; // 2^28 + lo *= 268435456.0; // 2^28 + } else { + temp = qd_splitter * input; + hi = temp - (temp - input); + lo = input - hi; + } } // Reference: QD / inline.h From bf736001d43386ef3b00b433d855fdb4811cf8b3 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 15:10:14 -0500 Subject: [PATCH 019/137] More implementation of sqr() untill call to quick_two_sum(). --- .../arithmetics/double-word-arithmetic.hpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index c7c94e9..a9ef5f3 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -338,9 +338,22 @@ inline T two_sqr(T input, T &err) { // Reference: QD / dd_inline.h template inline two sqr(const two &input) { + + T twopointzero; + if constexpr (std::is_same_v) { + twopointzero = 2.0f; + } else if constexpr (std::is_same_v) { + twopointzero = 2.0; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + T p1, p2; T s1, s2; p1 = two_sqr(input.h, p2); + p2 += twopointzero * input.h * input.l; + p2 += input.l * input.l; + s1 = quick_two_sum(p1, p2, s2); } // Reference: QD / dd_real.cpp From 91e42bb5f69eb894d9ff877b8ed8b96a68b92487 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 15:15:21 -0500 Subject: [PATCH 020/137] Implementation of quick_two_sum(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index a9ef5f3..20ed87b 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -335,6 +335,16 @@ inline T two_sqr(T input, T &err) { return q; } +// Reference: QD / inline.h +/* Computes fl(a+b) and err(a+b). + Assumes |a| >= |b|. */ +template +inline T quick_two_sum(T a, T b, T &err) { + T s = a + b; + err = b - (s - a); + return s; +} + // Reference: QD / dd_inline.h template inline two sqr(const two &input) { From a868ea1974c1e7702aafa6a623e1ba94f81a45c3 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 15:17:21 -0500 Subject: [PATCH 021/137] Completing sqr() implementation. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 20ed87b..038909c 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -364,6 +364,7 @@ inline two sqr(const two &input) { p2 += twopointzero * input.h * input.l; p2 += input.l * input.l; s1 = quick_two_sum(p1, p2, s2); + return two(s1, s2); } // Reference: QD / dd_real.cpp From e63594a35c600b3841cf592428f0379d95dde200 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 15:21:08 -0500 Subject: [PATCH 022/137] Extra statements. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 038909c..70ed1f8 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -395,6 +395,8 @@ static two sin_taylor(const two &input) { int i = 0; x = -sqr(input); + s = input; + r = input; } From 7e240cd4ad2e930945833008d34c33df7b327029 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 15:52:54 -0500 Subject: [PATCH 023/137] Completing preliminary implementation of sin_taylor(). --- .../arithmetics/double-word-arithmetic.hpp | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 70ed1f8..d79002d 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -253,6 +253,45 @@ static const double fp64_nan = std::numeric_limits::quiet_NaN(); const double fp64_eps = 4.93038065763132e-32; // 2^-104 const double fp64_qd_split_thresh = 6.69692879491417e+299 // = 2^996 +// Reference: QD / dd_real.cpp +static const int n_inv_fact = 15; + +static const float fp32_inv_fact[n_inv_fact][2] = { + { 1.6666666e-01, 9.2518585e-18}, + { 4.1666666e-02, 2.3129646e-18}, + { 8.3333333e-03, 1.1564823e-19}, + { 1.3888888e-03, -5.3005439e-20}, + { 1.9841269e-04, 1.7209558e-22}, + { 2.4801587e-05, 2.1511947e-23}, + { 2.7557319e-06, -1.8583932e-22}, + { 2.7557319e-07, 2.3767714e-23}, + { 2.5052108e-08, -1.4488140e-24}, + { 2.0876756e-09, -1.2073450e-25}, + { 1.6059043e-10, 1.2585294e-26}, + { 1.1470745e-11, 2.0655512e-28}, + { 7.6471637e-13, 7.0387287e-30}, + { 4.7794773e-14, 4.3992054e-31}, + { 2.8114572e-15, 1.6508842e-31} +}; + +static const double fp64_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 / inline.h /* Computes the nearest integer to input. */ template @@ -374,12 +413,15 @@ template static two sin_taylor(const two &input) { T pointfive, local_eps; + T* local_ptr_inv_fact; if constexpr (std::is_same_v) { pointfive = 0.5f; local_eps = fp32_eps; + local_ptr_inv_fact = &fp32_inv_fact[0][0]; } else if constexpr (std::is_same_v) { pointfive = 0.5; local_eps = fp64_eps; + local_ptr_inv_fact = &fp64_inv_fact[0][0]; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } @@ -398,6 +440,12 @@ static two sin_taylor(const two &input) { 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(to_double(t)) > thresh); } // Reference: QD / dd_real.cpp From ce0431936003e5a512b154fc8da8e966cc2b1107 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:02:14 -0500 Subject: [PATCH 024/137] Starting to incorporate cos_taylor() implementation. --- .../arithmetics/double-word-arithmetic.hpp | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index d79002d..4bfe9ab 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -448,6 +448,40 @@ static two sin_taylor(const two &input) { } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); } +template +static two cos_taylor(const two &input) { + + T pointfive, local_eps, onepointzero; + T* local_ptr_inv_fact; + if constexpr (std::is_same_v) { + pointfive = 0.5f; + local_eps = fp32_eps; + onepointzero = 1.0f + local_ptr_inv_fact = &fp32_inv_fact[0][0]; + } else if constexpr (std::is_same_v) { + pointfive = 0.5; + local_eps = fp64_eps; + onepointzero = 1.0; + local_ptr_inv_fact = &fp64_inv_fact[0][0]; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + + const T thresh = pointfive * local_eps; + + two r, s, t, x; + + // TODO: make sure this is already supported in twofloat + if (input.is_zero()) { + return onepointzero; + } + + x = -sqr(input); + r = x; + s = 1.0 + mul_pwr2(r, pointfive); + int i = 1; +} + // Reference: QD / dd_real.cpp template inline two sin(const two &input) { From 1cb051ca46134f021dd13e6486591e3b9f5e7cdf Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:02:58 -0500 Subject: [PATCH 025/137] Minor fix. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 4bfe9ab..d29964e 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -478,7 +478,7 @@ static two cos_taylor(const two &input) { x = -sqr(input); r = x; - s = 1.0 + mul_pwr2(r, pointfive); + s = add(onepointzero, mul_pwr2(r, pointfive)); int i = 1; } From 69d44d3175ab370587e943bfbca8abac47b5d077 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:07:42 -0500 Subject: [PATCH 026/137] Adding implementation of mul_pwr2(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index d29964e..354a6df 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -448,6 +448,13 @@ static two sin_taylor(const two &input) { } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); } +// Reference: QD / dd_inline.h +/* double-double * double, where double is a power of 2. */ +template +inline two mul_pwr2(const two &input, T b) { + return two(input.h * b, input.l * b); +} + template static two cos_taylor(const two &input) { From f6e9e58ed432c31b870019226a815f12de6dac40 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:10:28 -0500 Subject: [PATCH 027/137] Completing implementation of cos_taylor(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 354a6df..4e8ac20 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -444,7 +444,7 @@ static two sin_taylor(const two &input) { 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; + i += 2; } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); } @@ -487,6 +487,12 @@ static two cos_taylor(const two &input) { 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(to_double(t)) > thresh); } // Reference: QD / dd_real.cpp From 53f11cfdf43d418812ee4319b527a79dd07492ea Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:13:32 -0500 Subject: [PATCH 028/137] More of sin(). --- .../arithmetics/double-word-arithmetic.hpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 4e8ac20..60e4c2e 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -546,6 +546,18 @@ inline two sin(const two &input) { return local_nan; } + if (k == 0) { + switch(j) { + case 0: + return sin_taylor(t); + case 1: + return cos_taylor(t); + case -1: + return -cos_taylor(t); + default: + return -sin_taylor(t); + } + } } } // namespace doubleword } // namespace twofloat \ No newline at end of file From 8140ec2efe43cfdced8f542acc53d8ce70efadef Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:28:44 -0500 Subject: [PATCH 029/137] Adding sin and cos tables for sin() implementation. --- .../arithmetics/double-word-arithmetic.hpp | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 60e4c2e..28fa325 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -292,6 +292,38 @@ static const double fp64_inv_fact[n_inv_fact][2] = { { 2.81145725434552060e-15, 1.65088427308614326e-31} }; +/* Table of sin(k * pi/16) and cos(k * pi/16). */ + +// Reference: QD / dd_real.cpp +static const float fp32_cos_table [4][2] = { + {9.8078528e-01, 1.8546939e-17}, + {9.2387953e-01, 1.7645047e-17}, + {8.3146961e-01, 1.4073856e-18}, + {7.0710678e-01, -4.8336466e-17} +}; + +static const double fp64_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} +}; + +// Reference: QD / dd_real.cpp +static const single fp32_sin_table [4][2] = { + {1.9509032e-01, -7.9910790e-18}, + {3.8268343e-01, -1.0050772e-17}, + {5.5557023e-01, 4.7094109e-17}, + {7.0710678e-01, -4.8336466e-17} +}; + +static const double fp64_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} +}; + // Reference: QD / inline.h /* Computes the nearest integer to input. */ template @@ -505,18 +537,23 @@ inline two sin(const two &input) { } T local_2pi, local_pi2, pointfive, local_pi16, local_nan; + T* local_ptr_cos_table, local_ptr_sin_table; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; local_pi2 = fp32_pi2; pointfive = 0.5f; local_pi16 = fp32_pi16: local_nan = fp32_nan; + local_ptr_cos_table = &fp32_cos_table[0][0]; + local_ptr_sin_table = &fp32_sin_table[0][0]; } else if constexpr (std::is_same_v) { local_2pi = fp64_2pi; local_pi2 = fp64_pi2; pointfive = 0.5; local_pi16 = fp64_pi16; local_nan = fp64_nan; + local_ptr_cos_table = &fp64_cos_table[0][0]; + local_ptr_sin_table = &fp64_sin_table[0][0]; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } @@ -558,6 +595,10 @@ inline two sin(const two &input) { return -sin_taylor(t); } } + + two u(local_ptr_cos_table[abs_k - 1][0], local_ptr_cos_table[abs_k - 1][1]); + two u(local_ptr_sin_table[abs_k - 1][0], local_ptr_sin_table[abs_k - 1][1]); + two sin_t, cos_t; } } // namespace doubleword } // namespace twofloat \ No newline at end of file From 0a3bd0d8ff76b1f61eec1273e95fdd299f9b9d0a Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:36:21 -0500 Subject: [PATCH 030/137] Adding implementation of sincos_taylor(). --- .../arithmetics/double-word-arithmetic.hpp | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 28fa325..59744b6 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -527,6 +527,30 @@ static two cos_taylor(const two &input) { } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); } +template +static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { + + T zeropointzero, onepointzero; + if constexpr (std::is_same_v) { + zeropointzero = 0.0f; + onepointzero = 1.0f + } else if constexpr (std::is_same_v) { + zeropointzero = 0.0; + onepointzero = 1.0; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + + if (input.is_zero()) { + sin_a = zeropointzero; + cos_a = onepointzero; + return; + } + + sin_a = sin_taylor(input); + cos_a = sqrt(onepointzero - sqr(sin_a)); +} + // Reference: QD / dd_real.cpp template inline two sin(const two &input) { From 75d1293024e533eae06dddfc87b66e1eac7589f2 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 23 Jan 2024 16:44:06 -0500 Subject: [PATCH 031/137] Completing sin() implementation. It was not compiled yet. --- .../arithmetics/double-word-arithmetic.hpp | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 59744b6..c7c6aea 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -621,8 +621,36 @@ inline two sin(const two &input) { } two u(local_ptr_cos_table[abs_k - 1][0], local_ptr_cos_table[abs_k - 1][1]); - two u(local_ptr_sin_table[abs_k - 1][0], local_ptr_sin_table[abs_k - 1][1]); + two v(local_ptr_sin_table[abs_k - 1][0], local_ptr_sin_table[abs_k - 1][1]); two sin_t, cos_t; + sincos_taylor(t, sin_t, cos_t); + 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)); + } + } else if (j == -1) { + if (k > 0) { + r = sub(mul(v, sin_t), mul(u, cos_t)); + } else if (k < 0) { + r = sub(-mul(u, cos_t), mul(v, sin_t)); + } + } else { + if (k > 0) { + r = sub(-mul(u, sin_t), mul(v, cos_t)); + } else { + r = sub(mul(v, cos_t), mul(u, sin_t)); + } + } + + return r; } } // namespace doubleword } // namespace twofloat \ No newline at end of file From 01817faaa94d7f861729a54ecf642437d4839861 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 08:02:18 -0500 Subject: [PATCH 032/137] Improve/remove comments. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index c7c6aea..b72576b 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -234,7 +234,7 @@ inline two div(const two &x, const two &y) { static_assert(sizeof(T) == 0, "Unsupported mode"); } -// Reference: QD / dd_const.cpp / inline.h +// Reference: QD / dd_const.cpp / inline.h / dd_real.cpp // TODO: transform this to FP32 (this is FP64) // TODO: check ranges (i.e., exponents) const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); @@ -244,14 +244,12 @@ static const float fp32_nan = std::numeric_limits::quiet_NaN(); const float fp32_eps = 4.9303806e-32; // 2^-104 const double fp32_qd_split_thresh = 6.6969287e+299 // = 2^996 -// TODO: check it was correctly copied from QD -// TODO: check ranges (i.e., exponents) const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); static const double fp64_nan = std::numeric_limits::quiet_NaN(); const double fp64_eps = 4.93038065763132e-32; // 2^-104 -const double fp64_qd_split_thresh = 6.69692879491417e+299 // = 2^996 +const double fp64_qd_split_thresh = 6.69692879491417e+299 // = 2^996 // TODO: check! range is beyond double // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; From 5ca0c179d8c57175076f68b8fc0660281d47951e Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 08:05:18 -0500 Subject: [PATCH 033/137] Improving and correcting definitions of cos and sin tables. --- .../arithmetics/double-word-arithmetic.hpp | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index b72576b..cd5aedc 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -290,15 +290,20 @@ static const double fp64_inv_fact[n_inv_fact][2] = { { 2.81145725434552060e-15, 1.65088427308614326e-31} }; -/* Table of sin(k * pi/16) and cos(k * pi/16). */ - // Reference: QD / dd_real.cpp +/* Table of sin(k * pi/16) and cos(k * pi/16). */ static const float fp32_cos_table [4][2] = { {9.8078528e-01, 1.8546939e-17}, {9.2387953e-01, 1.7645047e-17}, {8.3146961e-01, 1.4073856e-18}, {7.0710678e-01, -4.8336466e-17} }; +static const float fp32_sin_table [4][2] = { + {1.9509032e-01, -7.9910790e-18}, + {3.8268343e-01, -1.0050772e-17}, + {5.5557023e-01, 4.7094109e-17}, + {7.0710678e-01, -4.8336466e-17} +}; static const double fp64_cos_table [4][2] = { {9.807852804032304306e-01, 1.854693999782500573e-17}, @@ -306,15 +311,6 @@ static const double fp64_cos_table [4][2] = { {8.314696123025452357e-01, 1.407385698472802389e-18}, {7.071067811865475727e-01, -4.833646656726456726e-17} }; - -// Reference: QD / dd_real.cpp -static const single fp32_sin_table [4][2] = { - {1.9509032e-01, -7.9910790e-18}, - {3.8268343e-01, -1.0050772e-17}, - {5.5557023e-01, 4.7094109e-17}, - {7.0710678e-01, -4.8336466e-17} -}; - static const double fp64_sin_table [4][2] = { {1.950903220161282758e-01, -7.991079068461731263e-18}, {3.826834323650897818e-01, -1.005077269646158761e-17}, From e7ba4f4b70c34872ac5a7897e482c5a8a7c7f994 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 08:34:46 -0500 Subject: [PATCH 034/137] Fixing till split(). --- .../arithmetics/double-word-arithmetic.hpp | 34 ++++++++++++------- 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index cd5aedc..4800d03 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -235,21 +235,27 @@ inline two div(const two &x, const two &y) { } // Reference: QD / dd_const.cpp / inline.h / dd_real.cpp -// TODO: transform this to FP32 (this is FP64) +// TODO: transform this to FP32 (the original code is FP64) // TODO: check ranges (i.e., exponents) const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); static const float fp32_nan = std::numeric_limits::quiet_NaN(); const float fp32_eps = 4.9303806e-32; // 2^-104 -const double fp32_qd_split_thresh = 6.6969287e+299 // = 2^996 +const float fp32_qd_splitter = 134217729.0f; // = 2^27 + 1 +const float fp32_qd_split_thresh = 6.6969287e+299; // = 2^996 // TODO: check! range is beyond float +const float fp32_split_factor = 3.7252902e-09; // 2^-28 +const float fp32_split_factor_2 = 268435456.0f; // 2^28 const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); static const double fp64_nan = std::numeric_limits::quiet_NaN(); const double fp64_eps = 4.93038065763132e-32; // 2^-104 -const double fp64_qd_split_thresh = 6.69692879491417e+299 // = 2^996 // TODO: check! range is beyond double +const double fp64_qd_splitter = 134217729.0; // = 2^27 + 1 +const double fp64_qd_split_thresh = 6.69692879491417e+299; // = 2^996 // TODO: check! range is beyond double +const double fp64_split_factor = 3.7252902984619140625e-09; // 2^-28 +const double fp64_split_factor_2 = 268435456.0; // 2^28 // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -351,29 +357,33 @@ template inline void split(T input, T &hi, T &lo) { T temp; + T local_qd_splitter; T local_qd_split_thresh; T local_split_factor; + T local_split_factor_2; if constexpr (std::is_same_v) { + local_qd_splitter = fp32_qd_splitter; local_qd_split_thresh = fp32_qd_split_thresh; - local_split_factor = 3.7252902e-09; // 2^-28 + local_split_factor = fp32_split_factor; + local_split_factor_2 = fp32_split_factor_2; } else if constexpr (std::is_same_v) { + local_qd_splitter = fp64_qd_splitter; local_qd_split_thresh = fp64_qd_split_thresh; - local_split_factor = 3.7252902984619140625e-09; // 2^-28 + local_split_factor = fp64_split_factor; + local_split_factor_2 = fp64_split_factor_2; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } - const T qd_splitter = 134217729.0; // = 2^27 + 1 - if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { - input *= local_split_factor; // 2^-28 - temp = qd_splitter * input; + input *= local_split_factor; + temp = local_qd_splitter * input; hi = temp - (temp - input); lo = input - hi; - hi *= 268435456.0; // 2^28 - lo *= 268435456.0; // 2^28 + hi *= local_split_factor_2; + lo *= local_split_factor_2; } else { - temp = qd_splitter * input; + temp = local_qd_splitter * input; hi = temp - (temp - input); lo = input - hi; } From 78e1f108d65336b8ab3af73282f4700bc5857160 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 08:47:51 -0500 Subject: [PATCH 035/137] Minor improvements. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 4800d03..9f6da74 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -410,9 +410,9 @@ inline T two_sqr(T input, T &err) { return q; } +// TODO: already implemented in twofloat? // Reference: QD / inline.h -/* Computes fl(a+b) and err(a+b). - Assumes |a| >= |b|. */ +/* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */ template inline T quick_two_sum(T a, T b, T &err) { T s = a + b; From ce47dfed7629bc4a3c785d435f3bbe33fa7a957c Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 09:09:39 -0500 Subject: [PATCH 036/137] Minor addition. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 9f6da74..e9c29e9 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -491,6 +491,7 @@ 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) { From c9dd711457fc04a6bacab787377c7c2cdc3e1489 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 09:56:33 -0500 Subject: [PATCH 037/137] Adding zeropointzero. --- .../arithmetics/double-word-arithmetic.hpp | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index e9c29e9..9b670b3 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -450,14 +450,17 @@ static two sin_taylor(const two &input) { T pointfive, local_eps; T* local_ptr_inv_fact; + T zeropointzero; if constexpr (std::is_same_v) { pointfive = 0.5f; local_eps = fp32_eps; local_ptr_inv_fact = &fp32_inv_fact[0][0]; + zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { pointfive = 0.5; local_eps = fp64_eps; local_ptr_inv_fact = &fp64_inv_fact[0][0]; + zeropointzero = 0.0; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } @@ -468,7 +471,7 @@ static two sin_taylor(const two &input) { // TODO: make sure this is already supported in twofloat if (input.is_zero()) { - return 0.0; + return zeropointzero; } int i = 0; @@ -560,13 +563,9 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { template inline two sin(const two &input) { - // TODO: make sure this is already supported in twofloat - if (input.is_zero()) { - return 0.0; - } - T local_2pi, local_pi2, pointfive, local_pi16, local_nan; T* local_ptr_cos_table, local_ptr_sin_table; + T zeropointzero; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; local_pi2 = fp32_pi2; @@ -575,6 +574,7 @@ inline two sin(const two &input) { local_nan = fp32_nan; local_ptr_cos_table = &fp32_cos_table[0][0]; local_ptr_sin_table = &fp32_sin_table[0][0]; + zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { local_2pi = fp64_2pi; local_pi2 = fp64_pi2; @@ -583,10 +583,16 @@ inline two sin(const two &input) { local_nan = fp64_nan; local_ptr_cos_table = &fp64_cos_table[0][0]; local_ptr_sin_table = &fp64_sin_table[0][0]; + zeropointzero = 0.0; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } + // TODO: make sure this is already supported in twofloat + if (input.is_zero()) { + return zeropointzero; + } + // Approximately reducing modulo 2*pi two z = nint(input / local_2pi); two r = input - mul(local_2pi, z); From fce6fd3b5b377224f36d6e3b0bd5bc1d892edfa5 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 09:59:32 -0500 Subject: [PATCH 038/137] Adding missing return to sin_taylor() and cos_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 9b670b3..fe5af48 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -485,6 +485,8 @@ static two sin_taylor(const two &input) { s = add(s, t); i += 2; } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); + + return s; } // Reference: QD / dd_inline.h @@ -533,6 +535,8 @@ static two cos_taylor(const two &input) { s = add(s, t); i += 2; } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); + + return s; } template From 3c19973c0e4da9c236a8410c06f3ba02265ce226 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 10:18:51 -0500 Subject: [PATCH 039/137] Starting to add implementation of sqrt(). --- .../arithmetics/double-word-arithmetic.hpp | 35 +++++++++++++++++-- 1 file changed, 33 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index fe5af48..40abd22 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -539,13 +539,44 @@ static two cos_taylor(const two &input) { return s; } +// Reference: QD / dd_real.cpp +/* Computes the square root of the double-double number dd. + NOTE: dd must be a non-negative number. */ +template +two sqrt(const two &input) { + /* 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. + */ + + T zeropointzero; + if constexpr (std::is_same_v) { + zeropointzero = 0.0f; + } else if constexpr (std::is_same_v) { + zeropointzero = 0.0; + } else { + std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + } + + if (input.is_zero()) { + return zeropointzero; + } + +} + +// Reference: QD / dd_real.cpp template static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { T zeropointzero, onepointzero; if constexpr (std::is_same_v) { zeropointzero = 0.0f; - onepointzero = 1.0f + onepointzero = 1.0f; } else if constexpr (std::is_same_v) { zeropointzero = 0.0; onepointzero = 1.0; @@ -560,7 +591,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { } sin_a = sin_taylor(input); - cos_a = sqrt(onepointzero - sqr(sin_a)); + cos_a = sqrt(sub(onepointzero, sqr(sin_a))); } // Reference: QD / dd_real.cpp From 1777f07477f58d21a7d3c057e5790b241572d824 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 10:42:41 -0500 Subject: [PATCH 040/137] Addind dd_add(). --- .../arithmetics/double-word-arithmetic.hpp | 31 ++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 40abd22..ec1ede2 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -539,6 +539,20 @@ static two cos_taylor(const two &input) { return s; } +// Reference QD / inline.h +template +inline T two_sum(T a, T b, T &err) { + +} + +// Reference: QD / dd_inline.h +template +inline two dd_add(T a, T b) { + T s, e; + s = two_sum(a, b, e); // TODO: make sure this has not been already implemented + return two (s, e); +} + // Reference: QD / dd_real.cpp /* Computes the square root of the double-double number dd. NOTE: dd must be a non-negative number. */ @@ -554,11 +568,17 @@ two sqrt(const two &input) { only half the precision. */ - T zeropointzero; + T zeropointzero, onepointzero, pointfive, local_nan; if constexpr (std::is_same_v) { + onepointzero = 1.0f; zeropointzero = 0.0f; + pointfive = 0.5f; + local_nan = fp32_nan; } else if constexpr (std::is_same_v) { + onepointzero = 1.0; zeropointzero = 0.0; + pointfive = 0.5; + local_nan = fp64_nan; } else { std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } @@ -567,6 +587,15 @@ two sqrt(const two &input) { return zeropointzero; } + if (input.is_negative()) { + std::error("LSV: sqrt(): negative argument"); // TODO: make sure std::error is best way to proceed here + return local_nan; + } + + T x = onepointzero / std::sqrt(input.h); + T ax = mul(input.h, x); + two temp = sub(input, sqr(ax)); + return dd_add(ax, temp.h * x * pointfive); } // Reference: QD / dd_real.cpp From d98da8025e543fc8056ea4fcc15afe8a04d9e423 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 10:46:42 -0500 Subject: [PATCH 041/137] Completing two_sum() implementation. --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index ec1ede2..fdb463c 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -540,16 +540,20 @@ static two cos_taylor(const two &input) { } // Reference QD / inline.h +// TODO: make sure this has not been already implemented 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 template inline two dd_add(T a, T b) { T s, e; - s = two_sum(a, b, e); // TODO: make sure this has not been already implemented + s = two_sum(a, b, e); return two (s, e); } From e9939c809d88b39abfb494c9b27926db805dca07 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 10:51:50 -0500 Subject: [PATCH 042/137] Using ops for two<> types. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index fdb463c..4759e81 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -663,7 +663,7 @@ inline two sin(const two &input) { // Approximately reducing modulo 2*pi two z = nint(input / local_2pi); - two r = input - mul(local_2pi, z); + two r = sub(input, mul(local_2pi, z)); // Approximately reducing modulo pi/2 and then modulo pi/16 From 057c25c4b9f5907c96670ce37929c6c2b3abdec3 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 10:57:30 -0500 Subject: [PATCH 043/137] Using ops for two<> types. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 4759e81..cab1260 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -672,7 +672,7 @@ inline two sin(const two &input) { two t = r - mul(local_pi2, q); int j = static_cast(q); q = std::floor(t.h / local_pi16.h + pointfive); - t -= mul(local_pi16, q); + t = sub(t, mul(local_pi16, q)); int k = static_cast(q); int abs_k = std:abs(k); From 0e3e8cc95e6c295f7864fb7f7325f8dead8685cb Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 15:01:45 -0500 Subject: [PATCH 044/137] Commenting out std::error + fixing syntax. --- .../arithmetics/double-word-arithmetic.hpp | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index cab1260..e929a69 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -335,7 +335,7 @@ inline T nint(T input) { } else if constexpr (std::is_same_v) { pointfive = 0.5; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } if(input = std::floor(input)) { @@ -372,7 +372,7 @@ inline void split(T input, T &hi, T &lo) { local_split_factor = fp64_split_factor; local_split_factor_2 = fp64_split_factor_2; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { @@ -400,7 +400,7 @@ inline T two_sqr(T input, T &err) { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } T hi, lo; @@ -430,7 +430,7 @@ inline two sqr(const two &input) { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } T p1, p2; @@ -462,7 +462,7 @@ static two sin_taylor(const two &input) { local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } const T thresh = pointfive * std::abs(to_double(input)) * local_eps; @@ -505,7 +505,7 @@ static two cos_taylor(const two &input) { if constexpr (std::is_same_v) { pointfive = 0.5f; local_eps = fp32_eps; - onepointzero = 1.0f + onepointzero = 1.0f; local_ptr_inv_fact = &fp32_inv_fact[0][0]; } else if constexpr (std::is_same_v) { pointfive = 0.5; @@ -513,7 +513,7 @@ static two cos_taylor(const two &input) { onepointzero = 1.0; local_ptr_inv_fact = &fp64_inv_fact[0][0]; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } const T thresh = pointfive * local_eps; @@ -584,7 +584,7 @@ two sqrt(const two &input) { pointfive = 0.5; local_nan = fp64_nan; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } if (input.is_zero()) { @@ -592,7 +592,7 @@ two sqrt(const two &input) { } if (input.is_negative()) { - std::error("LSV: sqrt(): negative argument"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: sqrt(): negative argument"); // TODO: make sure std::error is best way to proceed here return local_nan; } @@ -614,7 +614,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { zeropointzero = 0.0; onepointzero = 1.0; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } if (input.is_zero()) { @@ -638,7 +638,7 @@ inline two sin(const two &input) { local_2pi = fp32_2pi; local_pi2 = fp32_pi2; pointfive = 0.5f; - local_pi16 = fp32_pi16: + local_pi16 = fp32_pi16; local_nan = fp32_nan; local_ptr_cos_table = &fp32_cos_table[0][0]; local_ptr_sin_table = &fp32_sin_table[0][0]; @@ -653,7 +653,7 @@ inline two sin(const two &input) { local_ptr_sin_table = &fp64_sin_table[0][0]; zeropointzero = 0.0; } else { - std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here +// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here } // TODO: make sure this is already supported in twofloat @@ -674,15 +674,15 @@ inline two sin(const two &input) { 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); + int abs_k = std::abs(k); if (j < -2 || j > 2) { - std::error("LSV: cannot reduce modulo pi/2"); +// std::error("LSV: cannot reduce modulo pi/2"); return local_nan; } if (abs_k > 4) { - std::error("LSV: cannot reduce modulo pi/16"); +// std::error("LSV: cannot reduce modulo pi/16"); return local_nan; } From 10b124f249a749751fce441fd202d96f878b6a5c Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 16:03:58 -0500 Subject: [PATCH 045/137] Replacing most commented-out std::error with static_assert(). --- .../arithmetics/double-word-arithmetic.hpp | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index e929a69..41583a1 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -335,7 +335,7 @@ inline T nint(T input) { } else if constexpr (std::is_same_v) { pointfive = 0.5; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } if(input = std::floor(input)) { @@ -372,7 +372,7 @@ inline void split(T input, T &hi, T &lo) { local_split_factor = fp64_split_factor; local_split_factor_2 = fp64_split_factor_2; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { @@ -400,7 +400,7 @@ inline T two_sqr(T input, T &err) { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } T hi, lo; @@ -430,7 +430,7 @@ inline two sqr(const two &input) { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } T p1, p2; @@ -462,7 +462,7 @@ static two sin_taylor(const two &input) { local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } const T thresh = pointfive * std::abs(to_double(input)) * local_eps; @@ -513,7 +513,7 @@ static two cos_taylor(const two &input) { onepointzero = 1.0; local_ptr_inv_fact = &fp64_inv_fact[0][0]; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } const T thresh = pointfive * local_eps; @@ -584,7 +584,7 @@ two sqrt(const two &input) { pointfive = 0.5; local_nan = fp64_nan; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } if (input.is_zero()) { @@ -592,7 +592,7 @@ two sqrt(const two &input) { } if (input.is_negative()) { -// std::error("LSV: sqrt(): negative argument"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); return local_nan; } @@ -614,7 +614,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { zeropointzero = 0.0; onepointzero = 1.0; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } if (input.is_zero()) { @@ -653,7 +653,7 @@ inline two sin(const two &input) { local_ptr_sin_table = &fp64_sin_table[0][0]; zeropointzero = 0.0; } else { -// std::error("LSV: other types are unsupported"); // TODO: make sure std::error is best way to proceed here + static_assert(sizeof(T) == 0, "LSV: other types not supported"); } // TODO: make sure this is already supported in twofloat From 060d44c78a4f373ddb84f2ebb2f8832b7b76cebb Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 16:09:39 -0500 Subject: [PATCH 046/137] Remove --- test/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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) From d38ee32677bb380696e29a329906162346ff58a6 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 16:51:02 -0500 Subject: [PATCH 047/137] Starting to add a very simple test function for sin(). --- test/double-word-arithmetic.test.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index d85614b..fcc994a 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -255,6 +255,9 @@ TEST(DoubleWordArithmetic, DivFPTest) { divFPTest(); } TEST(DoubleWordArithmetic, DivFPFMATest) { divFPTest(); } +TEST(DoubleWordArithmetic, SinTest) {} + float x = 10; + float sinx = std::sin(x); } // namespace test } // namespace doubleword } // namespace twofloat \ No newline at end of file From bb18240ccde98df23872217e8603318c5cc8eb42 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 16:55:02 -0500 Subject: [PATCH 048/137] Starting to add very simple test to sin() --- test/double-word-arithmetic.test.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index fcc994a..5987f9b 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -258,6 +258,9 @@ TEST(DoubleWordArithmetic, DivFPFMATest) { divFPTest(); } TEST(DoubleWordArithmetic, SinTest) {} float x = 10; float sinx = std::sin(x); + + two two_x = {x, 0.0f}; + //float two_sinx = doubleword::sin(two_x); } // namespace test } // namespace doubleword } // namespace twofloat \ No newline at end of file From 4d86992858d7b54dac5549dfab1b2e665a3fb5ff Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 24 Jan 2024 16:57:09 -0500 Subject: [PATCH 049/137] Starting to call two sin(). It fails compilation. --- test/double-word-arithmetic.test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index 5987f9b..31576c1 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -260,7 +260,7 @@ TEST(DoubleWordArithmetic, SinTest) {} float sinx = std::sin(x); two two_x = {x, 0.0f}; - //float two_sinx = doubleword::sin(two_x); + two sin_two_x = doubleword::sin(two_x); } // namespace test } // namespace doubleword } // namespace twofloat \ No newline at end of file From 1daf666f3b97281f653dc91f67fda71154c2512f Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 09:01:29 -0500 Subject: [PATCH 050/137] Starting to fix implementation using compiler info. --- .../arithmetics/double-word-arithmetic.hpp | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 41583a1..e6a17fb 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -351,6 +351,7 @@ inline T to_double(const two &input) { return input.h; } +// TODO: make sure it has not been already implemented in twofloat // Reference: QD / inline.h /* Computes high word and low word of input */ template @@ -631,34 +632,43 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { template inline two sin(const two &input) { - T local_2pi, local_pi2, pointfive, local_pi16, local_nan; + two local_2pi; + T local_pi2, pointfive, local_pi16, local_nan; T* local_ptr_cos_table, local_ptr_sin_table; T zeropointzero; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; + /* local_pi2 = fp32_pi2; pointfive = 0.5f; local_pi16 = fp32_pi16; local_nan = fp32_nan; local_ptr_cos_table = &fp32_cos_table[0][0]; local_ptr_sin_table = &fp32_sin_table[0][0]; + */ zeropointzero = 0.0f; + } else if constexpr (std::is_same_v) { local_2pi = fp64_2pi; + /* local_pi2 = fp64_pi2; pointfive = 0.5; local_pi16 = fp64_pi16; local_nan = fp64_nan; local_ptr_cos_table = &fp64_cos_table[0][0]; local_ptr_sin_table = &fp64_sin_table[0][0]; + */ zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "LSV: other types not supported"); } // TODO: make sure this is already supported in twofloat - if (input.is_zero()) { - return zeropointzero; + //if (input.is_zero()) { + if (input.eval() == zeropointzero) { + T hi, lo; + split(zeropointzero, hi, lo); + return two{hi, lo}; } // Approximately reducing modulo 2*pi @@ -666,7 +676,7 @@ inline two sin(const two &input) { two r = sub(input, mul(local_2pi, z)); // Approximately reducing modulo pi/2 and then modulo pi/16 - +/* //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); two t = r - mul(local_pi2, q); @@ -730,6 +740,7 @@ inline two sin(const two &input) { } return r; +*/ } } // namespace doubleword } // namespace twofloat \ No newline at end of file From 33f420ad71d49ba688aa09ebafb8fcd070d1e029 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 09:50:25 -0500 Subject: [PATCH 051/137] Creating the namespace qd and moving some functions into it. --- .../arithmetics/double-word-arithmetic.hpp | 78 +++++++++++++++---- 1 file changed, 61 insertions(+), 17 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index e6a17fb..fc78079 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -324,24 +324,78 @@ static const double fp64_sin_table [4][2] = { {7.071067811865475727e-01, -4.833646656726456726e-17} }; -// Reference: QD / inline.h -/* Computes the nearest integer to input. */ +// Notice the definition of the qd namespace +namespace qd { + + // Reference: QD / inline.h + /* Computes the nearest integer to input. */ + template + inline T nint(T input) { + T pointfive; + + if constexpr (std::is_same_v) { + pointfive = 0.5f; + } else if constexpr (std::is_same_v) { + pointfive = 0.5; + } else { + static_assert(sizeof(T) == 0, "LSV: other types not supported"); + } + + if(input = std::floor(input)) { + return input; + } + return std::floor(input + pointfive); + } + + // TODO: already implemented in twofloat? + // Reference: QD / inline.h + /* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */ + template + inline T quick_two_sum(T a, T b, T &err) { + T s = a + b; + err = b - (s - a); + return s; + } +} + +// Reference: QD / dd_inline.h +/* Round to nearest integer */ template -inline T nint(T input) { - T pointfive; +inline two nint(const two &input) { + T hi = qd::nint(input.h); + T lo; + T zeropointzero, pointfive, onepointzero; if constexpr (std::is_same_v) { + zeropointzero = 0.0f; pointfive = 0.5f; + onepointzero = 1.0f; } else if constexpr (std::is_same_v) { + zeropointzero = 0.0; pointfive = 0.5; + onepointzero = 1.0; } else { static_assert(sizeof(T) == 0, "LSV: other types not supported"); } - if(input = std::floor(input)) { - return input; + if(hi == input.h) { + /* High word is an integer already. Round the low word. */ + lo = dq::nint(input.l); + + /* Renormalize. This is needed if h = some integer, l = 1/2. */ + hi = qd::quick_two_sum(hi, lo, lo); } - return std::floor(input + pointfive); + 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 @@ -411,16 +465,6 @@ inline T two_sqr(T input, T &err) { return q; } -// TODO: already implemented in twofloat? -// Reference: QD / inline.h -/* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */ -template -inline T quick_two_sum(T a, T b, T &err) { - T s = a + b; - err = b - (s - a); - return s; -} - // Reference: QD / dd_inline.h template inline two sqr(const two &input) { From cb43c9fb3e2a937666f7ac2824049237c4f5b786 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 09:51:38 -0500 Subject: [PATCH 052/137] Fixing namespace. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index fc78079..26d64c3 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -380,7 +380,7 @@ inline two nint(const two &input) { if(hi == input.h) { /* High word is an integer already. Round the low word. */ - lo = dq::nint(input.l); + lo = qd::nint(input.l); /* Renormalize. This is needed if h = some integer, l = 1/2. */ hi = qd::quick_two_sum(hi, lo, lo); From d149699a544da70528a4ed1518dddff0aaa6e3d0 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 09:56:55 -0500 Subject: [PATCH 053/137] Calling div<> with Accurate and FMA options enabled. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 26d64c3..97a41d1 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -677,8 +677,10 @@ template inline two sin(const two &input) { two local_2pi; +/* T local_pi2, pointfive, local_pi16, local_nan; T* local_ptr_cos_table, local_ptr_sin_table; +*/ T zeropointzero; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; @@ -716,7 +718,7 @@ inline two sin(const two &input) { } // Approximately reducing modulo 2*pi - two z = nint(input / local_2pi); + two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div two r = sub(input, mul(local_2pi, z)); // Approximately reducing modulo pi/2 and then modulo pi/16 From 96c6f0d9b88f4ac857c7372f9b9aaa6a87baa06b Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 10:37:56 -0500 Subject: [PATCH 054/137] Fixing sub by adding mode (Accurate). Code compiles successfully so far. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 97a41d1..bd5443b 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -718,8 +718,8 @@ inline two sin(const two &input) { } // Approximately reducing modulo 2*pi - two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div - two r = sub(input, mul(local_2pi, z)); + two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div (accurate + fma) + two r = sub(input, mul(local_2pi, z)); // Approximately reducing modulo pi/2 and then modulo pi/16 /* From 22f41ae6fdd06fc6e41b5f3daf51b4f5eaa319f5 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 10:41:59 -0500 Subject: [PATCH 055/137] Enabling more compiling code. --- .../arithmetics/double-word-arithmetic.hpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index bd5443b..fc25d9c 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -677,16 +677,17 @@ template inline two sin(const two &input) { two local_2pi; + two local_pi2; /* - T local_pi2, pointfive, local_pi16, local_nan; + T pointfive, local_pi16, local_nan; T* local_ptr_cos_table, local_ptr_sin_table; */ - T zeropointzero; + T zeropointzero, pointfive; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; - /* local_pi2 = fp32_pi2; pointfive = 0.5f; + /* local_pi16 = fp32_pi16; local_nan = fp32_nan; local_ptr_cos_table = &fp32_cos_table[0][0]; @@ -696,9 +697,9 @@ inline two sin(const two &input) { } else if constexpr (std::is_same_v) { local_2pi = fp64_2pi; - /* local_pi2 = fp64_pi2; pointfive = 0.5; + /* local_pi16 = fp64_pi16; local_nan = fp64_nan; local_ptr_cos_table = &fp64_cos_table[0][0]; @@ -722,9 +723,10 @@ inline two sin(const two &input) { two r = sub(input, mul(local_2pi, z)); // Approximately reducing modulo pi/2 and then modulo pi/16 -/* + //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); +/* two t = r - mul(local_pi2, q); int j = static_cast(q); q = std::floor(t.h / local_pi16.h + pointfive); From 8bbd5fb8de5eeb130cf85787ade6584b25b77c9b Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 10:52:43 -0500 Subject: [PATCH 056/137] Enabling more compiling code. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index fc25d9c..c547b00 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -726,8 +726,8 @@ inline two sin(const two &input) { //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); + two t = sub(r, mul(local_pi2, q)); /* - two t = 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)); From 918d5203a1d291e74d9d099552316ef6c28e5a54 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 10:53:20 -0500 Subject: [PATCH 057/137] Enabling more compiling code. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index c547b00..8daf67d 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -727,8 +727,8 @@ inline two sin(const two &input) { //TODO: original type is double, here it is templated 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); From d43745f43e1f95ff895e2e2ddfa6fd2158102a2a Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 10:55:42 -0500 Subject: [PATCH 058/137] Enabling more compiling code. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 8daf67d..f66a384 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -678,6 +678,7 @@ inline two sin(const two &input) { two local_2pi; two local_pi2; + two local_pi16; /* T pointfive, local_pi16, local_nan; T* local_ptr_cos_table, local_ptr_sin_table; @@ -687,8 +688,8 @@ inline two sin(const two &input) { local_2pi = fp32_2pi; local_pi2 = fp32_pi2; pointfive = 0.5f; - /* local_pi16 = fp32_pi16; + /* local_nan = fp32_nan; local_ptr_cos_table = &fp32_cos_table[0][0]; local_ptr_sin_table = &fp32_sin_table[0][0]; @@ -699,8 +700,8 @@ inline two sin(const two &input) { local_2pi = fp64_2pi; local_pi2 = fp64_pi2; pointfive = 0.5; - /* local_pi16 = fp64_pi16; + /* local_nan = fp64_nan; local_ptr_cos_table = &fp64_cos_table[0][0]; local_ptr_sin_table = &fp64_sin_table[0][0]; @@ -728,8 +729,8 @@ inline two sin(const two &input) { 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); From 679b266675da0082f14bc29555d639cf9a400f97 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 10:57:46 -0500 Subject: [PATCH 059/137] Enabling more compiling code. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index f66a384..5131013 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -730,8 +730,8 @@ inline two sin(const two &input) { 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)); /* - t = sub(t, mul(local_pi16, q)); int k = static_cast(q); int abs_k = std::abs(k); From 3c0861e0324510305ebd06f0064476ad65071b0f Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 10:58:49 -0500 Subject: [PATCH 060/137] Enabling more compiling code. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 5131013..738a959 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -731,10 +731,10 @@ inline two sin(const two &input) { 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); +/* if (j < -2 || j > 2) { // std::error("LSV: cannot reduce modulo pi/2"); return local_nan; From 57e85cfd2f613a58a4324a89b32a60eed9a6d9c2 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 11:18:43 -0500 Subject: [PATCH 061/137] Enabling more compiling code. --- .../arithmetics/double-word-arithmetic.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 738a959..d33d605 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -680,17 +680,17 @@ inline two sin(const two &input) { two local_pi2; two local_pi16; /* - T pointfive, local_pi16, local_nan; + T pointfive, local_pi16; T* local_ptr_cos_table, local_ptr_sin_table; */ - T zeropointzero, pointfive; + T zeropointzero, pointfive, local_nan; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; local_pi2 = fp32_pi2; pointfive = 0.5f; local_pi16 = fp32_pi16; - /* local_nan = fp32_nan; + /* local_ptr_cos_table = &fp32_cos_table[0][0]; local_ptr_sin_table = &fp32_sin_table[0][0]; */ @@ -701,8 +701,8 @@ inline two sin(const two &input) { local_pi2 = fp64_pi2; pointfive = 0.5; local_pi16 = fp64_pi16; - /* local_nan = fp64_nan; + /* local_ptr_cos_table = &fp64_cos_table[0][0]; local_ptr_sin_table = &fp64_sin_table[0][0]; */ @@ -734,17 +734,17 @@ inline two sin(const two &input) { int k = static_cast(q); int abs_k = std::abs(k); -/* if (j < -2 || j > 2) { // std::error("LSV: cannot reduce modulo pi/2"); - return local_nan; + return two{local_nan, local_nan}; } if (abs_k > 4) { // std::error("LSV: cannot reduce modulo pi/16"); - return local_nan; + return two{local_nan, local_nan}; } +/* if (k == 0) { switch(j) { case 0: From 144c5f3563563b7eb5b168c4dc750efc624f65c2 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 11:23:31 -0500 Subject: [PATCH 062/137] Enabling more code. This requires a fix within sin_taylor() and cos_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index d33d605..78cbb93 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -744,7 +744,6 @@ inline two sin(const two &input) { return two{local_nan, local_nan}; } -/* if (k == 0) { switch(j) { case 0: @@ -752,12 +751,12 @@ inline two sin(const two &input) { case 1: return cos_taylor(t); case -1: - return -cos_taylor(t); + return two{-cos_taylor(t).h, -cos_taylor(t).l}; // Reference: QD / dd_inline.h default: - return -sin_taylor(t); + return two{-sin_taylor(t).h, -sin_taylor(t).l}; } } - +/* two u(local_ptr_cos_table[abs_k - 1][0], local_ptr_cos_table[abs_k - 1][1]); two v(local_ptr_sin_table[abs_k - 1][0], local_ptr_sin_table[abs_k - 1][1]); two sin_t, cos_t; From 83a8c5dc8c1c5fb100c5ae55910006763caee07a Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 11:27:02 -0500 Subject: [PATCH 063/137] Adding const type for pointer. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 78cbb93..0422791 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -494,7 +494,7 @@ template static two sin_taylor(const two &input) { T pointfive, local_eps; - T* local_ptr_inv_fact; + const T* local_ptr_inv_fact; T zeropointzero; if constexpr (std::is_same_v) { pointfive = 0.5f; @@ -748,12 +748,14 @@ inline two sin(const two &input) { 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}; +*/ } } /* From 372efdc4e51a735916f12cef2717197f6bcc0693 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 11:29:21 -0500 Subject: [PATCH 064/137] Starting to fix sin_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 0422791..90037d6 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -515,8 +515,11 @@ static two sin_taylor(const two &input) { two r, s, t, x; // TODO: make sure this is already supported in twofloat - if (input.is_zero()) { - return zeropointzero; + //if (input.is_zero()) { + if (input.eval() == zeropointzero) { + T hi, lo; + split(zeropointzero, hi, lo); + return two{hi, lo}; } int i = 0; From c2df2cb41bd233a83602b1082ad333fd4510bb99 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 11:47:17 -0500 Subject: [PATCH 065/137] Starting to fix sin_taylor(). --- .../arithmetics/double-word-arithmetic.hpp | 130 +++++++++--------- 1 file changed, 67 insertions(+), 63 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 90037d6..45f083c 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -356,7 +356,68 @@ namespace qd { err = b - (s - a); return s; } -} + + // TODO: make sure it has not been already implemented in twofloat + // Reference: QD / inline.h + /* Computes high word and low word of input */ + template + inline void split(T input, T &hi, T &lo) { + T temp; + + T local_qd_splitter; + T local_qd_split_thresh; + T local_split_factor; + T local_split_factor_2; + if constexpr (std::is_same_v) { + local_qd_splitter = fp32_qd_splitter; + local_qd_split_thresh = fp32_qd_split_thresh; + local_split_factor = fp32_split_factor; + local_split_factor_2 = fp32_split_factor_2; + } else if constexpr (std::is_same_v) { + local_qd_splitter = fp64_qd_splitter; + local_qd_split_thresh = fp64_qd_split_thresh; + local_split_factor = fp64_split_factor; + local_split_factor_2 = fp64_split_factor_2; + } else { + static_assert(sizeof(T) == 0, "LSV: other types not supported"); + } + + if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { + input *= local_split_factor; + temp = local_qd_splitter * input; + hi = temp - (temp - input); + lo = input - hi; + hi *= local_split_factor_2; + lo *= local_split_factor_2; + } else { + temp = local_qd_splitter * input; + hi = temp - (temp - input); + lo = input - hi; + } + } + + // Reference: QD / inline.h + /* Computes fl(input * input) and err(input * input) */ + template + inline T two_sqr(T input, T &err) { + + T twopointzero; + if constexpr (std::is_same_v) { + twopointzero = 2.0f; + } else if constexpr (std::is_same_v) { + twopointzero = 2.0; + } else { + static_assert(sizeof(T) == 0, "LSV: other types not supported"); + } + + T hi, lo; + T q = input * input; + qd::split(input, hi, lo); + err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; + return q; + } + +} // End namespace qd // Reference: QD / dd_inline.h /* Round to nearest integer */ @@ -405,65 +466,8 @@ inline T to_double(const two &input) { return input.h; } -// TODO: make sure it has not been already implemented in twofloat -// Reference: QD / inline.h -/* Computes high word and low word of input */ -template -inline void split(T input, T &hi, T &lo) { - T temp; - - T local_qd_splitter; - T local_qd_split_thresh; - T local_split_factor; - T local_split_factor_2; - if constexpr (std::is_same_v) { - local_qd_splitter = fp32_qd_splitter; - local_qd_split_thresh = fp32_qd_split_thresh; - local_split_factor = fp32_split_factor; - local_split_factor_2 = fp32_split_factor_2; - } else if constexpr (std::is_same_v) { - local_qd_splitter = fp64_qd_splitter; - local_qd_split_thresh = fp64_qd_split_thresh; - local_split_factor = fp64_split_factor; - local_split_factor_2 = fp64_split_factor_2; - } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); - } - - if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { - input *= local_split_factor; - temp = local_qd_splitter * input; - hi = temp - (temp - input); - lo = input - hi; - hi *= local_split_factor_2; - lo *= local_split_factor_2; - } else { - temp = local_qd_splitter * input; - hi = temp - (temp - input); - lo = input - hi; - } -} - -// Reference: QD / inline.h -/* Computes fl(input * input) and err(input * input) */ -template -inline T two_sqr(T input, T &err) { - T twopointzero; - if constexpr (std::is_same_v) { - twopointzero = 2.0f; - } else if constexpr (std::is_same_v) { - twopointzero = 2.0; - } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); - } - T hi, lo; - T q = input * input; - split(input, hi, lo); - err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; - return q; -} // Reference: QD / dd_inline.h template @@ -480,11 +484,11 @@ inline two sqr(const two &input) { T p1, p2; T s1, s2; - p1 = two_sqr(input.h, p2); + p1 = qd::two_sqr(input.h, p2); p2 += twopointzero * input.h * input.l; p2 += input.l * input.l; - s1 = quick_two_sum(p1, p2, s2); - return two(s1, s2); + s1 = qd::quick_two_sum(p1, p2, s2); + return two{s1, s2}; } // Reference: QD / dd_real.cpp @@ -518,7 +522,7 @@ static two sin_taylor(const two &input) { //if (input.is_zero()) { if (input.eval() == zeropointzero) { T hi, lo; - split(zeropointzero, hi, lo); + qd::split(zeropointzero, hi, lo); return two{hi, lo}; } @@ -718,7 +722,7 @@ inline two sin(const two &input) { //if (input.is_zero()) { if (input.eval() == zeropointzero) { T hi, lo; - split(zeropointzero, hi, lo); + qd::split(zeropointzero, hi, lo); return two{hi, lo}; } From b90c1d86f31a4e5530d0baa570670dc4e8ae01dd Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 13:36:16 -0500 Subject: [PATCH 066/137] Improvements for sin_taylor(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 45f083c..4608088 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -527,15 +527,18 @@ static two sin_taylor(const two &input) { } int i = 0; - x = -sqr(input); + 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])); + 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(to_double(t)) > thresh); return s; From f4cc3f28b3810939b4b659d6bf5e58d6d8e6fa98 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 13:42:43 -0500 Subject: [PATCH 067/137] Improvements for sin_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 4608088..5ec933a 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -534,8 +534,8 @@ static two sin_taylor(const two &input) { do { r = mul(r, x); + t = mul(r, two{(local_ptr_inv_fact+i)[0], (local_ptr_inv_fact+i)[1]}); /* - t = mul(r, two{(*local_ptr_inv_fact)[i][0], (*local_ptr_inv_fact)[i][1]}); s = add(s, t); i += 2; */ From fe761537ce91b69e0352063b9deac505cd41c742 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 13:45:28 -0500 Subject: [PATCH 068/137] Improvements for sin_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 5ec933a..aa81094 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -535,10 +535,8 @@ static two sin_taylor(const two &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); + s = add(s, t); i += 2; - */ } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); return s; From e7da4fd22bcb6672cd8818d71496f67a400e5fc3 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 13:50:01 -0500 Subject: [PATCH 069/137] Improvements in cos_taylor(). --- .../arithmetics/double-word-arithmetic.hpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index aa81094..7d92dbe 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -554,17 +554,20 @@ template static two cos_taylor(const two &input) { T pointfive, local_eps, onepointzero; - T* local_ptr_inv_fact; + const T* local_ptr_inv_fact; + T zeropointzero; if constexpr (std::is_same_v) { pointfive = 0.5f; local_eps = fp32_eps; onepointzero = 1.0f; local_ptr_inv_fact = &fp32_inv_fact[0][0]; + zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { pointfive = 0.5; local_eps = fp64_eps; onepointzero = 1.0; local_ptr_inv_fact = &fp64_inv_fact[0][0]; + zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "LSV: other types not supported"); } @@ -574,8 +577,11 @@ static two cos_taylor(const two &input) { two r, s, t, x; // TODO: make sure this is already supported in twofloat - if (input.is_zero()) { - return onepointzero; + //if (input.is_zero()) { + if (input.eval() == zeropointzero) { + T hi, lo; + qd::split(zeropointzero, hi, lo); + return two{hi, lo}; } x = -sqr(input); @@ -756,9 +762,9 @@ inline two sin(const two &input) { 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: From 044981ce7b341abfc5d2e456681a23a39a13a0cd Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 13:55:29 -0500 Subject: [PATCH 070/137] Improvements in cos_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7d92dbe..2448f68 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -584,15 +584,18 @@ static two cos_taylor(const two &input) { return two{hi, lo}; } - x = -sqr(input); + 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(to_double(t)) > thresh); return s; From 1116d845c1edc27a06482e7fec496d34be1e9936 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 13:56:58 -0500 Subject: [PATCH 071/137] Improvements in cos_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 2448f68..62586ee 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -590,8 +590,8 @@ static two cos_taylor(const two &input) { s = add(onepointzero, mul_pwr2(r, pointfive)); int i = 1; do { + r = mul(r, x); /* - 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; From 46b7581ccf55bbe45edfc70093087b99b484c530 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 13:58:49 -0500 Subject: [PATCH 072/137] Improvements in cos_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 62586ee..79ba821 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -591,8 +591,8 @@ static two cos_taylor(const two &input) { int i = 1; do { r = mul(r, x); + t = mul(r, two{(local_ptr_inv_fact+i)[0], (local_ptr_inv_fact+i)[1]}); /* - t = mul(r, two(*local_ptr_inv_fact[i][0], *local_ptr_inv_fact[i][1])); s = add(s, t); i += 2; */ From 35a0b62666669a8ab6fce98d590ddf7d47df9489 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 25 Jan 2024 14:00:06 -0500 Subject: [PATCH 073/137] Improvements in cos_taylor(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 79ba821..6f29da6 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -592,10 +592,8 @@ static two cos_taylor(const two &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); + s = add(s, t); i += 2; -*/ } while (i < n_inv_fact && std::abs(to_double(t)) > thresh); return s; From 8f8c43568be0144c8c1f9b2b7b2a774718970bec Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 11:37:28 -0500 Subject: [PATCH 074/137] More fixes in sin(). --- .../arithmetics/double-word-arithmetic.hpp | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 6f29da6..6b98efc 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -696,8 +696,10 @@ inline two sin(const two &input) { two local_pi16; /* T pointfive, local_pi16; - T* local_ptr_cos_table, local_ptr_sin_table; */ + const T* local_ptr_cos_table; + const T* local_ptr_sin_table; + T zeropointzero, pointfive, local_nan; if constexpr (std::is_same_v) { local_2pi = fp32_2pi; @@ -705,10 +707,8 @@ inline two sin(const two &input) { pointfive = 0.5f; local_pi16 = fp32_pi16; local_nan = fp32_nan; - /* local_ptr_cos_table = &fp32_cos_table[0][0]; local_ptr_sin_table = &fp32_sin_table[0][0]; - */ zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { @@ -717,10 +717,8 @@ inline two sin(const two &input) { pointfive = 0.5; local_pi16 = fp64_pi16; local_nan = fp64_nan; - /* local_ptr_cos_table = &fp64_cos_table[0][0]; local_ptr_sin_table = &fp64_sin_table[0][0]; - */ zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "LSV: other types not supported"); @@ -765,16 +763,15 @@ inline two sin(const two &input) { 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}; -*/ } } + + two u{(local_ptr_cos_table + (abs_k - 1))[0], (local_ptr_cos_table + (abs_k - 1))[1]}; /* - two u(local_ptr_cos_table[abs_k - 1][0], local_ptr_cos_table[abs_k - 1][1]); two v(local_ptr_sin_table[abs_k - 1][0], local_ptr_sin_table[abs_k - 1][1]); two sin_t, cos_t; sincos_taylor(t, sin_t, cos_t); From 02e3009e3d9ee10b86792926e3a110b5377d2669 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 11:40:19 -0500 Subject: [PATCH 075/137] More fixes in sin(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 6b98efc..a3b36c0 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -771,8 +771,8 @@ inline two sin(const two &input) { } two u{(local_ptr_cos_table + (abs_k - 1))[0], (local_ptr_cos_table + (abs_k - 1))[1]}; + two v{(local_ptr_sin_table + (abs_k - 1))[0], (local_ptr_sin_table + (abs_k - 1))[1]}; /* - two v(local_ptr_sin_table[abs_k - 1][0], local_ptr_sin_table[abs_k - 1][1]); two sin_t, cos_t; sincos_taylor(t, sin_t, cos_t); if (j == 0) { From 2a3f6c994f4681bb4741e331c528acf401ae4028 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 11:40:49 -0500 Subject: [PATCH 076/137] More fixes in sin(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index a3b36c0..69ed523 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -772,8 +772,9 @@ inline two sin(const two &input) { two u{(local_ptr_cos_table + (abs_k - 1))[0], (local_ptr_cos_table + (abs_k - 1))[1]}; two v{(local_ptr_sin_table + (abs_k - 1))[0], (local_ptr_sin_table + (abs_k - 1))[1]}; -/* + two sin_t, cos_t; +/* sincos_taylor(t, sin_t, cos_t); if (j == 0) { if (k > 0) { From 74c0ff497b51c1cf469d85043c8569bed59ae3d7 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 11:56:44 -0500 Subject: [PATCH 077/137] Partly fixing sincos_taylor(). --- .../arithmetics/double-word-arithmetic.hpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 69ed523..5d94ffa 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -677,14 +677,20 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { static_assert(sizeof(T) == 0, "LSV: other types not supported"); } - if (input.is_zero()) { - sin_a = zeropointzero; - cos_a = onepointzero; + // Reference: QD / dd_inline.h + // Assignment is performed according to operator== + // TODO: make sure this is already supported in twofloat + //if (input.is_zero()) { + 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))); + //cos_a = sqrt(sub(onepointzero, sqr(sin_a))); } // Reference: QD / dd_real.cpp @@ -774,8 +780,9 @@ inline two sin(const two &input) { two v{(local_ptr_sin_table + (abs_k - 1))[0], (local_ptr_sin_table + (abs_k - 1))[1]}; two sin_t, cos_t; -/* + sincos_taylor(t, sin_t, cos_t); +/* if (j == 0) { if (k > 0) { r = add(mul(u, sin_t), mul(v, cos_t)); From 84e86426a906fdd60123461477378db3b45b7831 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 12:01:38 -0500 Subject: [PATCH 078/137] Partly fixing sqrt(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 5d94ffa..0f5afe0 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -647,8 +647,10 @@ two sqrt(const two &input) { static_assert(sizeof(T) == 0, "LSV: other types not supported"); } - if (input.is_zero()) { - return zeropointzero; + // TODO: make sure this is already supported in twofloat + //if (input.is_zero()) { + if (input.eval() == zeropointzero) { + return two{zeropointzero, zeropointzero}; } if (input.is_negative()) { @@ -690,7 +692,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { } sin_a = sin_taylor(input); - //cos_a = sqrt(sub(onepointzero, sqr(sin_a))); + cos_a = sqrt(sub(onepointzero, sqr(sin_a))); } // Reference: QD / dd_real.cpp From 52795d4a93e9029715ce32b69b2ef2156cdab6e9 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 12:09:37 -0500 Subject: [PATCH 079/137] Fixes in sqrt(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 0f5afe0..853d467 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -653,13 +653,14 @@ two sqrt(const two &input) { return two{zeropointzero, zeropointzero}; } - if (input.is_negative()) { + // if (input.is_negative()) { + if (input.h < zeropointzero) { static_assert(sizeof(T) == 0, "LSV: other types not supported"); - return local_nan; + return two{local_nan, local_nan}; } T x = onepointzero / std::sqrt(input.h); - T ax = mul(input.h, x); + T ax = input.h * x; two temp = sub(input, sqr(ax)); return dd_add(ax, temp.h * x * pointfive); } From e9e2159df6b70c69c4101d7ac8b71ace7a84d28f Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 13:11:55 -0500 Subject: [PATCH 080/137] More fixes in sqrt(). --- .../arithmetics/double-word-arithmetic.hpp | 56 +++++++++++-------- 1 file changed, 33 insertions(+), 23 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 853d467..7c5efb2 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -417,8 +417,39 @@ namespace qd { return q; } + // Reference QD / inline.h + // TODO: make sure this has not been already implemented + 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; + } + } // End namespace qd +namespace dd_real { + + // Reference: QD / dd_inline.h + template + inline two sqr(T input) { + T p1, p2; + p1 = qd::two_sqr(input, p2); + return two{p1, p2}; + } + + // Reference: QD / dd_inline.h + // TODO: make sure this has not been already implemented + template + inline two add(T a, T b) { + T s, e; + s = two_sum(a, b, e); + return two (s, e); + } + +} // End namespace dd_real + // Reference: QD / dd_inline.h /* Round to nearest integer */ template @@ -466,9 +497,6 @@ inline T to_double(const two &input) { return input.h; } - - - // Reference: QD / dd_inline.h template inline two sqr(const two &input) { @@ -599,24 +627,6 @@ static two cos_taylor(const two &input) { return s; } -// Reference QD / inline.h -// TODO: make sure this has not been already implemented -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 -template -inline two dd_add(T a, T b) { - T s, e; - s = two_sum(a, b, e); - return two (s, e); -} - // Reference: QD / dd_real.cpp /* Computes the square root of the double-double number dd. NOTE: dd must be a non-negative number. */ @@ -661,8 +671,8 @@ two sqrt(const two &input) { T x = onepointzero / std::sqrt(input.h); T ax = input.h * x; - two temp = sub(input, sqr(ax)); - return dd_add(ax, temp.h * x * pointfive); + two temp = sub(input, dd_real::sqr(ax)); + return dd_real::add(ax, temp.h * x * pointfive); } // Reference: QD / dd_real.cpp From 32eceefc2cc6d6833d2ecbe25efd41228efbf5f9 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 13:16:53 -0500 Subject: [PATCH 081/137] Fixing add(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7c5efb2..1eb77fb 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -444,7 +444,7 @@ namespace dd_real { template inline two add(T a, T b) { T s, e; - s = two_sum(a, b, e); + s = qd::two_sum(a, b, e); return two (s, e); } @@ -665,7 +665,8 @@ two sqrt(const two &input) { // if (input.is_negative()) { if (input.h < zeropointzero) { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + // TODO: replace assert with message printing + // static_assert(sizeof(T) == 0, "LSV: other types not supported"); return two{local_nan, local_nan}; } From 077eadc2f54403e1d6c04d8ab785186fbd5da71e Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 13:19:55 -0500 Subject: [PATCH 082/137] Fixing sin(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 1eb77fb..753c7c5 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -796,14 +796,14 @@ inline two sin(const two &input) { two sin_t, cos_t; sincos_taylor(t, sin_t, cos_t); -/* + if (j == 0) { if (k > 0) { - r = add(mul(u, sin_t), mul(v, cos_t)); + r = add(mul(u, sin_t), mul(v, cos_t)); } else { - r = sub(mul(u, sin_t), mul(v, cos_t)); + r = sub(mul(u, sin_t), mul(v, cos_t)); } - } else if (j == 1) { + } /*else if (j == 1) { if (k > 0) { r = sub(mul(u, cos_t), mul(v, sin_t)); } else { From 2f2c3469a0fed30523694b1d099649e3c1ce79c2 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 13:21:11 -0500 Subject: [PATCH 083/137] More fixes in sin(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 753c7c5..6bd3612 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -803,13 +803,13 @@ inline two sin(const two &input) { } else { r = sub(mul(u, sin_t), mul(v, cos_t)); } - } /*else if (j == 1) { + } else if (j == 1) { if (k > 0) { - r = sub(mul(u, cos_t), mul(v, sin_t)); + r = sub(mul(u, cos_t), mul(v, sin_t)); } else { - r = add(mul(u, cos_t), mul(v, sin_t)); + r = add(mul(u, cos_t), mul(v, sin_t)); } - } else if (j == -1) { + } /*else if (j == -1) { if (k > 0) { r = sub(mul(v, sin_t), mul(u, cos_t)); } else if (k < 0) { From 73c39764aa605c40ae28b3914545e647588916dc Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 13:24:13 -0500 Subject: [PATCH 084/137] More fixes in sin(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 6bd3612..761121d 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -809,13 +809,15 @@ inline two sin(const two &input) { } else { r = add(mul(u, cos_t), mul(v, sin_t)); } - } /*else if (j == -1) { + } else if (j == -1) { if (k > 0) { - r = sub(mul(v, sin_t), mul(u, cos_t)); + r = sub(mul(v, sin_t), mul(u, cos_t)); } else if (k < 0) { - r = sub(-mul(u, cos_t), mul(v, sin_t)); + two temp = mul(u, cos_t); + two neg_temp = two{-temp.h, -temp.l}; + r = sub(neg_temp, mul(v, sin_t)); } - } else { + } /*else { if (k > 0) { r = sub(-mul(u, sin_t), mul(v, cos_t)); } else { From 1c907a5608b7d95bd50a73d1ac845b86ee69ecae Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 13:26:30 -0500 Subject: [PATCH 085/137] More fixes in sin(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 761121d..d7baa72 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -817,16 +817,17 @@ inline two sin(const two &input) { two neg_temp = two{-temp.h, -temp.l}; r = sub(neg_temp, mul(v, sin_t)); } - } /*else { + } else { if (k > 0) { - r = sub(-mul(u, sin_t), mul(v, cos_t)); + 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)); + r = sub(mul(v, cos_t), mul(u, sin_t)); } } return r; -*/ } } // namespace doubleword } // namespace twofloat \ No newline at end of file From e987e8986258fe3b8a3455bc3c145a12f3774848 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 14:05:46 -0500 Subject: [PATCH 086/137] Completing simple SinTest(). --- test/double-word-arithmetic.test.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index 31576c1..3fa1a48 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -255,12 +255,16 @@ TEST(DoubleWordArithmetic, DivFPTest) { divFPTest(); } TEST(DoubleWordArithmetic, DivFPFMATest) { divFPTest(); } -TEST(DoubleWordArithmetic, SinTest) {} +TEST(DoubleWordArithmetic, SinTest) { float x = 10; float sinx = std::sin(x); two two_x = {x, 0.0f}; two sin_two_x = doubleword::sin(two_x); + + EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); +} + } // namespace test } // namespace doubleword } // namespace twofloat \ No newline at end of file From 466aaa7a08b5257d7cb29379fff34b0f7f791dd2 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 14:12:42 -0500 Subject: [PATCH 087/137] Additional hints for inspecting testing. --- test/double-word-arithmetic.test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index 3fa1a48..884e7de 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -262,6 +262,7 @@ TEST(DoubleWordArithmetic, SinTest) { 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_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); } From b3ce7efc4b954d1a84c6388933168d0168212869 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 14:52:49 -0500 Subject: [PATCH 088/137] Adding SinTest2 for specifically testing double-precision sin(). --- test/double-word-arithmetic.test.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index 884e7de..727b7c2 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -266,6 +266,17 @@ TEST(DoubleWordArithmetic, SinTest) { EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); } +TEST(DoubleWordArithmetic, SinTest2) { + double x = 10; + 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_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); +} + } // namespace test } // namespace doubleword } // namespace twofloat \ No newline at end of file From 7269ef542cb97beab79fef7e126c09f81896bf64 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 14:59:31 -0500 Subject: [PATCH 089/137] Parameterizing returning type of nint(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index d7baa72..5a77157 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -487,7 +487,7 @@ inline two nint(const two &input) { } } - return two{hi, lo}; + return two{hi, lo}; } // Reference: QD / dd_inline.h From abec02ffbb03cc2d1bde7ef01ee2886ec327d8a3 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 15:02:06 -0500 Subject: [PATCH 090/137] Minor cleanup. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 5a77157..729b6de 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -714,9 +714,6 @@ inline two sin(const two &input) { two local_2pi; two local_pi2; two local_pi16; -/* - T pointfive, local_pi16; -*/ const T* local_ptr_cos_table; const T* local_ptr_sin_table; From 1aa48f557420bf7279bf6f8c8f500ca94b522084 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 16:07:38 -0500 Subject: [PATCH 091/137] Cleaning up. --- .../arithmetics/double-word-arithmetic.hpp | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 729b6de..0e3fb4f 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -716,33 +716,30 @@ inline two sin(const two &input) { two local_pi16; const T* local_ptr_cos_table; const T* local_ptr_sin_table; - T zeropointzero, pointfive, local_nan; + if constexpr (std::is_same_v) { - local_2pi = fp32_2pi; - local_pi2 = fp32_pi2; - pointfive = 0.5f; + local_2pi = fp32_2pi; + local_pi2 = fp32_pi2; local_pi16 = fp32_pi16; - local_nan = fp32_nan; local_ptr_cos_table = &fp32_cos_table[0][0]; local_ptr_sin_table = &fp32_sin_table[0][0]; zeropointzero = 0.0f; - + pointfive = 0.5f; + local_nan = fp32_nan; } else if constexpr (std::is_same_v) { - local_2pi = fp64_2pi; - local_pi2 = fp64_pi2; - pointfive = 0.5; + local_2pi = fp64_2pi; + local_pi2 = fp64_pi2; local_pi16 = fp64_pi16; - local_nan = fp64_nan; local_ptr_cos_table = &fp64_cos_table[0][0]; local_ptr_sin_table = &fp64_sin_table[0][0]; zeropointzero = 0.0; + pointfive = 0.5; + local_nan = fp64_nan; } else { static_assert(sizeof(T) == 0, "LSV: other types not supported"); } - // TODO: make sure this is already supported in twofloat - //if (input.is_zero()) { if (input.eval() == zeropointzero) { T hi, lo; qd::split(zeropointzero, hi, lo); From 8353d1ceaf338143c7f49887fc5db597a7aa588a Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 16:41:12 -0500 Subject: [PATCH 092/137] Fixing return value. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 0e3fb4f..892d215 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -546,8 +546,6 @@ static two sin_taylor(const two &input) { two r, s, t, x; - // TODO: make sure this is already supported in twofloat - //if (input.is_zero()) { if (input.eval() == zeropointzero) { T hi, lo; qd::split(zeropointzero, hi, lo); @@ -604,11 +602,9 @@ static two cos_taylor(const two &input) { two r, s, t, x; - // TODO: make sure this is already supported in twofloat - //if (input.is_zero()) { if (input.eval() == zeropointzero) { T hi, lo; - qd::split(zeropointzero, hi, lo); + qd::split(onepointzero, hi, lo); return two{hi, lo}; } From ce7b4d5650f1becce4477dd186b7a998ee8552fc Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 30 Jan 2024 16:49:38 -0500 Subject: [PATCH 093/137] Uniform return value. --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 892d215..7f73cb3 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -653,10 +653,10 @@ two sqrt(const two &input) { static_assert(sizeof(T) == 0, "LSV: other types not supported"); } - // TODO: make sure this is already supported in twofloat - //if (input.is_zero()) { if (input.eval() == zeropointzero) { - return two{zeropointzero, zeropointzero}; + T hi, lo; + qd::split(zeropointzero, hi, lo); + return two{hi, lo}; } // if (input.is_negative()) { @@ -689,8 +689,6 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { // Reference: QD / dd_inline.h // Assignment is performed according to operator== - // TODO: make sure this is already supported in twofloat - //if (input.is_zero()) { if (input.eval() == zeropointzero) { sin_a.h = zeropointzero; sin_a.l = zeropointzero; From 34dd5f2b6d7a6f765866bdcbff18778240da98ae Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 1 Feb 2024 16:04:40 -0500 Subject: [PATCH 094/137] Adding print features to debug sin(). This changes are mirrowed from those in QD / sin(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 6 ++++++ include/libtwofloat/twofloat.hpp | 1 + test/double-word-arithmetic.test.cpp | 4 ++-- 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7f73cb3..5292b7f 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -744,6 +744,12 @@ inline two sin(const two &input) { two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div (accurate + fma) 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 //TODO: original type is double, here it is templated 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/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index 727b7c2..0293128 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -256,7 +256,7 @@ TEST(DoubleWordArithmetic, DivFPTest) { divFPTest(); } TEST(DoubleWordArithmetic, DivFPFMATest) { divFPTest(); } TEST(DoubleWordArithmetic, SinTest) { - float x = 10; + float x = 1; float sinx = std::sin(x); two two_x = {x, 0.0f}; @@ -267,7 +267,7 @@ TEST(DoubleWordArithmetic, SinTest) { } TEST(DoubleWordArithmetic, SinTest2) { - double x = 10; + double x = 1; double sinx = std::sin(x); two two_x = {x, 0.0}; From ef196cf84a32c765811b32f818c10bb4c6437abf Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 1 Feb 2024 16:17:06 -0500 Subject: [PATCH 095/137] More printing. --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 5292b7f..578d55a 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -761,6 +761,14 @@ inline two sin(const two &input) { 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; + if (j < -2 || j > 2) { // std::error("LSV: cannot reduce modulo pi/2"); return two{local_nan, local_nan}; From ca8274ce58fe9b855a8e16747f55b9a90e9986e1 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 1 Feb 2024 16:32:09 -0500 Subject: [PATCH 096/137] Fixing error messages. --- .../arithmetics/double-word-arithmetic.hpp | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 578d55a..9908493 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -338,7 +338,7 @@ namespace qd { } else if constexpr (std::is_same_v) { pointfive = 0.5; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } if(input = std::floor(input)) { @@ -379,7 +379,7 @@ namespace qd { local_split_factor = fp64_split_factor; local_split_factor_2 = fp64_split_factor_2; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { @@ -407,7 +407,7 @@ namespace qd { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } T hi, lo; @@ -467,7 +467,7 @@ inline two nint(const two &input) { pointfive = 0.5; onepointzero = 1.0; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } if(hi == input.h) { @@ -507,7 +507,7 @@ inline two sqr(const two &input) { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } T p1, p2; @@ -539,7 +539,7 @@ static two sin_taylor(const two &input) { local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } const T thresh = pointfive * std::abs(to_double(input)) * local_eps; @@ -595,7 +595,7 @@ static two cos_taylor(const two &input) { local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } const T thresh = pointfive * local_eps; @@ -650,7 +650,7 @@ two sqrt(const two &input) { pointfive = 0.5; local_nan = fp64_nan; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } if (input.eval() == zeropointzero) { @@ -661,8 +661,7 @@ two sqrt(const two &input) { // if (input.is_negative()) { if (input.h < zeropointzero) { - // TODO: replace assert with message printing - // static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Negative argument.\n"; return two{local_nan, local_nan}; } @@ -684,7 +683,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { zeropointzero = 0.0; onepointzero = 1.0; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } // Reference: QD / dd_inline.h @@ -731,7 +730,7 @@ inline two sin(const two &input) { pointfive = 0.5; local_nan = fp64_nan; } else { - static_assert(sizeof(T) == 0, "LSV: other types not supported"); + std::cerr << "Other types not supported\n"; } if (input.eval() == zeropointzero) { @@ -770,12 +769,12 @@ inline two sin(const two &input) { std::cout << "... LSV ends" << std::endl; if (j < -2 || j > 2) { -// std::error("LSV: cannot reduce modulo pi/2"); + std::cerr << "ERROR: cannot reduce modulo pi/2\n"; return two{local_nan, local_nan}; } if (abs_k > 4) { -// std::error("LSV: cannot reduce modulo pi/16"); + std::cerr << "ERROR: cannot reduce modulo pi/16\n"; return two{local_nan, local_nan}; } From 23ab1e3950192b13ee1eb85c45747d704b888fb4 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 1 Feb 2024 16:38:33 -0500 Subject: [PATCH 097/137] More printing. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 9908493..0bb1846 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -798,6 +798,11 @@ inline two sin(const two &input) { 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)); From 82a33c8577d68a2dba158411ab4a5e24c1be0ba8 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 1 Feb 2024 16:43:38 -0500 Subject: [PATCH 098/137] More printing. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 0bb1846..a3ca863 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -833,6 +833,10 @@ inline two sin(const two &input) { } } + std::cout << "LSV starts ..." << std::endl; + std::cout << "r = " << r.eval() << std::endl; + std::cout << "... LSV ends" << std::endl; + return r; } } // namespace doubleword From 865bd26b9aee6330b9384b49b598058ec85b8c4a Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Thu, 1 Feb 2024 16:51:25 -0500 Subject: [PATCH 099/137] More printing. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index a3ca863..4897667 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -814,6 +814,13 @@ inline two sin(const two &input) { 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) { From 61bdaf45f113a8e3c6885e5a4bc5b5a9df081d55 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Fri, 2 Feb 2024 11:34:32 -0500 Subject: [PATCH 100/137] Fixing implementation of sin(). --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 4897667..0f20979 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -791,8 +791,10 @@ inline two sin(const two &input) { } } - two u{(local_ptr_cos_table + (abs_k - 1))[0], (local_ptr_cos_table + (abs_k - 1))[1]}; - two v{(local_ptr_sin_table + (abs_k - 1))[0], (local_ptr_sin_table + (abs_k - 1))[1]}; + // 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; From 33636c7dbb9d4c6a77efd20f9ebf283af9000a8a Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Fri, 2 Feb 2024 11:44:18 -0500 Subject: [PATCH 101/137] REplacing EXPECT_NEAR with EXPECT_FLOAT_EQ and EXPECt_DOUBLE_EQ. --- test/double-word-arithmetic.test.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index 0293128..fb6597c 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -263,7 +263,8 @@ TEST(DoubleWordArithmetic, SinTest) { two sin_two_x = doubleword::sin(two_x); // Check the expected and actual values within /build/Testing/Temporary/LastTest.log - EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); + EXPECT_FLOAT_EQ(sinx, sin_two_x.eval()); + //EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); } TEST(DoubleWordArithmetic, SinTest2) { @@ -274,7 +275,8 @@ TEST(DoubleWordArithmetic, SinTest2) { two sin_two_x = doubleword::sin(two_x); // Check the expected and actual values within /build/Testing/Temporary/LastTest.log - EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); + EXPECT_DOUBLE_EQ(sinx, sin_two_x.eval()); + //EXPECT_NEAR(sinx, sin_two_x.eval(), sinx * std::numeric_limits>::epsilon()); } } // namespace test From 0a1713801aaa074e7e9ccc2e2eb90e93d8f1cdc0 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Fri, 2 Feb 2024 11:47:40 -0500 Subject: [PATCH 102/137] Commenting out some helper code. --- .../arithmetics/double-word-arithmetic.hpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 0f20979..2d6d5a7 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -743,11 +743,13 @@ inline two sin(const two &input) { two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div (accurate + fma) 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 @@ -760,6 +762,7 @@ inline two sin(const two &input) { 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; @@ -767,6 +770,7 @@ inline two sin(const two &input) { std::cout << "k = " << k << std::endl; std::cout << "abs_k = " << abs_k << std::endl; std::cout << "... LSV ends" << std::endl; + */ if (j < -2 || j > 2) { std::cerr << "ERROR: cannot reduce modulo pi/2\n"; @@ -800,10 +804,12 @@ inline two sin(const two &input) { 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) { @@ -816,6 +822,7 @@ inline two sin(const two &input) { 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; @@ -823,6 +830,7 @@ inline two sin(const two &input) { 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) { @@ -841,11 +849,11 @@ inline two sin(const two &input) { 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 From e5028370ba89844bc70b9844da805092184d5ec9 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 10:27:41 -0500 Subject: [PATCH 103/137] Replacing std::cerr with static_assert in some cases. --- .../arithmetics/double-word-arithmetic.hpp | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 2d6d5a7..c500072 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -338,7 +338,7 @@ namespace qd { } else if constexpr (std::is_same_v) { pointfive = 0.5; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } if(input = std::floor(input)) { @@ -379,7 +379,7 @@ namespace qd { local_split_factor = fp64_split_factor; local_split_factor_2 = fp64_split_factor_2; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { @@ -407,7 +407,7 @@ namespace qd { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } T hi, lo; @@ -467,7 +467,7 @@ inline two nint(const two &input) { pointfive = 0.5; onepointzero = 1.0; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } if(hi == input.h) { @@ -507,7 +507,7 @@ inline two sqr(const two &input) { } else if constexpr (std::is_same_v) { twopointzero = 2.0; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } T p1, p2; @@ -539,7 +539,7 @@ static two sin_taylor(const two &input) { local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } const T thresh = pointfive * std::abs(to_double(input)) * local_eps; @@ -595,7 +595,7 @@ static two cos_taylor(const two &input) { local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } const T thresh = pointfive * local_eps; @@ -650,7 +650,7 @@ two sqrt(const two &input) { pointfive = 0.5; local_nan = fp64_nan; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } if (input.eval() == zeropointzero) { @@ -683,7 +683,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { zeropointzero = 0.0; onepointzero = 1.0; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } // Reference: QD / dd_inline.h @@ -730,7 +730,7 @@ inline two sin(const two &input) { pointfive = 0.5; local_nan = fp64_nan; } else { - std::cerr << "Other types not supported\n"; + static_assert(sizeof(T) == 0, "Other types not supported"); } if (input.eval() == zeropointzero) { From 7b7fb912495647d0ec9dbde57b6df6d2124b0035 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 10:55:06 -0500 Subject: [PATCH 104/137] Removing remaining std::cerr instances. Instead adding warning comments on corresponding code block. --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index c500072..7cfcd0c 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -659,9 +659,9 @@ two sqrt(const two &input) { return two{hi, lo}; } - // if (input.is_negative()) { + // ERROR: Negative argument + // if (input.is_negative()) if (input.h < zeropointzero) { - std::cerr << "Negative argument.\n"; return two{local_nan, local_nan}; } @@ -772,13 +772,13 @@ inline two sin(const two &input) { std::cout << "... LSV ends" << std::endl; */ + // ERROR: cannot reduce modulo pi/2\n if (j < -2 || j > 2) { - std::cerr << "ERROR: cannot reduce modulo pi/2\n"; return two{local_nan, local_nan}; } + // ERROR: cannot reduce modulo pi/16 if (abs_k > 4) { - std::cerr << "ERROR: cannot reduce modulo pi/16\n"; return two{local_nan, local_nan}; } From dac3be13d000e5bbccd05f3f622a26d369b9aa19 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 11:41:52 -0500 Subject: [PATCH 105/137] Replacing calls to qd::quick_two_sum() with already implemented twofloat::algorithms::FastTwoSum(). --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7cfcd0c..99ee244 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -475,7 +475,10 @@ inline two nint(const two &input) { lo = qd::nint(input.l); /* Renormalize. This is needed if h = some integer, l = 1/2. */ - hi = qd::quick_two_sum(hi, lo, lo); + //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 */ @@ -515,7 +518,10 @@ inline two sqr(const two &input) { p1 = qd::two_sqr(input.h, p2); p2 += twopointzero * input.h * input.l; p2 += input.l * input.l; - s1 = qd::quick_two_sum(p1, p2, s2); + //s1 = qd::quick_two_sum(p1, p2, s2); + two temp = twofloat::algorithms::FastTwoSum(p1, p2); + s1 = temp.h; + s2 = temp.l; return two{s1, s2}; } From 322ac49c838d38916b61f7c145e11f334cc36fec Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 11:44:05 -0500 Subject: [PATCH 106/137] Removing implementation of qd::quick_two_sum(). This is not needed as it was already implemented similarly in twofloat. --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 99ee244..1189eaa 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -347,16 +347,6 @@ namespace qd { return std::floor(input + pointfive); } - // TODO: already implemented in twofloat? - // Reference: QD / inline.h - /* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */ - template - inline T quick_two_sum(T a, T b, T &err) { - T s = a + b; - err = b - (s - a); - return s; - } - // TODO: make sure it has not been already implemented in twofloat // Reference: QD / inline.h /* Computes high word and low word of input */ From e96936dfe28f3fcb674cbc4aa765a8cc0ee1a744 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 12:00:46 -0500 Subject: [PATCH 107/137] Replacing calls to qd::split() with already-implemented twofloat::algorithms::Split(). --- .../arithmetics/double-word-arithmetic.hpp | 25 +++++++++++++++---- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 1189eaa..65d1cd6 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -402,7 +402,10 @@ namespace qd { T hi, lo; T q = input * input; - qd::split(input, hi, lo); + //qd::split(input, hi, lo); + two temp = twofloat::algorithms::Split(input); + hi = temp.h; + lo = temp.l; err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; return q; } @@ -544,7 +547,10 @@ static two sin_taylor(const two &input) { if (input.eval() == zeropointzero) { T hi, lo; - qd::split(zeropointzero, hi, lo); + //qd::split(zeropointzero, hi, lo); + two temp = twofloat::algorithms::Split(zeropointzero); + hi = temp.h; + lo = temp.l; return two{hi, lo}; } @@ -600,7 +606,10 @@ static two cos_taylor(const two &input) { if (input.eval() == zeropointzero) { T hi, lo; - qd::split(onepointzero, hi, lo); + //qd::split(onepointzero, hi, lo); + two temp = twofloat::algorithms::Split(onepointzero); + hi = temp.h; + lo = temp.l; return two{hi, lo}; } @@ -651,7 +660,10 @@ two sqrt(const two &input) { if (input.eval() == zeropointzero) { T hi, lo; - qd::split(zeropointzero, hi, lo); + //qd::split(zeropointzero, hi, lo); + two temp = twofloat::algorithms::Split(zeropointzero); + hi = temp.h; + lo = temp.l; return two{hi, lo}; } @@ -731,7 +743,10 @@ inline two sin(const two &input) { if (input.eval() == zeropointzero) { T hi, lo; - qd::split(zeropointzero, hi, lo); + //qd::split(zeropointzero, hi, lo); + two temp = twofloat::algorithms::Split(zeropointzero); + hi = temp.h; + lo = temp.l; return two{hi, lo}; } From 5301ded2cdd5067228c9e22ccc710af206e0224b Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 12:03:32 -0500 Subject: [PATCH 108/137] Removing definition of qd::split(). This is not needed anymore as same functionality was already implemented in twofloat. --- .../arithmetics/double-word-arithmetic.hpp | 39 ------------------- 1 file changed, 39 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 65d1cd6..598ddfa 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -347,45 +347,6 @@ namespace qd { return std::floor(input + pointfive); } - // TODO: make sure it has not been already implemented in twofloat - // Reference: QD / inline.h - /* Computes high word and low word of input */ - template - inline void split(T input, T &hi, T &lo) { - T temp; - - T local_qd_splitter; - T local_qd_split_thresh; - T local_split_factor; - T local_split_factor_2; - if constexpr (std::is_same_v) { - local_qd_splitter = fp32_qd_splitter; - local_qd_split_thresh = fp32_qd_split_thresh; - local_split_factor = fp32_split_factor; - local_split_factor_2 = fp32_split_factor_2; - } else if constexpr (std::is_same_v) { - local_qd_splitter = fp64_qd_splitter; - local_qd_split_thresh = fp64_qd_split_thresh; - local_split_factor = fp64_split_factor; - local_split_factor_2 = fp64_split_factor_2; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } - - if(input > local_qd_split_thresh || input < -local_qd_split_thresh) { - input *= local_split_factor; - temp = local_qd_splitter * input; - hi = temp - (temp - input); - lo = input - hi; - hi *= local_split_factor_2; - lo *= local_split_factor_2; - } else { - temp = local_qd_splitter * input; - hi = temp - (temp - input); - lo = input - hi; - } - } - // Reference: QD / inline.h /* Computes fl(input * input) and err(input * input) */ template From 3dd021e300d4b6b1a5d6067f032f4bb5492bb2e9 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 12:18:22 -0500 Subject: [PATCH 109/137] Removing definition and replacing invocations of to_double(A) with A.eval(). --- .../arithmetics/double-word-arithmetic.hpp | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 598ddfa..8e702d2 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -447,13 +447,6 @@ inline two nint(const two &input) { return two{hi, lo}; } -// Reference: QD / dd_inline.h -/* Cast to double. */ -template -inline T to_double(const two &input) { - return input.h; -} - // Reference: QD / dd_inline.h template inline two sqr(const two &input) { @@ -502,7 +495,8 @@ static two sin_taylor(const two &input) { static_assert(sizeof(T) == 0, "Other types not supported"); } - const T thresh = pointfive * std::abs(to_double(input)) * local_eps; + //const T thresh = pointfive * std::abs(to_double(input)) * local_eps; + const T thresh = pointfive * std::abs(input.eval()) * local_eps; two r, s, t, x; @@ -526,7 +520,8 @@ static two sin_taylor(const two &input) { 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(to_double(t)) > thresh); + //} while (i < n_inv_fact && std::abs(to_double(t)) > thresh); + } while (i < n_inv_fact && std::abs(t.eval()) > thresh); return s; } @@ -584,7 +579,8 @@ static two cos_taylor(const two &input) { 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(to_double(t)) > thresh); + //} while (i < n_inv_fact && std::abs(to_double(t)) > thresh); + } while (i < n_inv_fact && std::abs(t.eval()) > thresh); return s; } From c29ad7b0c8f9ec2ace2c803fa6352c65e9f8d75d Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 12:55:34 -0500 Subject: [PATCH 110/137] Commenting out unnecessarily compute-intensive assignments to 0.0. --- .../arithmetics/double-word-arithmetic.hpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 8e702d2..b4a8e88 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -501,12 +501,15 @@ static two sin_taylor(const two &input) { two r, s, t, x; if (input.eval() == zeropointzero) { + /* T hi, lo; //qd::split(zeropointzero, hi, lo); two temp = twofloat::algorithms::Split(zeropointzero); hi = temp.h; lo = temp.l; return two{hi, lo}; + */ + return two{zeropointzero, zeropointzero}; } int i = 0; @@ -561,12 +564,15 @@ static two cos_taylor(const two &input) { two r, s, t, x; if (input.eval() == zeropointzero) { + /* T hi, lo; //qd::split(onepointzero, hi, lo); two temp = twofloat::algorithms::Split(onepointzero); hi = temp.h; lo = temp.l; return two{hi, lo}; + */ + return two{onepointzero, zeropointzero}; } two temp = sqr(input); @@ -616,12 +622,15 @@ two sqrt(const two &input) { } if (input.eval() == zeropointzero) { + /* T hi, lo; //qd::split(zeropointzero, hi, lo); two temp = twofloat::algorithms::Split(zeropointzero); hi = temp.h; lo = temp.l; return two{hi, lo}; + */ + return two{zeropointzero, zeropointzero}; } // ERROR: Negative argument @@ -699,12 +708,15 @@ inline two sin(const two &input) { } if (input.eval() == zeropointzero) { + /* T hi, lo; //qd::split(zeropointzero, hi, lo); two temp = twofloat::algorithms::Split(zeropointzero); hi = temp.h; lo = temp.l; return two{hi, lo}; + */ + return two{zeropointzero, zeropointzero}; } // Approximately reducing modulo 2*pi From 28d95dd609f0dd902331a47db5b276b7b50c0f51 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Tue, 6 Feb 2024 12:58:00 -0500 Subject: [PATCH 111/137] Removing commented-out code from previous commit. --- .../arithmetics/double-word-arithmetic.hpp | 32 ------------------- 1 file changed, 32 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index b4a8e88..7f31ee6 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -501,14 +501,6 @@ static two sin_taylor(const two &input) { two r, s, t, x; if (input.eval() == zeropointzero) { - /* - T hi, lo; - //qd::split(zeropointzero, hi, lo); - two temp = twofloat::algorithms::Split(zeropointzero); - hi = temp.h; - lo = temp.l; - return two{hi, lo}; - */ return two{zeropointzero, zeropointzero}; } @@ -564,14 +556,6 @@ static two cos_taylor(const two &input) { two r, s, t, x; if (input.eval() == zeropointzero) { - /* - T hi, lo; - //qd::split(onepointzero, hi, lo); - two temp = twofloat::algorithms::Split(onepointzero); - hi = temp.h; - lo = temp.l; - return two{hi, lo}; - */ return two{onepointzero, zeropointzero}; } @@ -622,14 +606,6 @@ two sqrt(const two &input) { } if (input.eval() == zeropointzero) { - /* - T hi, lo; - //qd::split(zeropointzero, hi, lo); - two temp = twofloat::algorithms::Split(zeropointzero); - hi = temp.h; - lo = temp.l; - return two{hi, lo}; - */ return two{zeropointzero, zeropointzero}; } @@ -708,14 +684,6 @@ inline two sin(const two &input) { } if (input.eval() == zeropointzero) { - /* - T hi, lo; - //qd::split(zeropointzero, hi, lo); - two temp = twofloat::algorithms::Split(zeropointzero); - hi = temp.h; - lo = temp.l; - return two{hi, lo}; - */ return two{zeropointzero, zeropointzero}; } From 97a998045d061f56472ae3d1bb090c1cc6f05013 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 07:33:56 -0500 Subject: [PATCH 112/137] Replacing fp32/fp64_sin_table and cos_table with templated struct members. --- .../arithmetics/double-word-arithmetic.hpp | 49 +++++++------------ 1 file changed, 19 insertions(+), 30 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7f31ee6..46c5dcb 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -296,32 +296,22 @@ static const double fp64_inv_fact[n_inv_fact][2] = { { 2.81145725434552060e-15, 1.65088427308614326e-31} }; -// Reference: QD / dd_real.cpp -/* Table of sin(k * pi/16) and cos(k * pi/16). */ -static const float fp32_cos_table [4][2] = { - {9.8078528e-01, 1.8546939e-17}, - {9.2387953e-01, 1.7645047e-17}, - {8.3146961e-01, 1.4073856e-18}, - {7.0710678e-01, -4.8336466e-17} -}; -static const float fp32_sin_table [4][2] = { - {1.9509032e-01, -7.9910790e-18}, - {3.8268343e-01, -1.0050772e-17}, - {5.5557023e-01, 4.7094109e-17}, - {7.0710678e-01, -4.8336466e-17} -}; - -static const double fp64_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 double fp64_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} +template +struct constants_trig { + // 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} + }; }; // Notice the definition of the qd namespace @@ -665,8 +655,6 @@ inline two sin(const two &input) { local_2pi = fp32_2pi; local_pi2 = fp32_pi2; local_pi16 = fp32_pi16; - local_ptr_cos_table = &fp32_cos_table[0][0]; - local_ptr_sin_table = &fp32_sin_table[0][0]; zeropointzero = 0.0f; pointfive = 0.5f; local_nan = fp32_nan; @@ -674,8 +662,6 @@ inline two sin(const two &input) { local_2pi = fp64_2pi; local_pi2 = fp64_pi2; local_pi16 = fp64_pi16; - local_ptr_cos_table = &fp64_cos_table[0][0]; - local_ptr_sin_table = &fp64_sin_table[0][0]; zeropointzero = 0.0; pointfive = 0.5; local_nan = fp64_nan; @@ -683,6 +669,9 @@ inline two sin(const two &input) { static_assert(sizeof(T) == 0, "Other types not supported"); } + local_ptr_cos_table = &constants_trig::cos_table[0][0]; + local_ptr_sin_table = &constants_trig::sin_table[0][0]; + if (input.eval() == zeropointzero) { return two{zeropointzero, zeropointzero}; } From 29d835986fab25fb82493b529a4de97c9727d554 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 07:45:35 -0500 Subject: [PATCH 113/137] Moving inv_factor into templated struct. --- .../arithmetics/double-word-arithmetic.hpp | 63 +++++++------------ 1 file changed, 23 insertions(+), 40 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 46c5dcb..54733ce 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -260,44 +260,27 @@ const double fp64_split_factor_2 = 268435456.0; // 2^28 // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; -static const float fp32_inv_fact[n_inv_fact][2] = { - { 1.6666666e-01, 9.2518585e-18}, - { 4.1666666e-02, 2.3129646e-18}, - { 8.3333333e-03, 1.1564823e-19}, - { 1.3888888e-03, -5.3005439e-20}, - { 1.9841269e-04, 1.7209558e-22}, - { 2.4801587e-05, 2.1511947e-23}, - { 2.7557319e-06, -1.8583932e-22}, - { 2.7557319e-07, 2.3767714e-23}, - { 2.5052108e-08, -1.4488140e-24}, - { 2.0876756e-09, -1.2073450e-25}, - { 1.6059043e-10, 1.2585294e-26}, - { 1.1470745e-11, 2.0655512e-28}, - { 7.6471637e-13, 7.0387287e-30}, - { 4.7794773e-14, 4.3992054e-31}, - { 2.8114572e-15, 1.6508842e-31} -}; - -static const double fp64_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} -}; - template struct constants_trig { + + 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] = { @@ -474,17 +457,17 @@ static two sin_taylor(const two &input) { if constexpr (std::is_same_v) { pointfive = 0.5f; local_eps = fp32_eps; - local_ptr_inv_fact = &fp32_inv_fact[0][0]; zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { pointfive = 0.5; local_eps = fp64_eps; - local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } + local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; + //const T thresh = pointfive * std::abs(to_double(input)) * local_eps; const T thresh = pointfive * std::abs(input.eval()) * local_eps; @@ -529,18 +512,18 @@ static two cos_taylor(const two &input) { pointfive = 0.5f; local_eps = fp32_eps; onepointzero = 1.0f; - local_ptr_inv_fact = &fp32_inv_fact[0][0]; zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { pointfive = 0.5; local_eps = fp64_eps; onepointzero = 1.0; - local_ptr_inv_fact = &fp64_inv_fact[0][0]; zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } + local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; + const T thresh = pointfive * local_eps; two r, s, t, x; From 652ccf288470184543bf43fd2d63202c0aff92ac Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 08:57:11 -0500 Subject: [PATCH 114/137] Templating 2pi constant. --- .../arithmetics/double-word-arithmetic.hpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 54733ce..bab3143 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -237,7 +237,6 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp / inline.h / dd_real.cpp // TODO: transform this to FP32 (the original code is FP64) // TODO: check ranges (i.e., exponents) -const two fp32_2pi = two(6.2831853e+00, 2.4492935e-16); const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); static const float fp32_nan = std::numeric_limits::quiet_NaN(); @@ -247,7 +246,7 @@ const float fp32_qd_split_thresh = 6.6969287e+299; // = 2^996 // TODO: check! ra const float fp32_split_factor = 3.7252902e-09; // 2^-28 const float fp32_split_factor_2 = 268435456.0f; // 2^28 -const two fp64_2pi = two(6.283185307179586232e+00, 2.449293598294706414e-16); + const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); static const double fp64_nan = std::numeric_limits::quiet_NaN(); @@ -257,12 +256,17 @@ const double fp64_qd_split_thresh = 6.69692879491417e+299; // = 2^996 // TODO: c const double fp64_split_factor = 3.7252902984619140625e-09; // 2^-28 const double fp64_split_factor_2 = 268435456.0; // 2^28 +template +static const two _2pi = {6.283185307179586232e+00, 2.449293598294706414e-16}; + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; template struct constants_trig { + + static const constexpr T inv_fact[n_inv_fact][2] = { { 1.66666666666666657e-01, 9.25185853854297066e-18}, { 4.16666666666666644e-02, 2.31296463463574266e-18}, @@ -627,7 +631,6 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { template inline two sin(const two &input) { - two local_2pi; two local_pi2; two local_pi16; const T* local_ptr_cos_table; @@ -635,14 +638,12 @@ inline two sin(const two &input) { T zeropointzero, pointfive, local_nan; if constexpr (std::is_same_v) { - local_2pi = fp32_2pi; local_pi2 = fp32_pi2; local_pi16 = fp32_pi16; zeropointzero = 0.0f; pointfive = 0.5f; local_nan = fp32_nan; } else if constexpr (std::is_same_v) { - local_2pi = fp64_2pi; local_pi2 = fp64_pi2; local_pi16 = fp64_pi16; zeropointzero = 0.0; @@ -652,6 +653,8 @@ inline two sin(const two &input) { static_assert(sizeof(T) == 0, "Other types not supported"); } + two local_2pi = _2pi; + local_ptr_cos_table = &constants_trig::cos_table[0][0]; local_ptr_sin_table = &constants_trig::sin_table[0][0]; From 80d36af5361d5067a5991a7a90d92579f9787365 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 09:02:46 -0500 Subject: [PATCH 115/137] Templating pi2. --- .../arithmetics/double-word-arithmetic.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index bab3143..519cae8 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -237,7 +237,6 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp / inline.h / dd_real.cpp // TODO: transform this to FP32 (the original code is FP64) // TODO: check ranges (i.e., exponents) -const two fp32_pi2 = two(1.5707963e+00, 6.1232339e-17); const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); static const float fp32_nan = std::numeric_limits::quiet_NaN(); const float fp32_eps = 4.9303806e-32; // 2^-104 @@ -247,7 +246,7 @@ const float fp32_split_factor = 3.7252902e-09; // 2^-28 const float fp32_split_factor_2 = 268435456.0f; // 2^28 -const two fp64_pi2 = two(1.570796326794896558e+00, 6.123233995736766036e-17); + const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); static const double fp64_nan = std::numeric_limits::quiet_NaN(); const double fp64_eps = 4.93038065763132e-32; // 2^-104 @@ -259,6 +258,9 @@ const double fp64_split_factor_2 = 268435456.0; // 2^28 template static const two _2pi = {6.283185307179586232e+00, 2.449293598294706414e-16}; +template +static const two _pi2 = {1.570796326794896558e+00, 6.123233995736766036e-17}; + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -631,20 +633,18 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { template inline two sin(const two &input) { - two local_pi2; + two local_pi16; const T* local_ptr_cos_table; const T* local_ptr_sin_table; T zeropointzero, pointfive, local_nan; if constexpr (std::is_same_v) { - local_pi2 = fp32_pi2; local_pi16 = fp32_pi16; zeropointzero = 0.0f; pointfive = 0.5f; local_nan = fp32_nan; } else if constexpr (std::is_same_v) { - local_pi2 = fp64_pi2; local_pi16 = fp64_pi16; zeropointzero = 0.0; pointfive = 0.5; @@ -654,6 +654,7 @@ inline two sin(const two &input) { } two local_2pi = _2pi; + two local_pi2 = _pi2; local_ptr_cos_table = &constants_trig::cos_table[0][0]; local_ptr_sin_table = &constants_trig::sin_table[0][0]; From dcfc61fa2fb0a973f7707bb0cef3390c9fe29860 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 09:05:30 -0500 Subject: [PATCH 116/137] Templating pi16. --- .../arithmetics/double-word-arithmetic.hpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 519cae8..ae4c97c 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -237,7 +237,6 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp / inline.h / dd_real.cpp // TODO: transform this to FP32 (the original code is FP64) // TODO: check ranges (i.e., exponents) -const two fp32_pi16 = two(1.9634954e-01, 7.6540424e-18); static const float fp32_nan = std::numeric_limits::quiet_NaN(); const float fp32_eps = 4.9303806e-32; // 2^-104 const float fp32_qd_splitter = 134217729.0f; // = 2^27 + 1 @@ -247,7 +246,7 @@ const float fp32_split_factor_2 = 268435456.0f; // 2^28 -const two fp64_pi16 = two(1.963495408493620697e-01, 7.654042494670957545e-18); + static const double fp64_nan = std::numeric_limits::quiet_NaN(); const double fp64_eps = 4.93038065763132e-32; // 2^-104 const double fp64_qd_splitter = 134217729.0; // = 2^27 + 1 @@ -261,6 +260,9 @@ 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}; + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -632,20 +634,15 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { // Reference: QD / dd_real.cpp template inline two sin(const two &input) { - - - two local_pi16; const T* local_ptr_cos_table; const T* local_ptr_sin_table; T zeropointzero, pointfive, local_nan; if constexpr (std::is_same_v) { - local_pi16 = fp32_pi16; zeropointzero = 0.0f; pointfive = 0.5f; local_nan = fp32_nan; } else if constexpr (std::is_same_v) { - local_pi16 = fp64_pi16; zeropointzero = 0.0; pointfive = 0.5; local_nan = fp64_nan; @@ -655,6 +652,7 @@ inline two sin(const two &input) { two local_2pi = _2pi; two local_pi2 = _pi2; + two local_pi16 = _pi16; local_ptr_cos_table = &constants_trig::cos_table[0][0]; local_ptr_sin_table = &constants_trig::sin_table[0][0]; From ed16c7967a8d6593fd338241783c736fc4ad8866 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 09:10:51 -0500 Subject: [PATCH 117/137] Templating nan. --- .../arithmetics/double-word-arithmetic.hpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index ae4c97c..ef43ee2 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -237,7 +237,6 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp / inline.h / dd_real.cpp // TODO: transform this to FP32 (the original code is FP64) // TODO: check ranges (i.e., exponents) -static const float fp32_nan = std::numeric_limits::quiet_NaN(); const float fp32_eps = 4.9303806e-32; // 2^-104 const float fp32_qd_splitter = 134217729.0f; // = 2^27 + 1 const float fp32_qd_split_thresh = 6.6969287e+299; // = 2^996 // TODO: check! range is beyond float @@ -247,7 +246,7 @@ const float fp32_split_factor_2 = 268435456.0f; // 2^28 -static const double fp64_nan = std::numeric_limits::quiet_NaN(); + const double fp64_eps = 4.93038065763132e-32; // 2^-104 const double fp64_qd_splitter = 134217729.0; // = 2^27 + 1 const double fp64_qd_split_thresh = 6.69692879491417e+299; // = 2^996 // TODO: check! range is beyond double @@ -263,6 +262,9 @@ 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(); + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -571,21 +573,21 @@ two sqrt(const two &input) { only half the precision. */ - T zeropointzero, onepointzero, pointfive, local_nan; + T zeropointzero, onepointzero, pointfive; if constexpr (std::is_same_v) { onepointzero = 1.0f; zeropointzero = 0.0f; pointfive = 0.5f; - local_nan = fp32_nan; } else if constexpr (std::is_same_v) { onepointzero = 1.0; zeropointzero = 0.0; pointfive = 0.5; - local_nan = fp64_nan; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } + T local_nan = _nan; + if (input.eval() == zeropointzero) { return two{zeropointzero, zeropointzero}; } @@ -636,16 +638,14 @@ template inline two sin(const two &input) { const T* local_ptr_cos_table; const T* local_ptr_sin_table; - T zeropointzero, pointfive, local_nan; + T zeropointzero, pointfive; if constexpr (std::is_same_v) { zeropointzero = 0.0f; pointfive = 0.5f; - local_nan = fp32_nan; } else if constexpr (std::is_same_v) { zeropointzero = 0.0; pointfive = 0.5; - local_nan = fp64_nan; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } @@ -653,6 +653,7 @@ inline two sin(const two &input) { two local_2pi = _2pi; two local_pi2 = _pi2; two local_pi16 = _pi16; + T local_nan = _nan; local_ptr_cos_table = &constants_trig::cos_table[0][0]; local_ptr_sin_table = &constants_trig::sin_table[0][0]; From 5d3de1c173b75de858bde44fb64d9cbcfc1709b1 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 09:14:58 -0500 Subject: [PATCH 118/137] Templating eps. --- .../arithmetics/double-word-arithmetic.hpp | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index ef43ee2..798039b 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -237,17 +237,11 @@ inline two div(const two &x, const two &y) { // Reference: QD / dd_const.cpp / inline.h / dd_real.cpp // TODO: transform this to FP32 (the original code is FP64) // TODO: check ranges (i.e., exponents) -const float fp32_eps = 4.9303806e-32; // 2^-104 const float fp32_qd_splitter = 134217729.0f; // = 2^27 + 1 const float fp32_qd_split_thresh = 6.6969287e+299; // = 2^996 // TODO: check! range is beyond float const float fp32_split_factor = 3.7252902e-09; // 2^-28 const float fp32_split_factor_2 = 268435456.0f; // 2^28 - - - - -const double fp64_eps = 4.93038065763132e-32; // 2^-104 const double fp64_qd_splitter = 134217729.0; // = 2^27 + 1 const double fp64_qd_split_thresh = 6.69692879491417e+299; // = 2^996 // TODO: check! range is beyond double const double fp64_split_factor = 3.7252902984619140625e-09; // 2^-28 @@ -265,6 +259,9 @@ 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 + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -461,16 +458,14 @@ inline two sqr(const two &input) { template static two sin_taylor(const two &input) { - T pointfive, local_eps; + T pointfive; const T* local_ptr_inv_fact; T zeropointzero; if constexpr (std::is_same_v) { pointfive = 0.5f; - local_eps = fp32_eps; zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { pointfive = 0.5; - local_eps = fp64_eps; zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); @@ -478,6 +473,8 @@ static two sin_taylor(const two &input) { local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; + T local_eps = _eps; + //const T thresh = pointfive * std::abs(to_double(input)) * local_eps; const T thresh = pointfive * std::abs(input.eval()) * local_eps; @@ -515,17 +512,15 @@ inline two mul_pwr2(const two &input, T b) { template static two cos_taylor(const two &input) { - T pointfive, local_eps, onepointzero; + T pointfive, onepointzero; const T* local_ptr_inv_fact; T zeropointzero; if constexpr (std::is_same_v) { pointfive = 0.5f; - local_eps = fp32_eps; onepointzero = 1.0f; zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { pointfive = 0.5; - local_eps = fp64_eps; onepointzero = 1.0; zeropointzero = 0.0; } else { @@ -534,6 +529,8 @@ static two cos_taylor(const two &input) { local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; + T local_eps = _eps; + const T thresh = pointfive * local_eps; two r, s, t, x; From 6efcb934eb14f92a94f5a77540431d27522b37c8 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 09:42:53 -0500 Subject: [PATCH 119/137] Removing unnecessary constants. --- .../arithmetics/double-word-arithmetic.hpp | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 798039b..da2699d 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -235,18 +235,6 @@ inline two div(const two &x, const two &y) { } // Reference: QD / dd_const.cpp / inline.h / dd_real.cpp -// TODO: transform this to FP32 (the original code is FP64) -// TODO: check ranges (i.e., exponents) -const float fp32_qd_splitter = 134217729.0f; // = 2^27 + 1 -const float fp32_qd_split_thresh = 6.6969287e+299; // = 2^996 // TODO: check! range is beyond float -const float fp32_split_factor = 3.7252902e-09; // 2^-28 -const float fp32_split_factor_2 = 268435456.0f; // 2^28 - -const double fp64_qd_splitter = 134217729.0; // = 2^27 + 1 -const double fp64_qd_split_thresh = 6.69692879491417e+299; // = 2^996 // TODO: check! range is beyond double -const double fp64_split_factor = 3.7252902984619140625e-09; // 2^-28 -const double fp64_split_factor_2 = 268435456.0; // 2^28 - template static const two _2pi = {6.283185307179586232e+00, 2.449293598294706414e-16}; @@ -268,8 +256,6 @@ static const int n_inv_fact = 15; template struct constants_trig { - - static const constexpr T inv_fact[n_inv_fact][2] = { { 1.66666666666666657e-01, 9.25185853854297066e-18}, { 4.16666666666666644e-02, 2.31296463463574266e-18}, From 3b3ab436f08d1d3d1e95eb8499a1c5ef7e1cd613 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 09:48:23 -0500 Subject: [PATCH 120/137] Compacting code. --- .../arithmetics/double-word-arithmetic.hpp | 22 ++++++------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index da2699d..cff0bcc 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -444,9 +444,7 @@ inline two sqr(const two &input) { template static two sin_taylor(const two &input) { - T pointfive; - const T* local_ptr_inv_fact; - T zeropointzero; + T pointfive, zeropointzero; if constexpr (std::is_same_v) { pointfive = 0.5f; zeropointzero = 0.0f; @@ -457,8 +455,7 @@ static two sin_taylor(const two &input) { static_assert(sizeof(T) == 0, "Other types not supported"); } - local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; - + const T* local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; T local_eps = _eps; //const T thresh = pointfive * std::abs(to_double(input)) * local_eps; @@ -498,9 +495,7 @@ inline two mul_pwr2(const two &input, T b) { template static two cos_taylor(const two &input) { - T pointfive, onepointzero; - const T* local_ptr_inv_fact; - T zeropointzero; + T pointfive, onepointzero, zeropointzero; if constexpr (std::is_same_v) { pointfive = 0.5f; onepointzero = 1.0f; @@ -513,8 +508,7 @@ static two cos_taylor(const two &input) { static_assert(sizeof(T) == 0, "Other types not supported"); } - local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; - + const T* local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; T local_eps = _eps; const T thresh = pointfive * local_eps; @@ -619,10 +613,8 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { // Reference: QD / dd_real.cpp template inline two sin(const two &input) { - const T* local_ptr_cos_table; - const T* local_ptr_sin_table; - T zeropointzero, pointfive; + T zeropointzero, pointfive; if constexpr (std::is_same_v) { zeropointzero = 0.0f; pointfive = 0.5f; @@ -638,8 +630,8 @@ inline two sin(const two &input) { two local_pi16 = _pi16; T local_nan = _nan; - local_ptr_cos_table = &constants_trig::cos_table[0][0]; - local_ptr_sin_table = &constants_trig::sin_table[0][0]; + const T* local_ptr_cos_table = &constants_trig::cos_table[0][0]; + const T* local_ptr_sin_table = &constants_trig::sin_table[0][0]; if (input.eval() == zeropointzero) { return two{zeropointzero, zeropointzero}; From e48dc0d611786e4c3ef04c4cdc156b290fde8d5c Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 09:50:51 -0500 Subject: [PATCH 121/137] Renaming struct. --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index cff0bcc..1698c1f 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -254,7 +254,7 @@ static const constexpr T _eps = 4.93038065763132e-32; // 2^-104 static const int n_inv_fact = 15; template -struct constants_trig { +struct constants_trig_tables { static const constexpr T inv_fact[n_inv_fact][2] = { { 1.66666666666666657e-01, 9.25185853854297066e-18}, @@ -455,7 +455,7 @@ static two sin_taylor(const two &input) { static_assert(sizeof(T) == 0, "Other types not supported"); } - const T* local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; + const T* local_ptr_inv_fact = &constants_trig_tables::inv_fact[0][0]; T local_eps = _eps; //const T thresh = pointfive * std::abs(to_double(input)) * local_eps; @@ -508,7 +508,7 @@ static two cos_taylor(const two &input) { static_assert(sizeof(T) == 0, "Other types not supported"); } - const T* local_ptr_inv_fact = &constants_trig::inv_fact[0][0]; + const T* local_ptr_inv_fact = &constants_trig_tables::inv_fact[0][0]; T local_eps = _eps; const T thresh = pointfive * local_eps; @@ -630,8 +630,8 @@ inline two sin(const two &input) { two local_pi16 = _pi16; T local_nan = _nan; - const T* local_ptr_cos_table = &constants_trig::cos_table[0][0]; - const T* local_ptr_sin_table = &constants_trig::sin_table[0][0]; + 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}; From 3ffd66cb8a9664e77d20c19f999d134168c57968 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 10:31:03 -0500 Subject: [PATCH 122/137] Defining and using templated variable pointfive. --- .../arithmetics/double-word-arithmetic.hpp | 50 +++++++------------ 1 file changed, 18 insertions(+), 32 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 1698c1f..127ea3a 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -250,6 +250,10 @@ static const constexpr T _nan = std::numeric_limits::quiet_NaN(); template static const constexpr T _eps = 4.93038065763132e-32; // 2^-104 +// Some recurrently-used contants +template +static const constexpr T pointfive = 0.5; + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -297,20 +301,12 @@ namespace qd { /* Computes the nearest integer to input. */ template inline T nint(T input) { - T pointfive; - - if constexpr (std::is_same_v) { - pointfive = 0.5f; - } else if constexpr (std::is_same_v) { - pointfive = 0.5; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + 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); + return std::floor(input + pointfive); } // Reference: QD / inline.h @@ -377,14 +373,12 @@ inline two nint(const two &input) { T hi = qd::nint(input.h); T lo; - T zeropointzero, pointfive, onepointzero; + T zeropointzero, onepointzero; if constexpr (std::is_same_v) { zeropointzero = 0.0f; - pointfive = 0.5f; onepointzero = 1.0f; } else if constexpr (std::is_same_v) { zeropointzero = 0.0; - pointfive = 0.5; onepointzero = 1.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); @@ -403,7 +397,7 @@ inline two nint(const two &input) { else { /* High word is not an integer */ lo = zeropointzero; - if(std::abs(hi - input.h) == pointfive && input.l < 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. */ @@ -444,12 +438,10 @@ inline two sqr(const two &input) { template static two sin_taylor(const two &input) { - T pointfive, zeropointzero; + T zeropointzero; if constexpr (std::is_same_v) { - pointfive = 0.5f; zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { - pointfive = 0.5; zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); @@ -459,7 +451,7 @@ static two sin_taylor(const two &input) { T local_eps = _eps; //const T thresh = pointfive * std::abs(to_double(input)) * local_eps; - const T thresh = pointfive * std::abs(input.eval()) * local_eps; + const T thresh = pointfive * std::abs(input.eval()) * local_eps; two r, s, t, x; @@ -495,13 +487,11 @@ inline two mul_pwr2(const two &input, T b) { template static two cos_taylor(const two &input) { - T pointfive, onepointzero, zeropointzero; + T onepointzero, zeropointzero; if constexpr (std::is_same_v) { - pointfive = 0.5f; onepointzero = 1.0f; zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { - pointfive = 0.5; onepointzero = 1.0; zeropointzero = 0.0; } else { @@ -511,7 +501,7 @@ static two cos_taylor(const two &input) { const T* local_ptr_inv_fact = &constants_trig_tables::inv_fact[0][0]; T local_eps = _eps; - const T thresh = pointfive * local_eps; + const T thresh = pointfive * local_eps; two r, s, t, x; @@ -522,7 +512,7 @@ static two cos_taylor(const two &input) { two temp = sqr(input); x = two{-temp.h, -temp.l}; r = x; - s = add(onepointzero, mul_pwr2(r, pointfive)); + s = add(onepointzero, mul_pwr2(r, pointfive)); int i = 1; do { r = mul(r, x); @@ -550,15 +540,13 @@ two sqrt(const two &input) { only half the precision. */ - T zeropointzero, onepointzero, pointfive; + T zeropointzero, onepointzero; if constexpr (std::is_same_v) { onepointzero = 1.0f; zeropointzero = 0.0f; - pointfive = 0.5f; } else if constexpr (std::is_same_v) { onepointzero = 1.0; zeropointzero = 0.0; - pointfive = 0.5; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } @@ -578,7 +566,7 @@ two sqrt(const two &input) { T x = onepointzero / std::sqrt(input.h); T ax = input.h * x; two temp = sub(input, dd_real::sqr(ax)); - return dd_real::add(ax, temp.h * x * pointfive); + return dd_real::add(ax, temp.h * x * pointfive); } // Reference: QD / dd_real.cpp @@ -614,13 +602,11 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { template inline two sin(const two &input) { - T zeropointzero, pointfive; + T zeropointzero; if constexpr (std::is_same_v) { zeropointzero = 0.0f; - pointfive = 0.5f; } else if constexpr (std::is_same_v) { zeropointzero = 0.0; - pointfive = 0.5; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } @@ -652,10 +638,10 @@ inline two sin(const two &input) { // Approximately reducing modulo pi/2 and then modulo pi/16 //TODO: original type is double, here it is templated - T q = std::floor(r.h / local_pi2.h + pointfive); + 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); + 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); From 800b777d24967acc2369cbff64d38228cbf07af9 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 10:34:30 -0500 Subject: [PATCH 123/137] Defining and using templated variable twopointzero. --- .../arithmetics/double-word-arithmetic.hpp | 27 +++++-------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 127ea3a..fafe211 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -254,6 +254,9 @@ static const constexpr T _eps = 4.93038065763132e-32; // 2^-104 template static const constexpr T pointfive = 0.5; +template +static const constexpr T twopointzero = 2.0; + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -313,15 +316,7 @@ namespace qd { /* Computes fl(input * input) and err(input * input) */ template inline T two_sqr(T input, T &err) { - - T twopointzero; - if constexpr (std::is_same_v) { - twopointzero = 2.0f; - } else if constexpr (std::is_same_v) { - twopointzero = 2.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); T hi, lo; T q = input * input; @@ -329,7 +324,7 @@ namespace qd { two temp = twofloat::algorithms::Split(input); hi = temp.h; lo = temp.l; - err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; + err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; return q; } @@ -410,20 +405,12 @@ inline two nint(const two &input) { // Reference: QD / dd_inline.h template inline two sqr(const two &input) { - - T twopointzero; - if constexpr (std::is_same_v) { - twopointzero = 2.0f; - } else if constexpr (std::is_same_v) { - twopointzero = 2.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); T p1, p2; T s1, s2; p1 = qd::two_sqr(input.h, p2); - p2 += twopointzero * input.h * input.l; + p2 += twopointzero * input.h * input.l; p2 += input.l * input.l; //s1 = qd::quick_two_sum(p1, p2, s2); two temp = twofloat::algorithms::FastTwoSum(p1, p2); From aaf8057b6ac01ab887a4a4a68afb0050377d3275 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 10:40:19 -0500 Subject: [PATCH 124/137] Defining and using templated variable zeropointzero. --- .../arithmetics/double-word-arithmetic.hpp | 69 +++++++------------ 1 file changed, 24 insertions(+), 45 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index fafe211..dbb6ac1 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -257,6 +257,9 @@ static const constexpr T pointfive = 0.5; template static const constexpr T twopointzero = 2.0; +template +static const constexpr T zeropointzero = 0.0; + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -368,12 +371,10 @@ inline two nint(const two &input) { T hi = qd::nint(input.h); T lo; - T zeropointzero, onepointzero; + T onepointzero; if constexpr (std::is_same_v) { - zeropointzero = 0.0f; onepointzero = 1.0f; } else if constexpr (std::is_same_v) { - zeropointzero = 0.0; onepointzero = 1.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); @@ -391,8 +392,8 @@ inline two nint(const two &input) { } else { /* High word is not an integer */ - lo = zeropointzero; - if(std::abs(hi - input.h) == pointfive && input.l < zeropointzero) { + 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. */ @@ -424,15 +425,7 @@ inline two sqr(const two &input) { Assumes |a| <= pi/32 */ template static two sin_taylor(const two &input) { - - T zeropointzero; - if constexpr (std::is_same_v) { - zeropointzero = 0.0f; - } else if constexpr (std::is_same_v) { - zeropointzero = 0.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + 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; @@ -442,8 +435,8 @@ static two sin_taylor(const two &input) { two r, s, t, x; - if (input.eval() == zeropointzero) { - return two{zeropointzero, zeropointzero}; + if (input.eval() == zeropointzero) { + return two{zeropointzero, zeropointzero}; } int i = 0; @@ -474,13 +467,11 @@ inline two mul_pwr2(const two &input, T b) { template static two cos_taylor(const two &input) { - T onepointzero, zeropointzero; + T onepointzero; if constexpr (std::is_same_v) { onepointzero = 1.0f; - zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { onepointzero = 1.0; - zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } @@ -492,8 +483,8 @@ static two cos_taylor(const two &input) { two r, s, t, x; - if (input.eval() == zeropointzero) { - return two{onepointzero, zeropointzero}; + if (input.eval() == zeropointzero) { + return two{onepointzero, zeropointzero}; } two temp = sqr(input); @@ -527,26 +518,24 @@ two sqrt(const two &input) { only half the precision. */ - T zeropointzero, onepointzero; + T onepointzero; if constexpr (std::is_same_v) { onepointzero = 1.0f; - zeropointzero = 0.0f; } else if constexpr (std::is_same_v) { onepointzero = 1.0; - zeropointzero = 0.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); } T local_nan = _nan; - if (input.eval() == zeropointzero) { - return two{zeropointzero, zeropointzero}; + if (input.eval() == zeropointzero) { + return two{zeropointzero, zeropointzero}; } // ERROR: Negative argument // if (input.is_negative()) - if (input.h < zeropointzero) { + if (input.h < zeropointzero) { return two{local_nan, local_nan}; } @@ -560,12 +549,10 @@ two sqrt(const two &input) { template static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { - T zeropointzero, onepointzero; + T onepointzero; if constexpr (std::is_same_v) { - zeropointzero = 0.0f; onepointzero = 1.0f; } else if constexpr (std::is_same_v) { - zeropointzero = 0.0; onepointzero = 1.0; } else { static_assert(sizeof(T) == 0, "Other types not supported"); @@ -573,11 +560,11 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { // Reference: QD / dd_inline.h // Assignment is performed according to operator== - if (input.eval() == zeropointzero) { - sin_a.h = zeropointzero; - sin_a.l = zeropointzero; + if (input.eval() == zeropointzero) { + sin_a.h = zeropointzero; + sin_a.l = zeropointzero; cos_a.h = onepointzero; - cos_a.l = zeropointzero; + cos_a.l = zeropointzero; return; } @@ -588,15 +575,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { // Reference: QD / dd_real.cpp template inline two sin(const two &input) { - - T zeropointzero; - if constexpr (std::is_same_v) { - zeropointzero = 0.0f; - } else if constexpr (std::is_same_v) { - zeropointzero = 0.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); two local_2pi = _2pi; two local_pi2 = _pi2; @@ -606,8 +585,8 @@ inline two sin(const two &input) { 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}; + if (input.eval() == zeropointzero) { + return two{zeropointzero, zeropointzero}; } // Approximately reducing modulo 2*pi From 9e1e0be7403cab689209d28400aef46366224a6d Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 10:44:58 -0500 Subject: [PATCH 125/137] Defining and using templated variable onepointzero. --- .../arithmetics/double-word-arithmetic.hpp | 54 +++++-------------- 1 file changed, 14 insertions(+), 40 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index dbb6ac1..35bf177 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -260,6 +260,9 @@ static const constexpr T twopointzero = 2.0; template static const constexpr T zeropointzero = 0.0; +template +static const constexpr T onepointzero = 1.0; + // Reference: QD / dd_real.cpp static const int n_inv_fact = 15; @@ -368,18 +371,11 @@ namespace dd_real { /* Round to nearest integer */ template inline two nint(const two &input) { + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); + T hi = qd::nint(input.h); T lo; - T onepointzero; - if constexpr (std::is_same_v) { - onepointzero = 1.0f; - } else if constexpr (std::is_same_v) { - onepointzero = 1.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } - if(hi == input.h) { /* High word is an integer already. Round the low word. */ lo = qd::nint(input.l); @@ -396,7 +392,7 @@ inline two nint(const two &input) { 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. */ + hi -= onepointzero; /* NOTE: This does not cause INEXACT. */ } } @@ -467,14 +463,7 @@ inline two mul_pwr2(const two &input, T b) { template static two cos_taylor(const two &input) { - T onepointzero; - if constexpr (std::is_same_v) { - onepointzero = 1.0f; - } else if constexpr (std::is_same_v) { - onepointzero = 1.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + 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; @@ -484,13 +473,13 @@ static two cos_taylor(const two &input) { two r, s, t, x; if (input.eval() == zeropointzero) { - return two{onepointzero, zeropointzero}; + return two{onepointzero, zeropointzero}; } two temp = sqr(input); x = two{-temp.h, -temp.l}; r = x; - s = add(onepointzero, mul_pwr2(r, pointfive)); + s = add(onepointzero, mul_pwr2(r, pointfive)); int i = 1; do { r = mul(r, x); @@ -518,14 +507,7 @@ two sqrt(const two &input) { only half the precision. */ - T onepointzero; - if constexpr (std::is_same_v) { - onepointzero = 1.0f; - } else if constexpr (std::is_same_v) { - onepointzero = 1.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); T local_nan = _nan; @@ -539,7 +521,7 @@ two sqrt(const two &input) { return two{local_nan, local_nan}; } - T x = onepointzero / std::sqrt(input.h); + T x = onepointzero / std::sqrt(input.h); T ax = input.h * x; two temp = sub(input, dd_real::sqr(ax)); return dd_real::add(ax, temp.h * x * pointfive); @@ -548,28 +530,20 @@ two sqrt(const two &input) { // Reference: QD / dd_real.cpp template static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { - - T onepointzero; - if constexpr (std::is_same_v) { - onepointzero = 1.0f; - } else if constexpr (std::is_same_v) { - onepointzero = 1.0; - } else { - static_assert(sizeof(T) == 0, "Other types not supported"); - } + static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); // Reference: QD / dd_inline.h // Assignment is performed according to operator== if (input.eval() == zeropointzero) { sin_a.h = zeropointzero; sin_a.l = zeropointzero; - cos_a.h = onepointzero; + cos_a.h = onepointzero; cos_a.l = zeropointzero; return; } sin_a = sin_taylor(input); - cos_a = sqrt(sub(onepointzero, sqr(sin_a))); + cos_a = sqrt(sub(onepointzero, sqr(sin_a))); } // Reference: QD / dd_real.cpp From 8b940f46961fafccc0d321a800de1c441e6ffa1d Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 14:21:17 -0500 Subject: [PATCH 126/137] Moving definition of nint() outside qd namespace. Note that nint() is now overloaded. --- .../arithmetics/double-word-arithmetic.hpp | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 35bf177..942a7de 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -306,18 +306,6 @@ struct constants_trig_tables { // Notice the definition of the qd namespace namespace qd { - // Reference: QD / inline.h - /* Computes the nearest integer to input. */ - 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); - } - // Reference: QD / inline.h /* Computes fl(input * input) and err(input * input) */ template @@ -367,18 +355,29 @@ namespace dd_real { } // End namespace dd_real +// Reference: QD / inline.h +/* Computes the nearest integer to input. */ +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); +} + // Reference: QD / dd_inline.h /* Round to nearest integer */ template inline two nint(const two &input) { static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); - T hi = qd::nint(input.h); + T hi = /*qd::*/nint(input.h); T lo; if(hi == input.h) { /* High word is an integer already. Round the low word. */ - lo = qd::nint(input.l); + lo = /*qd::*/nint(input.l); /* Renormalize. This is needed if h = some integer, l = 1/2. */ //hi = qd::quick_two_sum(hi, lo, lo); From 2d7c5e4f28907cec3b2f15ae946f07844e8d1767 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 14:26:40 -0500 Subject: [PATCH 127/137] Moving definition of two_sqr() outside qd namespace. --- .../arithmetics/double-word-arithmetic.hpp | 37 +++++++++---------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 942a7de..dc29525 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -305,23 +305,6 @@ struct constants_trig_tables { // Notice the definition of the qd namespace namespace qd { - - // Reference: QD / inline.h - /* Computes fl(input * input) and err(input * input) */ - 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; - //qd::split(input, hi, lo); - 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 // TODO: make sure this has not been already implemented template @@ -334,13 +317,29 @@ namespace qd { } // End namespace qd +// Reference: QD / inline.h +/* Computes fl(input * input) and err(input * input) */ +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; + //qd::split(input, hi, lo); + two temp = twofloat::algorithms::Split(input); + hi = temp.h; + lo = temp.l; + err = ((hi * hi - q) + twopointzero * hi * lo) + lo * lo; + return q; +} + namespace dd_real { // Reference: QD / dd_inline.h template inline two sqr(T input) { T p1, p2; - p1 = qd::two_sqr(input, p2); + p1 = two_sqr(input, p2); return two{p1, p2}; } @@ -405,7 +404,7 @@ inline two sqr(const two &input) { T p1, p2; T s1, s2; - p1 = qd::two_sqr(input.h, p2); + p1 = two_sqr(input.h, p2); p2 += twopointzero * input.h * input.l; p2 += input.l * input.l; //s1 = qd::quick_two_sum(p1, p2, s2); From 84ef26661867f3c55faf1737726c2aa50d7d7a13 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 14:29:29 -0500 Subject: [PATCH 128/137] Moving definition of two_sum() outside qd namespace. --- .../arithmetics/double-word-arithmetic.hpp | 26 ++++++++----------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index dc29525..19fe64a 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -303,20 +303,6 @@ struct constants_trig_tables { }; }; -// Notice the definition of the qd namespace -namespace qd { - // Reference QD / inline.h - // TODO: make sure this has not been already implemented - 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; - } - -} // End namespace qd - // Reference: QD / inline.h /* Computes fl(input * input) and err(input * input) */ template @@ -333,6 +319,16 @@ inline T two_sqr(T input, T &err) { return q; } +// Reference QD / inline.h +// TODO: make sure this has not been already implemented +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; +} + namespace dd_real { // Reference: QD / dd_inline.h @@ -348,7 +344,7 @@ namespace dd_real { template inline two add(T a, T b) { T s, e; - s = qd::two_sum(a, b, e); + s = two_sum(a, b, e); return two (s, e); } From ab0c884667d064464e15cb40acccb8bbaa25c5c7 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 14:31:31 -0500 Subject: [PATCH 129/137] Additional comment. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 19fe64a..3b96cc9 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -303,7 +303,7 @@ struct constants_trig_tables { }; }; -// Reference: QD / inline.h +// Reference: QD / inline.h (within qd namespace) /* Computes fl(input * input) and err(input * input) */ template inline T two_sqr(T input, T &err) { @@ -319,7 +319,7 @@ inline T two_sqr(T input, T &err) { return q; } -// Reference QD / inline.h +// Reference QD / inline.h (within qd namespace) // TODO: make sure this has not been already implemented template inline T two_sum(T a, T b, T &err) { @@ -350,7 +350,7 @@ namespace dd_real { } // End namespace dd_real -// Reference: QD / inline.h +// Reference: QD / inline.h (within qd namespace) /* Computes the nearest integer to input. */ template inline T nint(T input) { From 50b5edae82ee51650f66c54a9ad7fada5504d00a Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 15:22:33 -0500 Subject: [PATCH 130/137] Moving definition of sqr() outside dd_real namespace. --- .../arithmetics/double-word-arithmetic.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 3b96cc9..2c5b6fe 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -331,14 +331,6 @@ inline T two_sum(T a, T b, T &err) { namespace dd_real { - // Reference: QD / dd_inline.h - template - inline two sqr(T input) { - T p1, p2; - p1 = two_sqr(input, p2); - return two{p1, p2}; - } - // Reference: QD / dd_inline.h // TODO: make sure this has not been already implemented template @@ -393,6 +385,14 @@ inline two nint(const two &input) { 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) { @@ -517,7 +517,7 @@ two sqrt(const two &input) { T x = onepointzero / std::sqrt(input.h); T ax = input.h * x; - two temp = sub(input, dd_real::sqr(ax)); + two temp = sub(input, sqr(ax)); return dd_real::add(ax, temp.h * x * pointfive); } From 242118023c7c5b21ca0d29ab4f3db7e1e0e5d517 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 15:30:47 -0500 Subject: [PATCH 131/137] Moving definition of add() outside dd_real namespace. Note that such add() is overloaded. --- .../arithmetics/double-word-arithmetic.hpp | 22 ++++++++----------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 2c5b6fe..52272cc 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -329,18 +329,14 @@ inline T two_sum(T a, T b, T &err) { return s; } -namespace dd_real { - - // Reference: QD / dd_inline.h - // TODO: make sure this has not been already implemented - template - inline two add(T a, T b) { - T s, e; - s = two_sum(a, b, e); - return two (s, e); - } - -} // End namespace dd_real +// Reference: QD / dd_inline.h (within dd_real namespace) +// TODO: make sure this has not been already implemented +template +inline two add(T a, T b) { + T s, e; + s = two_sum(a, b, e); + return two (s, e); +} // Reference: QD / inline.h (within qd namespace) /* Computes the nearest integer to input. */ @@ -518,7 +514,7 @@ two sqrt(const two &input) { T x = onepointzero / std::sqrt(input.h); T ax = input.h * x; two temp = sub(input, sqr(ax)); - return dd_real::add(ax, temp.h * x * pointfive); + return add(ax, temp.h * x * pointfive); } // Reference: QD / dd_real.cpp From 608bbdaf3e4fb97164de40f5c87f05180de158d9 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 15:43:57 -0500 Subject: [PATCH 132/137] Minor removals. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 52272cc..2783764 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -355,12 +355,12 @@ template inline two nint(const two &input) { static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); - T hi = /*qd::*/nint(input.h); + T hi = nint(input.h); T lo; if(hi == input.h) { /* High word is an integer already. Round the low word. */ - lo = /*qd::*/nint(input.l); + lo = nint(input.l); /* Renormalize. This is needed if h = some integer, l = 1/2. */ //hi = qd::quick_two_sum(hi, lo, lo); From 274696b93d5025e26d7fc05ebc1cbc0bdfc63f87 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 15:51:55 -0500 Subject: [PATCH 133/137] Cleaning up. --- .../arithmetics/double-word-arithmetic.hpp | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 2783764..f3b4c96 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -308,10 +308,8 @@ struct constants_trig_tables { 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; - //qd::split(input, hi, lo); two temp = twofloat::algorithms::Split(input); hi = temp.h; lo = temp.l; @@ -354,10 +352,8 @@ inline T nint(T input) { 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); @@ -393,13 +389,11 @@ inline two sqr(T input) { 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; - //s1 = qd::quick_two_sum(p1, p2, s2); two temp = twofloat::algorithms::FastTwoSum(p1, p2); s1 = temp.h; s2 = temp.l; @@ -415,10 +409,7 @@ static two sin_taylor(const two &input) { const T* local_ptr_inv_fact = &constants_trig_tables::inv_fact[0][0]; T local_eps = _eps; - - //const T thresh = pointfive * std::abs(to_double(input)) * local_eps; const T thresh = pointfive * std::abs(input.eval()) * local_eps; - two r, s, t, x; if (input.eval() == zeropointzero) { @@ -452,14 +443,11 @@ inline two mul_pwr2(const two &input, T 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) { @@ -476,7 +464,6 @@ static two cos_taylor(const two &input) { 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(to_double(t)) > thresh); } while (i < n_inv_fact && std::abs(t.eval()) > thresh); return s; @@ -498,7 +485,6 @@ 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) { @@ -523,7 +509,6 @@ 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"); // Reference: QD / dd_inline.h - // Assignment is performed according to operator== if (input.eval() == zeropointzero) { sin_a.h = zeropointzero; sin_a.l = zeropointzero; From 59e215d7a83855e671e88a5db174fec1c679e7fe Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 16:16:00 -0500 Subject: [PATCH 134/137] Parameterizing the mode and useFMA for sin(). --- .../arithmetics/double-word-arithmetic.hpp | 66 +++++++++---------- test/double-word-arithmetic.test.cpp | 4 +- 2 files changed, 35 insertions(+), 35 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index f3b4c96..7b9e3dd 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -403,7 +403,7 @@ inline two sqr(const two &input) { // Reference: QD / dd_real.cpp /* Computes sin(a) using Taylor series. Assumes |a| <= pi/32 */ -template +template static two sin_taylor(const two &input) { static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); @@ -423,9 +423,9 @@ static two sin_taylor(const two &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); + 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(to_double(t)) > thresh); } while (i < n_inv_fact && std::abs(t.eval()) > thresh); @@ -441,7 +441,7 @@ inline two mul_pwr2(const two &input, T b) { } // Reference: QD / dd_real.cpp -template +template static two cos_taylor(const two &input) { static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); @@ -460,9 +460,9 @@ static two cos_taylor(const two &input) { 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); + 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); @@ -472,7 +472,7 @@ static two cos_taylor(const two &input) { // Reference: QD / dd_real.cpp /* Computes the square root of the double-double number dd. NOTE: dd must be a non-negative number. */ -template +template two sqrt(const two &input) { /* Strategy: Use Karp's trick: if x is an approximation to sqrt(a), then @@ -499,12 +499,12 @@ two sqrt(const two &input) { T x = onepointzero / std::sqrt(input.h); T ax = input.h * x; - two temp = sub(input, sqr(ax)); + two temp = sub

(input, sqr(ax)); return add(ax, temp.h * x * pointfive); } // Reference: QD / dd_real.cpp -template +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"); @@ -517,12 +517,12 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { return; } - sin_a = sin_taylor(input); - cos_a = sqrt(sub(onepointzero, sqr(sin_a))); + sin_a = sin_taylor(input); + cos_a = sqrt

(sub(onepointzero, sqr(sin_a))); } // Reference: QD / dd_real.cpp -template +template inline two sin(const two &input) { static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); @@ -539,8 +539,8 @@ inline two sin(const two &input) { } // Approximately reducing modulo 2*pi - two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div (accurate + fma) - two r = sub(input, mul(local_2pi, z)); + two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div (accurate + fma) + two r = sub

(input, mul(local_2pi, z)); /* std::cout << "LSV starts ..." << std::endl; @@ -554,10 +554,10 @@ inline two sin(const two &input) { //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); - two t = sub(r, mul(local_pi2, q)); + 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)); + t = sub

(t, mul(local_pi16, q)); int k = static_cast(q); int abs_k = std::abs(k); @@ -584,13 +584,13 @@ inline two sin(const two &input) { if (k == 0) { switch(j) { case 0: - return sin_taylor(t); + return sin_taylor(t); case 1: - return cos_taylor(t); + return cos_taylor(t); case -1: - return two{-cos_taylor(t).h, -cos_taylor(t).l}; // Reference: QD / dd_inline.h + 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}; + return two{-sin_taylor(t).h, -sin_taylor(t).l}; } } @@ -601,7 +601,7 @@ inline two sin(const two &input) { two sin_t, cos_t; - sincos_taylor(t, sin_t, cos_t); + sincos_taylor(t, sin_t, cos_t); /* std::cout << "LSV starts ..." << std::endl; @@ -612,15 +612,15 @@ inline two sin(const two &input) { if (j == 0) { if (k > 0) { - r = add(mul(u, sin_t), mul(v, cos_t)); + r = add

(mul(u, sin_t), mul(v, cos_t)); } else { - r = sub(mul(u, sin_t), mul(v, cos_t)); + 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)); + r = sub

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

(mul(u, cos_t), mul(v, sin_t)); /* std::cout << "LSV starts ..." << std::endl; std::cout << "u = " << u.eval() << std::endl; @@ -633,19 +633,19 @@ inline two sin(const two &input) { } } else if (j == -1) { if (k > 0) { - r = sub(mul(v, sin_t), mul(u, cos_t)); + r = sub

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

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

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

(mul(v, cos_t), mul(u, sin_t)); } } /* diff --git a/test/double-word-arithmetic.test.cpp b/test/double-word-arithmetic.test.cpp index fb6597c..ed0b690 100644 --- a/test/double-word-arithmetic.test.cpp +++ b/test/double-word-arithmetic.test.cpp @@ -260,7 +260,7 @@ TEST(DoubleWordArithmetic, SinTest) { float sinx = std::sin(x); two two_x = {x, 0.0f}; - two sin_two_x = doubleword::sin(two_x); + 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()); @@ -272,7 +272,7 @@ TEST(DoubleWordArithmetic, SinTest2) { double sinx = std::sin(x); two two_x = {x, 0.0}; - two sin_two_x = doubleword::sin(two_x); + 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()); From 912281b3fe321362efa66dc6980b1ed80ae1d588 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 16:17:24 -0500 Subject: [PATCH 135/137] Minor cleanup. --- include/libtwofloat/arithmetics/double-word-arithmetic.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 7b9e3dd..0daa595 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -427,7 +427,6 @@ static two sin_taylor(const two &input) { 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(to_double(t)) > thresh); } while (i < n_inv_fact && std::abs(t.eval()) > thresh); return s; From abe766412513131cb8cb12d4f0c062fd3bcdbcae Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 16:27:01 -0500 Subject: [PATCH 136/137] More clarifications. --- .../libtwofloat/arithmetics/double-word-arithmetic.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index 0daa595..c515cf2 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -318,7 +318,7 @@ inline T two_sqr(T input, T &err) { } // Reference QD / inline.h (within qd namespace) -// TODO: make sure this has not been already implemented +// 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; @@ -328,7 +328,7 @@ inline T two_sum(T a, T b, T &err) { } // Reference: QD / dd_inline.h (within dd_real namespace) -// TODO: make sure this has not been already implemented +// TODO: make sure this has not been already implemented in twofloat template inline two add(T a, T b) { T s, e; @@ -538,7 +538,7 @@ inline two sin(const two &input) { } // Approximately reducing modulo 2*pi - two z = nint(div(input, local_2pi)); // TODO: check mode chosen for div (accurate + fma) + two z = nint(div(input, local_2pi)); two r = sub

(input, mul(local_2pi, z)); /* @@ -550,8 +550,6 @@ inline two sin(const two &input) { */ // Approximately reducing modulo pi/2 and then modulo pi/16 - - //TODO: original type is double, here it is templated T q = std::floor(r.h / local_pi2.h + pointfive); two t = sub

(r, mul(local_pi2, q)); int j = static_cast(q); From f1811d896350f681b4055b9813c9b666d1cfcce3 Mon Sep 17 00:00:00 2001 From: Leonardo Solis Vasquez Date: Wed, 7 Feb 2024 16:41:31 -0500 Subject: [PATCH 137/137] Documenting with /// \brief as much as possible. --- .../arithmetics/double-word-arithmetic.hpp | 89 +++++++++---------- 1 file changed, 44 insertions(+), 45 deletions(-) diff --git a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp index c515cf2..4aa4e62 100644 --- a/include/libtwofloat/arithmetics/double-word-arithmetic.hpp +++ b/include/libtwofloat/arithmetics/double-word-arithmetic.hpp @@ -234,7 +234,8 @@ inline two div(const two &x, const two &y) { static_assert(sizeof(T) == 0, "Unsupported mode"); } -// Reference: QD / dd_const.cpp / inline.h / dd_real.cpp +/// \brief QD constants +/// Reference: QD / dd_const.cpp / inline.h / dd_real.cpp template static const two _2pi = {6.283185307179586232e+00, 2.449293598294706414e-16}; @@ -250,7 +251,9 @@ static const constexpr T _nan = std::numeric_limits::quiet_NaN(); template static const constexpr T _eps = 4.93038065763132e-32; // 2^-104 -// Some recurrently-used contants +static const int n_inv_fact = 15; + +/// \brief Some recurrently-used contants template static const constexpr T pointfive = 0.5; @@ -263,9 +266,7 @@ static const constexpr T zeropointzero = 0.0; template static const constexpr T onepointzero = 1.0; -// Reference: QD / dd_real.cpp -static const int n_inv_fact = 15; - +/// \brief Tables containing constant data for trigonometric functions template struct constants_trig_tables { @@ -288,7 +289,7 @@ struct constants_trig_tables { }; // Reference: QD / dd_real.cpp - /* Table of sin(k * pi/16) and cos(k * pi/16). */ + // 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}, @@ -303,8 +304,8 @@ struct constants_trig_tables { }; }; -// Reference: QD / inline.h (within qd namespace) -/* Computes fl(input * input) and err(input * input) */ +/// \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"); @@ -317,8 +318,8 @@ inline T two_sqr(T input, T &err) { return q; } -// Reference QD / inline.h (within qd namespace) -// TODO: make sure this has not been already implemented in twofloat +/// 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; @@ -327,8 +328,8 @@ inline T two_sum(T a, T b, T &err) { return s; } -// Reference: QD / dd_inline.h (within dd_real namespace) -// TODO: make sure this has not been already implemented in twofloat +/// 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; @@ -336,8 +337,8 @@ inline two add(T a, T b) { return two (s, e); } -// Reference: QD / inline.h (within qd namespace) -/* Computes the nearest integer to input. */ +/// \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"); @@ -347,29 +348,28 @@ inline T nint(T input) { return std::floor(input + pointfive); } -// Reference: QD / dd_inline.h -/* Round to nearest integer */ +/// \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. */ + // 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. */ + // 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 */ + // 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. */ + // There is a tie in the high word, consult the low word to break the tie hi -= onepointzero; /* NOTE: This does not cause INEXACT. */ } } @@ -377,7 +377,7 @@ inline two nint(const two &input) { return two{hi, lo}; } -// Reference: QD / dd_inline.h (within dd_real namespace) +/// Reference: QD / dd_inline.h (within dd_real namespace) template inline two sqr(T input) { T p1, p2; @@ -385,7 +385,7 @@ inline two sqr(T input) { return two{p1, p2}; } -// Reference: QD / dd_inline.h +/// 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"); @@ -400,9 +400,9 @@ inline two sqr(const two &input) { return two{s1, s2}; } -// Reference: QD / dd_real.cpp -/* Computes sin(a) using Taylor series. - Assumes |a| <= pi/32 */ +/// \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"); @@ -432,14 +432,14 @@ static two sin_taylor(const two &input) { return s; } -// Reference: QD / dd_inline.h -/* double-double * double, where double is a power of 2. */ +/// \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 +/// 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"); @@ -468,21 +468,21 @@ static two cos_taylor(const two &input) { return s; } -// Reference: QD / dd_real.cpp -/* Computes the square root of the double-double number dd. - NOTE: dd must be a non-negative number. */ +/// \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) { - /* 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. - */ - static_assert((std::is_same_v || std::is_same_v), "Other types not supported"); T local_nan = _nan; @@ -502,12 +502,11 @@ two sqrt(const two &input) { return add(ax, temp.h * x * pointfive); } -// Reference: QD / dd_real.cpp +/// 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"); - // Reference: QD / dd_inline.h if (input.eval() == zeropointzero) { sin_a.h = zeropointzero; sin_a.l = zeropointzero; @@ -520,7 +519,7 @@ static void sincos_taylor(const two &input, two &sin_a, two &cos_a) { cos_a = sqrt

(sub(onepointzero, sqr(sin_a))); } -// Reference: QD / dd_real.cpp +/// 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");