From a7b6c5f9b51c879b222bc2985dbe7638268bc45f Mon Sep 17 00:00:00 2001 From: ZERICO2005 <71151164+ZERICO2005@users.noreply.github.com> Date: Fri, 9 Aug 2024 13:09:32 -0600 Subject: [PATCH] v1.2.2 | 2024/08/09 13:09 | Added more precision options to ABS-Fractal explorer, increasing the maximum zoom along with improving performance. --- .gitignore | 5 +- CMakeLists.txt | 22 +- docs/Build-Instructions.md | 10 +- docs/Libraries-Used.md | 13 +- docs/README.txt | 8 +- docs/Update-Log.txt | 5 + src/Program_Def.cpp | 14 +- src/Program_Def.h | 6 +- src/floats/double_Float32.hpp | 246 ++++--- src/floats/double_Float64.hpp | 240 ++++--- src/floats/double_Float64_AVX.h | 478 ++++++++++++- src/floats/double_Float64_SSE2.h | 708 -------------------- src/floats/double_Float80.hpp | 211 ++++-- src/floats/double_FloatN_snprintf.hpp | 75 ++- src/floats/double_FloatN_stringTo.hpp | 9 +- src/render_CPU/frac_Multi_Float64x2_AVX.cpp | 132 ++-- 16 files changed, 1120 insertions(+), 1062 deletions(-) delete mode 100644 src/floats/double_Float64_SSE2.h diff --git a/.gitignore b/.gitignore index ae3de2c..4f642e8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,7 @@ desktop.ini .vscode +.vs imgui.ini -.cache \ No newline at end of file +.cache +out +CMakeSettings.json \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 7820076..4a33872 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,10 +18,10 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "./bin") # GCC quadmath __float128 set(FEATURE_Float128 "") - # Requires libqd.a + # (Deprecated) Requires libqd.a # set(FEATURE_LIBQD "") - # MPFR converts between dekker-floats and strings + # MPFR converts between coordinates and strings set(FEATURE_FloatMPFR "") # GPU Rendering (OpenCL 1.2 or later) @@ -35,9 +35,6 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "./bin") # Requires libquadmath.a # set(STATIC_OpenCV_Scaler "") - # Defines the NDEBUG macro to disable -# target_compile_definitions(${PROJECT_NAME} PRIVATE NDEBUG) - # Long compile times! Optimizes and reduces the size of the binary for release # set(COMPILE_FLTO "") @@ -54,7 +51,7 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "./bin") set(LIBRARY_OPT_FLAG -O3) # -Ofast appears to have a noticable speedup compared to -O3 for SSE2 and AVX rendering. set(RENDER_CPU_OPT_FLAG -Ofast) - set(MARCH_FLAGS -msse2) + set(MARCH_FLAGS -msse3) # Source Files file(GLOB_RECURSE SRC_FILES @@ -70,11 +67,13 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "./bin") endif() - - file(GLOB SRC_AVX_FILES "${SRC_DIR}/render_CPU/frac_Multi_AVX.cpp") - set_source_files_properties(SOURCE ${SRC_AVX_FILES} PROPERTY COMPILE_FLAGS "-mavx ${RENDER_CPU_OPT_FLAG}") file(GLOB SRC_AVX_FILES "${SRC_DIR}/render_CPU/frac_Multi_SSE2.cpp") set_source_files_properties(SOURCE ${SRC_AVX_FILES} PROPERTY COMPILE_FLAGS "${RENDER_CPU_OPT_FLAG}") + file(GLOB SRC_AVX_FILES "${SRC_DIR}/render_CPU/frac_Multi_AVX.cpp") + set_source_files_properties(SOURCE ${SRC_AVX_FILES} PROPERTY COMPILE_FLAGS "-mavx ${RENDER_CPU_OPT_FLAG}") + file(GLOB SRC_AVX_FILES "${SRC_DIR}/render_CPU/frac_Multi_Float64x2_AVX.cpp") + set_source_files_properties(SOURCE ${SRC_AVX_FILES} PROPERTY COMPILE_FLAGS "-mavx ${PROGRAM_OPT_FLAG}") + # Icon set(ICON_FILE "./icons/ABS-Icon.o") @@ -89,6 +88,9 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "./bin") # Packages and Macros find_package(SDL2 REQUIRED) + # Defines the NDEBUG macro to disable +# target_compile_definitions(${PROJECT_NAME} PRIVATE NDEBUG) + if(DEFINED FEATURE_Float80) target_compile_definitions(${PROJECT_NAME} PRIVATE Enable_Float80) endif() @@ -170,8 +172,6 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "./bin") # Compile Options - set(MARCH_FLAGS -msse3 -mavx) - set(WARNING_FLAGS -Wall -Wextra -Wshadow -Werror=infinite-recursion diff --git a/docs/Build-Instructions.md b/docs/Build-Instructions.md index 384ddb6..f8911cc 100644 --- a/docs/Build-Instructions.md +++ b/docs/Build-Instructions.md @@ -1,7 +1,7 @@ Build Instructions
-ABS-Fractal-Explorer v1.2.1 rev-2 +ABS-Fractal-Explorer v1.2.2 -Updated: 2024/06/09 +Updated: 2024/08/09 You can read this guide: https://github.com/ZERICO2005/ABS-Fractal-Explorer/wiki/How-to-Compile @@ -16,9 +16,15 @@ Enables 80bit floats (`long double` or `__float80`) for extra precision. Used to ## Float128 Enables 128bit floats (quadmath.h) for extra precision. Used to store cordinates if it is the highest available precision. +## libqd +(Deprecated) Experimental double-double/quad-double float support. May require some manual configuration to enable it as a rendering mode. + ## OpenCL Enables OpenCL GPU rendering +## MPFR +Used to convert between coordinates and strings + ## WindowsFileIO (Not implemented in v1.2.1 rev-2)
Uses the Windows file dialogs instead of the generic ones diff --git a/docs/Libraries-Used.md b/docs/Libraries-Used.md index 8d7fb19..fcca8a6 100644 --- a/docs/Libraries-Used.md +++ b/docs/Libraries-Used.md @@ -50,11 +50,22 @@ Used for GPU rendering https://gcc.gnu.org/onlinedocs/libquadmath/ Used for high precision 128bit float calculations +## MPFR +Used for converting between coordinates and strings with greater accuracy. + +## libQD +https://www.davidhbailey.com/dhbsoftware/ +Used for experimental quad-double and double-double arithmetic support. + +## LIB-Dekker-Float +https://github.com/ZERICO2005/LIB-Dekker-Float +A header-only library I wrote for the Dekker float implementations used in ABS-Fractal-Explorer. Float64x2 is roughly the same speed as dd_real in libQD, while Float64x2 AVX is ~3.3x faster for rendering fractals. + ## STB Image https://github.com/nothings/stb Used for writing PNG/JPG images ## OpenCV https://opencv.org/ -Used for frame scaling in custom versions of ABS-Fractal-Explorer +Used for frame scaling in some versions of ABS-Fractal-Explorer diff --git a/docs/README.txt b/docs/README.txt index 94b635b..3e432bb 100644 --- a/docs/README.txt +++ b/docs/README.txt @@ -1,5 +1,5 @@ -ABS-Fractal-Explorer v1.2.1 (Revision-2) (Windows 10 Edition) -zerico2005 | 2024/06/09 +ABS-Fractal-Explorer v1.2.2 (Windows 10 Edition) +zerico2005 | 2024/08/09 ABS-Fractal-Explorer will allow you to explore hundereds of different Mandelbrot variants. From Quadratic to Sextic fractals, there is always something new and novel to explore! @@ -30,7 +30,9 @@ Controls (Full controls/keybinds are in the Key-binds Menu): ;": Inc/dec supersampling (Renders at higher quality but is much slower) M : fp32 GPU rendering (10^5.7 max zoom, fastest) N : fp64 CPU rendering (10^14.4 max zoom) - B : fp80 CPU rendering (10^17.7 max zoom) + B : fp80 CPU rendering (10^17.8 max zoom) + Rctrl: fp64x2 CPU rendering (10^30.7 max zoom) + View the rendering menu for more precision options. Julia-Sets: C: Toggle Julia Sets V: Toggle Starting Z-Value diff --git a/docs/Update-Log.txt b/docs/Update-Log.txt index 11b4a14..0717a85 100644 --- a/docs/Update-Log.txt +++ b/docs/Update-Log.txt @@ -1,5 +1,10 @@ ABS-Fractal-Explorer Update-Log +# v1.2.2 (2024/08/09) +* Added Dekker floats, increasing the maximum precision from 10^32.5 to 10^36.7 +* Precisions added: Float32x2 (10^12.6), Float64x2 (10^30.7), Float80x2 (10^36.7) +* General improvements and bug fixes. + # v1.2.1 (Revision-2) (2024/06/09) * Corrected color averaging algorithm diff --git a/src/Program_Def.cpp b/src/Program_Def.cpp index 19a3c03..3a0f9f2 100644 --- a/src/Program_Def.cpp +++ b/src/Program_Def.cpp @@ -140,10 +140,12 @@ #else int FloatCoordinate_snprintf(char* buf, size_t len, const char* format, fpCord cord) { - #if defined(Enable_Float128) + #if defined(Enable_Float80) + return Float80x2_snprintf(buf, len, format, cord); + #elif defined(Enable_Float128) return quadmath_snprintf(buf, len, format, (fp128)cord); #else - return snprintf(buf, len, format, cord); + return Float64x2_snprintf(buf, len, format, cord); #endif } @@ -174,12 +176,12 @@ } fpCord stringTo_FloatCoordinate(const char* nPtr, char** endPtr) { - #if defined(Enable_Float128) + #if defined(Enable_Float80) + return (fpCord)stringTo_Float80x2(nPtr, endPtr); + #elif defined(Enable_Float128) return (fpCord)stringTo_Float128(nPtr, endPtr); - #elif defined(Enable_Float80) - return (fpCord)stringTo_Float80(nPtr, endPtr); #else - return (fpCord)stringTo_Float64(nPtr, endPtr); + return (fpCord)stringTo_Float64x2(nPtr, endPtr); #endif } diff --git a/src/Program_Def.h b/src/Program_Def.h index e0b463b..5f16ffe 100644 --- a/src/Program_Def.h +++ b/src/Program_Def.h @@ -37,11 +37,11 @@ /* Version */ #define PROGRAM_NAME "ABS-Fractal-Explorer" -#define PROGRAM_DATE "2024/08/08" /* YYYY/MM/DD */ +#define PROGRAM_DATE "2024/08/09" /* YYYY/MM/DD */ #define PROGRAM_V_MAJOR 1 #define PROGRAM_V_MINOR 2 -#define PROGRAM_V_PATCH 1 -#define PROGRAM_V_TAG "(Revision-12)" +#define PROGRAM_V_PATCH 2 +#define PROGRAM_V_TAG "" #define PROGRAM_VERSION STR_N(PROGRAM_V_MAJOR) "." STR_N(PROGRAM_V_MINOR) "." STR_N(PROGRAM_V_PATCH) " " PROGRAM_V_TAG /* Float80 and Float128 */ diff --git a/src/floats/double_Float32.hpp b/src/floats/double_Float32.hpp index 7bc182d..aa1df05 100644 --- a/src/floats/double_Float32.hpp +++ b/src/floats/double_Float32.hpp @@ -1,6 +1,6 @@ /* -** Author: zerico2005 (2023 - 2024) -** Project: ABS-Fractal-Explorer +** Author: zerico2005 (2024) +** Project: LIB-Dekker-Float ** License: MIT License ** A copy of the MIT License should be included with ** this project. If not, see https://opensource.org/license/MIT @@ -19,26 +19,27 @@ typedef float fp32; typedef double fp64; +// Can be changed to other types for better accuracy typedef fp64 fp32x2_Math; /** * @brief Double-Float32 Dekker Float implementation. * Source: Creel "Double it Like Dekker" on YouTube. + * + * @warning -Ofast may break this library. -O3 compiles okay on gcc and clang. */ -class Float32x2 { -public: +struct Float32x2 { + fp32 hi; fp32 lo; -private: - /* Arithmetic */ - inline Float32x2 Dekker_Add( + static inline Float32x2 Dekker_Add( const Float32x2& x, const Float32x2& y - ) const { + ) { fp32 r_hi = x.hi + y.hi; - fp32 r_lo = 0.0f; + fp32 r_lo = static_cast(0.0); if (fabsf(x.hi) > fabsf(y.hi)) { r_lo = x.hi - r_hi + y.hi + y.lo + x.lo; } else { @@ -51,11 +52,11 @@ class Float32x2 { return c; } - inline Float32x2 Dekker_Sub( + static inline Float32x2 Dekker_Sub( const Float32x2& x, const Float32x2& y - ) const { + ) { fp32 r_hi = x.hi - y.hi; - fp32 r_lo = 0.0f; + fp32 r_lo = static_cast(0.0); if (fabsf(x.hi) > fabsf(y.hi)) { r_lo = x.hi - r_hi - y.hi - y.lo + x.lo; } else { @@ -70,7 +71,7 @@ class Float32x2 { static constexpr fp32 Dekker_Scale = 4097.0f; // (2^ceil(24 / 2) + 1) - inline Float32x2 Dekker_Split(const fp32& x) const { + static inline Float32x2 Dekker_Split(const fp32& x) { fp32 p = x * Dekker_Scale; Float32x2 r; r.hi = (x - p) + p; @@ -80,7 +81,7 @@ class Float32x2 { // static constexpr uint64_t Dekker_Split_Mask = ~((uint64_t)0xFFF); - // inline Float32x2 Dekker_Split(const fp32& x) const { + // static inline Float32x2 Dekker_Split(const fp32& x) { // Float32x2 r; // uint64_t temp = (*(uint64_t*)((void*)&x)) & Dekker_Split_Mask; // r.hi = (*(fp32*)((void*)&temp)); @@ -88,9 +89,9 @@ class Float32x2 { // return r; // } - inline Float32x2 Dekker_Mul12( + static inline Float32x2 Dekker_Mul12( const fp32& x, const fp32& y - ) const { + ) { Float32x2 a = Dekker_Split(x); Float32x2 b = Dekker_Split(y); fp32 p = a.hi * b.hi; @@ -102,9 +103,9 @@ class Float32x2 { return r; } - inline Float32x2 Dekker_Mul( + static inline Float32x2 Dekker_Mul( const Float32x2& x, const Float32x2& y - ) const { + ) { Float32x2 t = Dekker_Mul12(x.hi, y.hi); fp32 c = x.hi * y.lo + x.lo * y.hi + t.lo; @@ -114,9 +115,9 @@ class Float32x2 { return r; } - inline Float32x2 Dekker_Div( + static inline Float32x2 Dekker_Div( const Float32x2& x, const Float32x2& y - ) const { + ) { Float32x2 u; u.hi = x.hi / y.hi; Float32x2 t = Dekker_Mul12(u.hi, y.hi); @@ -128,13 +129,38 @@ class Float32x2 { return r; } + static inline Float32x2 Dekker_Sqr12( + const fp32& x + ) { + Float32x2 a = Dekker_Split(x); + fp32 p = a.hi * a.hi; + fp32 q = static_cast(2.0) * a.hi * a.lo; + + Float32x2 r; + r.hi = p + q; + r.lo = p - r.hi + q + a.lo * a.lo; + return r; + } + + static inline Float32x2 Dekker_Sqr( + const Float32x2& x + ) { + Float32x2 t = Dekker_Sqr12(x.hi); + fp32 c = static_cast(2.0) * x.hi * x.lo + t.lo; + + Float32x2 r; + r.hi = t.hi + c; + r.lo = t.hi - r.hi + c; + return r; + } + /* Double-Single Arithmetic */ - inline Float32x2 Dekker_Add_Float32( + static inline Float32x2 Dekker_Add_Float32( const Float32x2& x, const fp32& y - ) const { + ) { fp32 r_hi = x.hi + y; - fp32 r_lo = 0.0f; + fp32 r_lo = static_cast(0.0); if (fabsf(x.hi) > fabsf(y)) { r_lo = x.hi - r_hi + y + x.lo; } else { @@ -147,11 +173,11 @@ class Float32x2 { return c; } - inline Float32x2 Dekker_Sub_Float32( + static inline Float32x2 Dekker_Sub_Float32( const Float32x2& x, const fp32& y - ) const { + ) { fp32 r_hi = x.hi - y; - fp32 r_lo = 0.0f; + fp32 r_lo = static_cast(0.0); if (fabsf(x.hi) > fabsf(y)) { r_lo = x.hi - r_hi - y + x.lo; } else { @@ -164,9 +190,9 @@ class Float32x2 { return c; } - inline Float32x2 Dekker_Mul_Float32( + static inline Float32x2 Dekker_Mul_Float32( const Float32x2& x, const fp32& y - ) const { + ) { Float32x2 t = Dekker_Mul12(x.hi, y); fp32 c = x.lo * y + t.lo; @@ -176,9 +202,9 @@ class Float32x2 { return r; } - inline Float32x2 Dekker_Div_Float32( + static inline Float32x2 Dekker_Div_Float32( const Float32x2& x, const fp32& y - ) const { + ) { Float32x2 u; u.hi = x.hi / y; Float32x2 t = Dekker_Mul12(u.hi, y); @@ -190,9 +216,9 @@ class Float32x2 { return r; } - inline Float32x2 Float32_Div_Dekker( + static inline Float32x2 Float32_Div_Dekker( const fp32& x, const Float32x2& y - ) const { + ) { Float32x2 u; u.hi = x / y.hi; Float32x2 t = Dekker_Mul12(u.hi, y.hi); @@ -204,8 +230,6 @@ class Float32x2 { return r; } -public: - /* Arithmetic */ inline Float32x2 operator+(const Float32x2 &value) const { @@ -360,26 +384,17 @@ class Float32x2 { /* Constructors */ - inline Float32x2() { - // this->hi = 0.0; - // this->lo = 0.0; - } + constexpr inline Float32x2() : hi(), lo() {} - inline Float32x2(const fp32& value) { - this->hi = value; - this->lo = 0.0; - } + constexpr inline Float32x2(const fp32& value) : + hi(value), lo(0.0f) {} - inline Float32x2(const fp64& value) { - this->hi = (fp32)value; - this->lo = (fp32)(value - (fp64)this->hi); - } + constexpr inline Float32x2(const fp64& value) : + hi((fp32)value), lo((fp32)(value - (fp64)this->hi)) {} template - inline Float32x2(const fpX& value) { - this->hi = (fp32)value; - this->lo = (fp32)(value - (fpX)this->hi); - } + constexpr inline Float32x2(const fpX& value) : + hi((fp32)value), lo((fp32)(value - (fpX)this->hi)) {} /* Casts */ @@ -402,12 +417,23 @@ typedef Float32x2 fp32x2; /* Math functions (Natively implemented) */ /* Arithmetic */ + inline fp32x2 fmax(fp32x2 x, fp32x2 y) { return (x > y) ? x : y; } + // inline fp32x2 fmax(fp32x2 x, fp32x2 y, fp32x2 z) { + // return (x > y) ? + // ((x > z) ? x : z) : + // ((y > z) ? y : z); + // } inline fp32x2 fmin(fp32x2 x, fp32x2 y) { return (x < y) ? x : y; } + // inline fp32x2 fmin(fp32x2 x, fp32x2 y, fp32x2 z) { + // return (x < y) ? + // ((x < z) ? x : z) : + // ((y < z) ? y : z); + // } inline fp32x2 fabs(fp32x2 x) { return (x < static_cast(0.0)) ? -x : x; } @@ -419,13 +445,44 @@ typedef Float32x2 fp32x2; } inline fp32x2 copysign(fp32x2 x, fp32x2 y) { return ( - (x < static_cast(0.0)) != (y < (static_cast(0.0))) + (x.hi < static_cast(0.0)) != (y.hi < (static_cast(0.0))) ) ? -x : x; } + /** @note This function name may change to square() or etc */ + inline fp32x2 sqr(fp32x2 x) { + return Float32x2::Dekker_Sqr(x); + } + inline fp32x2 sqrt(fp32x2 x) { + if (x == static_cast(0.0)) { + return x; + } + fp32x2 guess = (fp32x2)sqrt(x.hi); + return (guess + x / guess) * static_cast(0.5); + } + inline fp32x2 cbrt(fp32x2 x) { + if (x == static_cast(0.0)) { + return x; + } + fp32x2 guess = (fp32x2)cbrt(x.hi); + return ( + guess * static_cast(2.0) + (x) / Float32x2::Dekker_Sqr(guess) + ) / static_cast(3.0); + } + inline fp32x2 hypot(fp32x2 x, fp32x2 y) { + return sqrt( + Float32x2::Dekker_Sqr(x) + Float32x2::Dekker_Sqr(y) + ); + } + // inline fp32x2 hypot(fp32x2 x, fp32x2 y, fp32x2 z) { + // return sqrt( + // Float32x2::Dekker_Sqr(x) + Float32x2::Dekker_Sqr(y) + Float32x2::Dekker_Sqr(z) + // ); + // } /* Tests */ + inline bool signbit(fp32x2 x) { - return (x < static_cast(0.0)) ? true : false; + return (x.hi < static_cast(0.0)) ? true : false; } /** Returns true if both x.hi and x.lo are finite */ inline bool isfinite(fp32x2 x) { @@ -456,6 +513,7 @@ typedef Float32x2 fp32x2; } /* Comparison */ + inline bool isgreater(fp32x2 x, fp32x2 y) { return (x > y); } @@ -472,26 +530,27 @@ typedef Float32x2 fp32x2; return (x < y) || (x > y); } - /* Rounding */ - inline fp32x2 trunc(fp32x2 x) { - fp32 frac_hi = x.hi - trunc(x.hi); - fp32 frac_lo = x.lo - trunc(x.lo); - fp32x2 int_hi = trunc(x.hi); - fp32x2 int_lo = trunc(x.lo); - // Sum in increasing order - fp32x2 trunc_all = static_cast(0.0); - trunc_all += ( + /* Rounding */ + + inline fp32x2 trunc(fp32x2 x) { + fp32x2 int_hi = trunc(x.hi); + fp32x2 int_lo = trunc(x.lo); + fp32 frac_hi = x.hi - int_hi.hi; + fp32 frac_lo = x.lo - int_lo.hi; + // Sum in increasing order + fp32x2 trunc_all = static_cast(0.0); + trunc_all += ( (fp32x2)frac_hi + (fp32x2)frac_lo >= static_cast(1.0) ) ? static_cast(1.0) : static_cast(0.0); - trunc_all += int_lo; - trunc_all += int_hi; - return trunc_all; - } + trunc_all += int_lo; + trunc_all += int_hi; + return trunc_all; + } inline fp32x2 floor(fp32x2 x) { fp32x2 int_part = trunc(x); return ( x < static_cast(0.0) && int_part != x - ) ? int_part : int_part - static_cast(1.0); + ) ? int_part - static_cast(1.0) : int_part; } inline fp32x2 ceil(fp32x2 x) { fp32x2 int_part = trunc(x); @@ -540,10 +599,36 @@ typedef Float32x2 fp32x2; } /* Integer and Remainder */ + + inline fp32x2 modf(fp32x2 x, fp32x2* int_part) { + fp32x2 trunc_part = trunc(x); + if (int_part != nullptr) { + *int_part = trunc_part; + } + return x - trunc_part; + } inline fp32x2 nearbyint(fp32x2 x) { return rint(x); } + /* Float Exponents */ + + inline fp32x2 ldexp(fp32x2 x, int exp) { + x.hi = ldexp(x.hi, exp); + x.lo = isfinite(x.hi) ? ldexp(x.lo, exp) : x.hi; + return x; + } + inline fp32x2 scalbn(fp32x2 x, int exp) { + x.hi = scalbn(x.hi, exp); + x.lo = isfinite(x.hi) ? scalbn(x.lo, exp) : x.hi; + return x; + } + inline fp32x2 scalbln(fp32x2 x, long exp) { + x.hi = scalbln(x.hi, exp); + x.lo = isfinite(x.hi) ? scalbln(x.lo, exp) : x.hi; + return x; + } + /* Math overloads (Casts to other types) */ /* Arithmetic */ @@ -553,9 +638,9 @@ typedef Float32x2 fp32x2; // inline fp32x2 fdim(fp32x2 x, fp32x2 y) { return (fp32x2)fdim((fp32x2_Math)x, (fp32x2_Math)y); } // inline fp32x2 fma(fp32x2 x, fp32x2 y, fp32x2 z) { return (fp32x2)fma((fp32x2_Math)x, (fp32x2_Math)y, (fp32x2_Math)z); } // inline fp32x2 copysign(fp32x2 x, fp32x2 y) { return (fp32x2)copysign((fp32x2_Math)x, (fp32x2_Math)y); } - inline fp32x2 sqrt(fp32x2 x) { return (fp32x2)sqrt((fp32x2_Math)x); } - inline fp32x2 cbrt(fp32x2 x) { return (fp32x2)cbrt((fp32x2_Math)x); } - inline fp32x2 hypot(fp32x2 x, fp32x2 y) { return (fp32x2)hypot((fp32x2_Math)x, (fp32x2_Math)y); } + // inline fp32x2 sqrt(fp32x2 x) { return (fp32x2)sqrt((fp32x2_Math)x); } + // inline fp32x2 cbrt(fp32x2 x) { return (fp32x2)cbrt((fp32x2_Math)x); } + // inline fp32x2 hypot(fp32x2 x, fp32x2 y) { return (fp32x2)hypot((fp32x2_Math)x, (fp32x2_Math)y); } /* Trigonometry */ inline fp32x2 sin (fp32x2 x) { return (fp32x2) sin ((fp32x2_Math)x); } inline fp32x2 cos (fp32x2 x) { return (fp32x2) cos ((fp32x2_Math)x); } @@ -574,6 +659,10 @@ typedef Float32x2 fp32x2; *p_sin = sin(x); *p_cos = cos(x); } + inline void sinhcosh(fp32x2 x, fp32x2* p_sinh, fp32x2* p_cosh) { + *p_sinh = sinh(x); + *p_cosh = cosh(x); + } /* Logarithms and Exponents */ inline fp32x2 log (fp32x2 x) { return (fp32x2)log ((fp32x2_Math)x); } inline fp32x2 log1p(fp32x2 x) { return (fp32x2)log1p((fp32x2_Math)x); } @@ -583,6 +672,7 @@ typedef Float32x2 fp32x2; inline fp32x2 exp (fp32x2 x) { return (fp32x2)exp ((fp32x2_Math)x); } inline fp32x2 expm1(fp32x2 x) { return (fp32x2)expm1((fp32x2_Math)x); } inline fp32x2 exp2 (fp32x2 x) { return (fp32x2)exp2 ((fp32x2_Math)x); } + inline fp32x2 exp10(fp32x2 x) { return (fp32x2)exp((fp32x2_Math)x * (fp32x2_Math)log(static_cast(10.0))); } inline fp32x2 pow(fp32x2 x, fp32x2 y) { return (fp32x2)pow((fp32x2_Math)x, (fp32x2_Math)y); } /* Rounding */ // inline fp32x2 trunc(fp32x2 x) { return (fp32x2)trunc((fp32x2_Math)x); } @@ -596,12 +686,12 @@ typedef Float32x2 fp32x2; // inline long long llround(fp32x2 x) { return llround((fp32x2_Math)x); } /* Integer and Remainder */ inline fp32x2 fmod(fp32x2 x, fp32x2 y) { return (fp32x2)fmod((fp32x2_Math)x, (fp32x2_Math)y); } - inline fp32x2 modf(fp32x2 x, fp32x2* y) { - fp32x2_Math y_temp; - fp32x2 result = modf((fp32x2_Math)x, &y_temp); - *y = (fp32x2)y_temp; - return result; - } + // inline fp32x2 modf(fp32x2 x, fp32x2* y) { + // fp32x2_Math y_temp; + // fp32x2 result = modf((fp32x2_Math)x, &y_temp); + // *y = (fp32x2)y_temp; + // return result; + // } // inline fp32x2 nearbyint(fp32x2 x) { return (fp32x2)nearbyint((fp32x2_Math)x); } // Incorrect Function // inline fp32x2 nextafter(fp32x2 x) { return (fp32x2)nextafter((fp32x2_Math)x); } inline fp32x2 remainder(fp32x2 x, fp32x2 y) { return (fp32x2)remainder((fp32x2_Math)x, (fp32x2_Math)y); } @@ -609,9 +699,9 @@ typedef Float32x2 fp32x2; /* Float Exponents */ inline int ilogb(fp32x2 x) { return ilogb((fp32x2_Math)x); } inline fp32x2 frexp (fp32x2 x, int* exp) { return (fp32x2)frexp ((fp32x2_Math)x, exp); } - inline fp32x2 ldexp (fp32x2 x, int exp) { return (fp32x2)ldexp ((fp32x2_Math)x, exp); } - inline fp32x2 scalbn (fp32x2 x, int exp) { return (fp32x2)scalbn ((fp32x2_Math)x, exp); } - inline fp32x2 scalbln(fp32x2 x, long exp) { return (fp32x2)scalbln((fp32x2_Math)x, exp); } + // inline fp32x2 ldexp (fp32x2 x, int exp) { return (fp32x2)ldexp ((fp32x2_Math)x, exp); } + // inline fp32x2 scalbn (fp32x2 x, int exp) { return (fp32x2)scalbn ((fp32x2_Math)x, exp); } + // inline fp32x2 scalbln(fp32x2 x, long exp) { return (fp32x2)scalbln((fp32x2_Math)x, exp); } /* Tests */ // inline bool signbit(fp32x2 x) { return (signbit((fp32x2_Math)x) != 0) ? true : false; } // inline bool isfinite(fp32x2 x) { return (isfinite((fp32x2_Math)x) != 0) ? true : false; } diff --git a/src/floats/double_Float64.hpp b/src/floats/double_Float64.hpp index 0d3ec1d..e3a7ffa 100644 --- a/src/floats/double_Float64.hpp +++ b/src/floats/double_Float64.hpp @@ -1,6 +1,6 @@ /* -** Author: zerico2005 (2023 - 2024) -** Project: ABS-Fractal-Explorer +** Author: zerico2005 (2024) +** Project: LIB-Dekker-Float ** License: MIT License ** A copy of the MIT License should be included with ** this project. If not, see https://opensource.org/license/MIT @@ -16,6 +16,7 @@ typedef float fp32; typedef double fp64; +// Can be changed to other types for better accuracy #if defined(Enable_Float128) #include "Float128.hpp" typedef fp128 fp64x2_Math; @@ -26,25 +27,23 @@ typedef double fp64; typedef fp64 fp64x2_Math; #endif - /** * @brief Double-Float64 Dekker Float implementation. * Source: Creel "Double it Like Dekker" on YouTube. + * + * @warning -Ofast may break this library. -O3 compiles okay on gcc and clang. */ -class Float64x2 { -public: +struct Float64x2 { fp64 hi; fp64 lo; - -private: /* Arithmetic */ - inline Float64x2 Dekker_Add( + static inline Float64x2 Dekker_Add( const Float64x2& x, const Float64x2& y - ) const { + ) { fp64 r_hi = x.hi + y.hi; - fp64 r_lo = 0.0; + fp64 r_lo = static_cast(0.0); if (fabs(x.hi) > fabs(y.hi)) { r_lo = x.hi - r_hi + y.hi + y.lo + x.lo; } else { @@ -57,11 +56,11 @@ class Float64x2 { return c; } - inline Float64x2 Dekker_Sub( + static inline Float64x2 Dekker_Sub( const Float64x2& x, const Float64x2& y - ) const { + ) { fp64 r_hi = x.hi - y.hi; - fp64 r_lo = 0.0; + fp64 r_lo = static_cast(0.0); if (fabs(x.hi) > fabs(y.hi)) { r_lo = x.hi - r_hi - y.hi - y.lo + x.lo; } else { @@ -76,7 +75,7 @@ class Float64x2 { static constexpr fp64 Dekker_Scale = 134217729.0; // (2^ceil(53 / 2) + 1) - inline Float64x2 Dekker_Split(const fp64& x) const { + static inline Float64x2 Dekker_Split(const fp64& x) { fp64 p = x * Dekker_Scale; Float64x2 r; r.hi = (x - p) + p; @@ -86,7 +85,7 @@ class Float64x2 { // static constexpr uint64_t Dekker_Split_Mask = ~((uint64_t)0x3FFFFFF); - // inline Float64x2 Dekker_Split(const fp64& x) const { + // static inline Float64x2 Dekker_Split(const fp64& x) { // Float64x2 r; // uint64_t temp = (*(uint64_t*)((void*)&x)) & Dekker_Split_Mask; // r.hi = (*(fp64*)((void*)&temp)); @@ -94,9 +93,9 @@ class Float64x2 { // return r; // } - inline Float64x2 Dekker_Mul12( + static inline Float64x2 Dekker_Mul12( const fp64& x, const fp64& y - ) const { + ) { Float64x2 a = Dekker_Split(x); Float64x2 b = Dekker_Split(y); fp64 p = a.hi * b.hi; @@ -108,9 +107,9 @@ class Float64x2 { return r; } - inline Float64x2 Dekker_Mul( + static inline Float64x2 Dekker_Mul( const Float64x2& x, const Float64x2& y - ) const { + ) { Float64x2 t = Dekker_Mul12(x.hi, y.hi); fp64 c = x.hi * y.lo + x.lo * y.hi + t.lo; @@ -120,9 +119,9 @@ class Float64x2 { return r; } - inline Float64x2 Dekker_Div( + static inline Float64x2 Dekker_Div( const Float64x2& x, const Float64x2& y - ) const { + ) { Float64x2 u; u.hi = x.hi / y.hi; Float64x2 t = Dekker_Mul12(u.hi, y.hi); @@ -134,13 +133,38 @@ class Float64x2 { return r; } + static inline Float64x2 Dekker_Sqr12( + const fp64& x + ) { + Float64x2 a = Dekker_Split(x); + fp64 p = a.hi * a.hi; + fp64 q = static_cast(2.0) * a.hi * a.lo; + + Float64x2 r; + r.hi = p + q; + r.lo = p - r.hi + q + a.lo * a.lo; + return r; + } + + static inline Float64x2 Dekker_Sqr( + const Float64x2& x + ) { + Float64x2 t = Dekker_Sqr12(x.hi); + fp64 c = static_cast(2.0) * x.hi * x.lo + t.lo; + + Float64x2 r; + r.hi = t.hi + c; + r.lo = t.hi - r.hi + c; + return r; + } + /* Double-Single Arithmetic */ - inline Float64x2 Dekker_Add_Float64( + static inline Float64x2 Dekker_Add_Float64( const Float64x2& x, const fp64& y - ) const { + ) { fp64 r_hi = x.hi + y; - fp64 r_lo = 0.0f; + fp64 r_lo = static_cast(0.0); if (fabs(x.hi) > fabs(y)) { r_lo = x.hi - r_hi + y + x.lo; } else { @@ -153,11 +177,11 @@ class Float64x2 { return c; } - inline Float64x2 Dekker_Sub_Float64( + static inline Float64x2 Dekker_Sub_Float64( const Float64x2& x, const fp64& y - ) const { + ) { fp64 r_hi = x.hi - y; - fp64 r_lo = 0.0f; + fp64 r_lo = static_cast(0.0); if (fabs(x.hi) > fabs(y)) { r_lo = x.hi - r_hi - y + x.lo; } else { @@ -170,9 +194,9 @@ class Float64x2 { return c; } - inline Float64x2 Dekker_Mul_Float64( + static inline Float64x2 Dekker_Mul_Float64( const Float64x2& x, const fp64& y - ) const { + ) { Float64x2 t = Dekker_Mul12(x.hi, y); fp64 c = x.lo * y + t.lo; @@ -182,9 +206,9 @@ class Float64x2 { return r; } - inline Float64x2 Dekker_Div_Float64( + static inline Float64x2 Dekker_Div_Float64( const Float64x2& x, const fp64& y - ) const { + ) { Float64x2 u; u.hi = x.hi / y; Float64x2 t = Dekker_Mul12(u.hi, y); @@ -196,9 +220,9 @@ class Float64x2 { return r; } - inline Float64x2 Float64_Div_Dekker( + static inline Float64x2 Float64_Div_Dekker( const fp64& x, const Float64x2& y - ) const { + ) { Float64x2 u; u.hi = x / y.hi; Float64x2 t = Dekker_Mul12(u.hi, y.hi); @@ -210,8 +234,6 @@ class Float64x2 { return r; } -public: - /* Arithmetic */ inline Float64x2 operator+(const Float64x2 &value) const { @@ -366,25 +388,17 @@ class Float64x2 { /* Constructors */ - inline Float64x2() { - // this->hi = 0.0; - // this->lo = 0.0; - } + constexpr inline Float64x2() : hi(), lo() {} - inline Float64x2(const fp32& value) { - this->hi = (fp64)value; - this->lo = 0.0; - } - inline Float64x2(const fp64& value) { - this->hi = value; - this->lo = 0.0; - } + constexpr inline Float64x2(const fp32& value) : + hi((fp64)value), lo(0.0) {} + + constexpr inline Float64x2(const fp64& value) : + hi(value), lo(0.0) {} template - inline Float64x2(const fpX& value) { - this->hi = (fp64)value; - this->lo = (fp64)(value - (fpX)this->hi); - } + constexpr inline Float64x2(const fpX& value) : + hi((fp64)value), lo((fp64)(value - (fpX)this->hi)) {} /* Casts */ @@ -407,12 +421,23 @@ typedef Float64x2 fp64x2; /* Math functions (Natively implemented) */ /* Arithmetic */ + inline fp64x2 fmax(fp64x2 x, fp64x2 y) { return (x > y) ? x : y; } + // inline fp64x2 fmax(fp64x2 x, fp64x2 y, fp64x2 z) { + // return (x > y) ? + // ((x > z) ? x : z) : + // ((y > z) ? y : z); + // } inline fp64x2 fmin(fp64x2 x, fp64x2 y) { return (x < y) ? x : y; } + // inline fp64x2 fmin(fp64x2 x, fp64x2 y, fp64x2 z) { + // return (x < y) ? + // ((x < z) ? x : z) : + // ((y < z) ? y : z); + // } inline fp64x2 fabs(fp64x2 x) { return (x < static_cast(0.0)) ? -x : x; } @@ -424,13 +449,44 @@ typedef Float64x2 fp64x2; } inline fp64x2 copysign(fp64x2 x, fp64x2 y) { return ( - (x < static_cast(0.0)) != (y < (static_cast(0.0))) + (x.hi < static_cast(0.0)) != (y.hi < (static_cast(0.0))) ) ? -x : x; } + /** @note This function name may change to square() or etc */ + inline fp64x2 sqr(fp64x2 x) { + return Float64x2::Dekker_Sqr(x); + } + inline fp64x2 sqrt(fp64x2 x) { + if (x == static_cast(0.0)) { + return x; + } + fp64x2 guess = (fp64x2)sqrt(x.hi); + return (guess + x / guess) * static_cast(0.5); + } + inline fp64x2 cbrt(fp64x2 x) { + if (x == static_cast(0.0)) { + return x; + } + fp64x2 guess = (fp64x2)cbrt(x.hi); + return ( + guess * static_cast(2.0) + (x) / Float64x2::Dekker_Sqr(guess) + ) / static_cast(3.0); + } + inline fp64x2 hypot(fp64x2 x, fp64x2 y) { + return sqrt( + Float64x2::Dekker_Sqr(x) + Float64x2::Dekker_Sqr(y) + ); + } + // inline fp64x2 hypot(fp64x2 x, fp64x2 y, fp64x2 z) { + // return sqrt( + // Float64x2::Dekker_Sqr(x) + Float64x2::Dekker_Sqr(y) + Float64x2::Dekker_Sqr(z) + // ); + // } /* Tests */ + inline bool signbit(fp64x2 x) { - return (x < static_cast(0.0)) ? true : false; + return (x.hi < 0.0) ? true : false; } /** Returns true if both x.hi and x.lo are finite */ inline bool isfinite(fp64x2 x) { @@ -461,6 +517,7 @@ typedef Float64x2 fp64x2; } /* Comparison */ + inline bool isgreater(fp64x2 x, fp64x2 y) { return (x > y); } @@ -478,25 +535,26 @@ typedef Float64x2 fp64x2; } /* Rounding */ - inline fp64x2 trunc(fp64x2 x) { - fp64 frac_hi = x.hi - trunc(x.hi); - fp64 frac_lo = x.lo - trunc(x.lo); - fp64x2 int_hi = trunc(x.hi); - fp64x2 int_lo = trunc(x.lo); - // Sum in increasing order - fp64x2 trunc_all = static_cast(0.0); - trunc_all += ( + + inline fp64x2 trunc(fp64x2 x) { + fp64x2 int_hi = trunc(x.hi); + fp64x2 int_lo = trunc(x.lo); + fp64 frac_hi = x.hi - int_hi.hi; + fp64 frac_lo = x.lo - int_lo.hi; + // Sum in increasing order + fp64x2 trunc_all = static_cast(0.0); + trunc_all += ( (fp64x2)frac_hi + (fp64x2)frac_lo >= static_cast(1.0) ) ? static_cast(1.0) : static_cast(0.0); - trunc_all += int_lo; - trunc_all += int_hi; - return trunc_all; - } + trunc_all += int_lo; + trunc_all += int_hi; + return trunc_all; + } inline fp64x2 floor(fp64x2 x) { fp64x2 int_part = trunc(x); return ( x < static_cast(0.0) && int_part != x - ) ? int_part : int_part - static_cast(1.0); + ) ? int_part - static_cast(1.0) : int_part; } inline fp64x2 ceil(fp64x2 x) { fp64x2 int_part = trunc(x); @@ -545,9 +603,35 @@ typedef Float64x2 fp64x2; } /* Integer and Remainder */ + + inline fp64x2 modf(fp64x2 x, fp64x2* int_part) { + fp64x2 trunc_part = trunc(x); + if (int_part != nullptr) { + *int_part = trunc_part; + } + return x - trunc_part; + } inline fp64x2 nearbyint(fp64x2 x) { return rint(x); } + + /* Float Exponents */ + + inline fp64x2 ldexp(fp64x2 x, int exp) { + x.hi = ldexp(x.hi, exp); + x.lo = isfinite(x.hi) ? ldexp(x.lo, exp) : x.hi; + return x; + } + inline fp64x2 scalbn(fp64x2 x, int exp) { + x.hi = scalbn(x.hi, exp); + x.lo = isfinite(x.hi) ? scalbn(x.lo, exp) : x.hi; + return x; + } + inline fp64x2 scalbln(fp64x2 x, long exp) { + x.hi = scalbln(x.hi, exp); + x.lo = isfinite(x.hi) ? scalbln(x.lo, exp) : x.hi; + return x; + } /* Math overloads (Casts to other types) */ @@ -558,9 +642,9 @@ typedef Float64x2 fp64x2; // inline fp64x2 fdim(fp64x2 x, fp64x2 y) { return (fp64x2)fdim((fp64x2_Math)x, (fp64x2_Math)y); } // inline fp64x2 fma(fp64x2 x, fp64x2 y, fp64x2 z) { return (fp64x2)fma((fp64x2_Math)x, (fp64x2_Math)y, (fp64x2_Math)z); } // inline fp64x2 copysign(fp64x2 x, fp64x2 y) { return (fp64x2)copysign((fp64x2_Math)x, (fp64x2_Math)y); } - inline fp64x2 sqrt(fp64x2 x) { return (fp64x2)sqrt((fp64x2_Math)x); } - inline fp64x2 cbrt(fp64x2 x) { return (fp64x2)cbrt((fp64x2_Math)x); } - inline fp64x2 hypot(fp64x2 x, fp64x2 y) { return (fp64x2)hypot((fp64x2_Math)x, (fp64x2_Math)y); } + // inline fp64x2 sqrt(fp64x2 x) { return (fp64x2)sqrt((fp64x2_Math)x); } + // inline fp64x2 cbrt(fp64x2 x) { return (fp64x2)cbrt((fp64x2_Math)x); } + // inline fp64x2 hypot(fp64x2 x, fp64x2 y) { return (fp64x2)hypot((fp64x2_Math)x, (fp64x2_Math)y); } /* Trigonometry */ inline fp64x2 sin (fp64x2 x) { return (fp64x2) sin ((fp64x2_Math)x); } inline fp64x2 cos (fp64x2 x) { return (fp64x2) cos ((fp64x2_Math)x); } @@ -601,12 +685,12 @@ typedef Float64x2 fp64x2; // inline long long llround(fp64x2 x) { return llround((fp64x2_Math)x); } /* Integer and Remainder */ inline fp64x2 fmod(fp64x2 x, fp64x2 y) { return (fp64x2)fmod((fp64x2_Math)x, (fp64x2_Math)y); } - inline fp64x2 modf(fp64x2 x, fp64x2* y) { - fp64x2_Math y_temp; - fp64x2 result = modf((fp64x2_Math)x, &y_temp); - *y = (fp64x2)y_temp; - return result; - } + // inline fp64x2 modf(fp64x2 x, fp64x2* y) { + // fp64x2_Math y_temp; + // fp64x2 result = modf((fp64x2_Math)x, &y_temp); + // *y = (fp64x2)y_temp; + // return result; + // } // inline fp64x2 nearbyint(fp64x2 x) { return (fp64x2)nearbyint((fp64x2_Math)x); } // Incorrect Function // inline fp64x2 nextafter(fp64x2 x, fp64x2 y) { return (fp64x2)nextafter((fp64x2_Math)x, (fp64x2_Math)y); } inline fp64x2 remainder(fp64x2 x, fp64x2 y) { return (fp64x2)remainder((fp64x2_Math)x, (fp64x2_Math)y); } @@ -614,9 +698,9 @@ typedef Float64x2 fp64x2; /* Float Exponents */ inline int ilogb(fp64x2 x) { return ilogb((fp64x2_Math)x); } inline fp64x2 frexp (fp64x2 x, int* exp) { return (fp64x2)frexp ((fp64x2_Math)x, exp); } - inline fp64x2 ldexp (fp64x2 x, int exp) { return (fp64x2)ldexp ((fp64x2_Math)x, exp); } - inline fp64x2 scalbn (fp64x2 x, int exp) { return (fp64x2)scalbn ((fp64x2_Math)x, exp); } - inline fp64x2 scalbln(fp64x2 x, long exp) { return (fp64x2)scalbln((fp64x2_Math)x, exp); } + // inline fp64x2 ldexp (fp64x2 x, int exp) { return (fp64x2)ldexp ((fp64x2_Math)x, exp); } + // inline fp64x2 scalbn (fp64x2 x, int exp) { return (fp64x2)scalbn ((fp64x2_Math)x, exp); } + // inline fp64x2 scalbln(fp64x2 x, long exp) { return (fp64x2)scalbln((fp64x2_Math)x, exp); } /* Tests */ // inline bool signbit(fp64x2 x) { return (signbit((fp64x2_Math)x) != 0) ? true : false; } // inline bool isfinite(fp64x2 x) { return (isfinite((fp64x2_Math)x) != 0) ? true : false; } diff --git a/src/floats/double_Float64_AVX.h b/src/floats/double_Float64_AVX.h index 5205e3c..bca7c03 100644 --- a/src/floats/double_Float64_AVX.h +++ b/src/floats/double_Float64_AVX.h @@ -1,6 +1,6 @@ /* -** Author: zerico2005 (2023 - 2024) -** Project: ABS-Fractal-Explorer +** Author: zerico2005 (2024) +** Project: LIB-Dekker-Float ** License: MIT License ** A copy of the MIT License should be included with ** this project. If not, see https://opensource.org/license/MIT @@ -21,8 +21,6 @@ */ #include -#include -#include #ifndef __AVX__ #error "__AVX__ is not enabled in your compiler. Try -mavx" @@ -48,11 +46,30 @@ inline __m256d _mm256_fabs_pd(__m256d x) { #endif #ifndef _mm256_trunc_pd +/** + * @brief _mm256_trunc_pd replacement function. + */ inline __m256d _mm256_trunc_pd(__m256d x) { return _mm256_round_pd(x, _MM_FROUND_TO_ZERO |_MM_FROUND_NO_EXC); } #endif +#ifndef _mm256_cbrt_pd +#include +/** + * @brief _mm256_cbrt_pd replacement function. May be slow/inefficient. + */ +inline __m256d _mm256_cbrt_pd(__m256d x) { + double val[4]; + _mm256_storeu_pd(val, x); + val[0] = cbrt(val[0]); + val[1] = cbrt(val[1]); + val[2] = cbrt(val[2]); + val[3] = cbrt(val[3]); + return _mm256_loadu_pd(val); +} +#endif + //------------------------------------------------------------------------------ // __m256dx2 struct //------------------------------------------------------------------------------ @@ -65,14 +82,10 @@ typedef struct __m256dx2 { __m256d lo; } __m256dx2; - //------------------------------------------------------------------------------ // __m256dx2 basic arithmetic //------------------------------------------------------------------------------ -__m256dx2 _mm256x2_max_pdx2(__m256dx2 x, __m256dx2 y); -__m256dx2 _mm256x2_min_pdx2(__m256dx2 x, __m256dx2 y); - inline __m256dx2 _mm256x2_add_pdx2(__m256dx2 x, __m256dx2 y) { // __m256d r_hi = _mm256_add_pd(x.hi, y.hi); // /* if (fabs(x.hi) < fabs(y.hi)) */ { @@ -221,6 +234,336 @@ inline __m256dx2 _mm256x2_div_pdx2(__m256dx2 x, __m256dx2 y) { return r; } +/** + * @brief returns 0 on division by 0 + */ +inline __m256dx2 _mm256x2_div_zero_pdx2(__m256dx2 x, __m256dx2 y) { + __m256dx2 u; + u.hi = _mm256_div_pd(x.hi, y.hi); + __m256dx2 t = _mm256x2_dekker_mul12_pd(u.hi, y.hi); + __m256d l = _mm256_div_pd(_mm256_sub_pd( + _mm256_add_pd(_mm256_sub_pd(_mm256_sub_pd(x.hi, t.hi), t.lo), x.lo), + _mm256_mul_pd(u.hi, y.lo) + ), y.hi); + + __m256dx2 r; + r.hi = _mm256_add_pd(u.hi, l); + r.lo = _mm256_add_pd(_mm256_sub_pd(u.hi, r.hi), l); + + __m256d cmp_zero = _mm256_cmp_pd(y.hi, _mm256_setzero_pd(), _CMP_EQ_OS); + r.hi = _mm256_andnot_pd(r.hi, cmp_zero); + r.lo = _mm256_andnot_pd(r.lo, cmp_zero); + return r; +} + +inline __m256dx2 _mm256x2_dekker_square12_pd(__m256d x) { + __m256dx2 a = _mm256x2_dekker_split_pd(x); + __m256d p = _mm256_mul_pd(a.hi, a.hi); + __m256d q = _mm256_mul_pd( + _mm256_set1_pd(2.0), _mm256_mul_pd(a.hi, a.lo) + ); + + __m256dx2 r; + r.hi = _mm256_add_pd(p, q); + r.lo = _mm256_add_pd( + _mm256_add_pd(_mm256_sub_pd(p, r.hi), q), + _mm256_mul_pd(a.lo, a.lo) + ); + return r; +} + +inline __m256dx2 _mm256x2_square_pdx2(__m256dx2 x) { + __m256dx2 t = _mm256x2_dekker_square12_pd(x.hi); + __m256d c = _mm256_add_pd( + _mm256_mul_pd(_mm256_set1_pd(2.0), _mm256_mul_pd(x.hi, x.lo)), t.lo + ); + + __m256dx2 r; + r.hi = _mm256_add_pd(t.hi, c); + r.lo = _mm256_add_pd(_mm256_sub_pd(t.hi, r.hi), c); + return r; +} + +//------------------------------------------------------------------------------ +// __m256dx2 optimized arithmetic +//------------------------------------------------------------------------------ + +inline __m256dx2 _mm256x2_add_pdx2_pd(__m256dx2 x, __m256d y) { + __m256d r_hi = _mm256_add_pd(x.hi, y); + + __m256d rx_lo = _mm256_add_pd( + _mm256_add_pd(_mm256_sub_pd(x.hi, r_hi), y), x.lo + ); + __m256d ry_lo = _mm256_add_pd( + _mm256_add_pd(_mm256_sub_pd(y, r_hi), x.hi), + x.lo); + + const __m256d cmp_result = _mm256_cmp_pd( + _mm256_fabs_pd(x.hi), _mm256_fabs_pd(y), + _CMP_LE_OQ); + __m256d r_lo = _mm256_blendv_pd(rx_lo, ry_lo, cmp_result); + + __m256dx2 c; + c.hi = _mm256_add_pd(r_hi, r_lo); + c.lo = _mm256_add_pd(_mm256_sub_pd(r_hi, c.hi), r_lo); + return c; +} + +inline __m256dx2 _mm256x2_add_pd_pdx2(__m256d x, __m256dx2 y) { + __m256d r_hi = _mm256_add_pd(x, y.hi); + + __m256d rx_lo = _mm256_add_pd( + _mm256_add_pd(_mm256_sub_pd(x, r_hi), y.hi), + y.lo); + __m256d ry_lo = _mm256_add_pd( + _mm256_add_pd(_mm256_sub_pd(y.hi, r_hi), x), y.lo + ); + + const __m256d cmp_result = _mm256_cmp_pd( + _mm256_fabs_pd(x), _mm256_fabs_pd(y.hi), + _CMP_LE_OQ); + __m256d r_lo = _mm256_blendv_pd(rx_lo, ry_lo, cmp_result); + + __m256dx2 c; + c.hi = _mm256_add_pd(r_hi, r_lo); + c.lo = _mm256_add_pd(_mm256_sub_pd(r_hi, c.hi), r_lo); + return c; +} + +/** + * @brief Adds two __m256d values with the result stored as a __m256dx2 + */ +inline __m256dx2 _mm256x2_add_pd_pd(__m256d x, __m256d y) { + __m256d r_hi = _mm256_add_pd(x, y); + + __m256d rx_lo = _mm256_add_pd(_mm256_sub_pd(x, r_hi), y); + __m256d ry_lo = _mm256_add_pd(_mm256_sub_pd(y, r_hi), x); + + const __m256d cmp_result = _mm256_cmp_pd( + _mm256_fabs_pd(x), _mm256_fabs_pd(y), + _CMP_LE_OQ); + __m256d r_lo = _mm256_blendv_pd(rx_lo, ry_lo, cmp_result); + + __m256dx2 c; + c.hi = _mm256_add_pd(r_hi, r_lo); + c.lo = _mm256_add_pd(_mm256_sub_pd(r_hi, c.hi), r_lo); + return c; +} + +inline __m256dx2 _mm256x2_sub_pdx2_pd(__m256dx2 x, __m256d y) { + __m256d r_hi = _mm256_sub_pd(x.hi, y); + + __m256d rx_lo = _mm256_add_pd( + _mm256_sub_pd(_mm256_sub_pd(x.hi, r_hi), y), x.lo + ); + + y = _mm256_mul_pd(y, _mm256_set1_pd(-1.0)); + __m256d ry_lo = _mm256_add_pd( + _mm256_add_pd(_mm256_sub_pd(/* negative */ y, r_hi), x.hi), x.lo + ); + + const __m256d cmp_result = _mm256_cmp_pd( + _mm256_fabs_pd(x.hi), _mm256_fabs_pd(y), + _CMP_LE_OQ); + __m256d r_lo = _mm256_blendv_pd(rx_lo, ry_lo, cmp_result); + + __m256dx2 c; + c.hi = _mm256_add_pd(r_hi, r_lo); + c.lo = _mm256_add_pd(_mm256_sub_pd(r_hi, c.hi), r_lo); + return c; +} + +inline __m256dx2 _mm256x2_sub_pd_pdx2(__m256d x, __m256dx2 y) { + __m256d r_hi = _mm256_sub_pd(x, y.hi); + + __m256d rx_lo = _mm256_sub_pd( + _mm256_sub_pd(_mm256_sub_pd(x, r_hi), y.hi), + y.lo); + + y.hi = _mm256_mul_pd(y.hi, _mm256_set1_pd(-1.0)); + __m256d ry_lo = _mm256_sub_pd( + _mm256_add_pd(_mm256_sub_pd(/* negative */ y.hi, r_hi), x), y.lo + ); + + const __m256d cmp_result = _mm256_cmp_pd( + _mm256_fabs_pd(x), _mm256_fabs_pd(y.hi), + _CMP_LE_OQ); + __m256d r_lo = _mm256_blendv_pd(rx_lo, ry_lo, cmp_result); + + __m256dx2 c; + c.hi = _mm256_add_pd(r_hi, r_lo); + c.lo = _mm256_add_pd(_mm256_sub_pd(r_hi, c.hi), r_lo); + return c; +} + +/** + * @brief Subtracts two __m256d values with the result stored as a __m256dx2 + */ +inline __m256dx2 _mm256x2_sub_pd_pd(__m256d x, __m256d y) { + __m256d r_hi = _mm256_sub_pd(x, y); + + __m256d rx_lo = _mm256_sub_pd(_mm256_sub_pd(x, r_hi), y); + + y = _mm256_mul_pd(y, _mm256_set1_pd(-1.0)); + __m256d ry_lo = _mm256_add_pd(_mm256_sub_pd(/* negative */ y, r_hi), x); + + const __m256d cmp_result = _mm256_cmp_pd( + _mm256_fabs_pd(x), _mm256_fabs_pd(y), + _CMP_LE_OQ); + __m256d r_lo = _mm256_blendv_pd(rx_lo, ry_lo, cmp_result); + + __m256dx2 c; + c.hi = _mm256_add_pd(r_hi, r_lo); + c.lo = _mm256_add_pd(_mm256_sub_pd(r_hi, c.hi), r_lo); + return c; +} + +inline __m256dx2 _mm256x2_mul_pdx2_pd(__m256dx2 x, __m256d y) { + __m256dx2 t = _mm256x2_dekker_mul12_pd(x.hi, y); + __m256d c = _mm256_add_pd(_mm256_mul_pd(x.lo, y), t.lo); + + __m256dx2 r; + r.hi = _mm256_add_pd(t.hi, c); + r.lo = _mm256_add_pd(_mm256_sub_pd(t.hi, r.hi), c); + return r; +} + +inline __m256dx2 _mm256x2_mul_pd_pdx2(__m256d x, __m256dx2 y) { + __m256dx2 t = _mm256x2_dekker_mul12_pd(x, y.hi); + __m256d c = _mm256_add_pd(_mm256_mul_pd(x, y.lo), t.lo); + + __m256dx2 r; + r.hi = _mm256_add_pd(t.hi, c); + r.lo = _mm256_add_pd(_mm256_sub_pd(t.hi, r.hi), c); + return r; +} + +/** + * @brief Multiplies two __m256d values with the result stored as a __m256dx2 + */ +inline __m256dx2 _mm256x2_mul_pd_pd(__m256d x, __m256d y) { + return _mm256x2_dekker_mul12_pd(x, y); +} + +inline __m256dx2 _mm256x2_div_pdx2_pd(__m256dx2 x, __m256d y) { + __m256dx2 u; + u.hi = _mm256_div_pd(x.hi, y); + __m256dx2 t = _mm256x2_dekker_mul12_pd(u.hi, y); + __m256d l = _mm256_div_pd( + _mm256_add_pd(_mm256_sub_pd(_mm256_sub_pd(x.hi, t.hi), t.lo), x.lo), y + ); + + __m256dx2 r; + r.hi = _mm256_add_pd(u.hi, l); + r.lo = _mm256_add_pd(_mm256_sub_pd(u.hi, r.hi), l); + return r; +} + +inline __m256dx2 _mm256x2_div_pd_pdx2(__m256d x, __m256dx2 y) { + __m256dx2 u; + u.hi = _mm256_div_pd(x, y.hi); + __m256dx2 t = _mm256x2_dekker_mul12_pd(u.hi, y.hi); + __m256d l = _mm256_div_pd(_mm256_sub_pd( + _mm256_sub_pd(_mm256_sub_pd(x, t.hi), t.lo), + _mm256_mul_pd(u.hi, y.lo) + ), y.hi); + + __m256dx2 r; + r.hi = _mm256_add_pd(u.hi, l); + r.lo = _mm256_add_pd(_mm256_sub_pd(u.hi, r.hi), l); + return r; +} + +/** + * @brief Divides two __m256d values with the result stored as a __m256dx2 + */ +inline __m256dx2 _mm256x2_div_pd_pd(__m256d x, __m256d y) { + __m256dx2 u; + u.hi = _mm256_div_pd(x, y); + __m256dx2 t = _mm256x2_dekker_mul12_pd(u.hi, y); + __m256d l = _mm256_div_pd( + _mm256_sub_pd(_mm256_sub_pd(x, t.hi), t.lo), y + ); + + __m256dx2 r; + r.hi = _mm256_add_pd(u.hi, l); + r.lo = _mm256_add_pd(_mm256_sub_pd(u.hi, r.hi), l); + return r; +} + +/** + * @brief Returns 0 on division by 0 + */ +inline __m256dx2 _mm256x2_div_zero_pdx2_pd(__m256dx2 x, __m256d y) { + __m256dx2 u; + u.hi = _mm256_div_pd(x.hi, y); + __m256dx2 t = _mm256x2_dekker_mul12_pd(u.hi, y); + __m256d l = _mm256_div_pd( + _mm256_add_pd(_mm256_sub_pd(_mm256_sub_pd(x.hi, t.hi), t.lo), x.lo), y + ); + + __m256dx2 r; + r.hi = _mm256_add_pd(u.hi, l); + r.lo = _mm256_add_pd(_mm256_sub_pd(u.hi, r.hi), l); + + __m256d cmp_zero = _mm256_cmp_pd(y, _mm256_setzero_pd(), _CMP_EQ_OS); + r.hi = _mm256_andnot_pd(r.hi, cmp_zero); + r.lo = _mm256_andnot_pd(r.lo, cmp_zero); + return r; +} + +/** + * @brief Returns 0 on division by 0 + */ +inline __m256dx2 _mm256x2_div_zero_pd_pdx2(__m256d x, __m256dx2 y) { + __m256dx2 u; + u.hi = _mm256_div_pd(x, y.hi); + __m256dx2 t = _mm256x2_dekker_mul12_pd(u.hi, y.hi); + __m256d l = _mm256_div_pd(_mm256_sub_pd( + _mm256_sub_pd(_mm256_sub_pd(x, t.hi), t.lo), + _mm256_mul_pd(u.hi, y.lo) + ), y.hi); + + __m256dx2 r; + r.hi = _mm256_add_pd(u.hi, l); + r.lo = _mm256_add_pd(_mm256_sub_pd(u.hi, r.hi), l); + + __m256d cmp_zero = _mm256_cmp_pd(y.hi, _mm256_setzero_pd(), _CMP_EQ_OS); + r.hi = _mm256_andnot_pd(r.hi, cmp_zero); + r.lo = _mm256_andnot_pd(r.lo, cmp_zero); + return r; +} + + +/** + * @brief Returns 0 on division by 0. Divides two __m256d values with the + * result stored as a __m256dx2 + */ +inline __m256dx2 _mm256x2_div_zero_pd_pd(__m256d x, __m256d y) { + __m256dx2 u; + u.hi = _mm256_div_pd(x, y); + __m256dx2 t = _mm256x2_dekker_mul12_pd(u.hi, y); + __m256d l = _mm256_div_pd( + _mm256_sub_pd(_mm256_sub_pd(x, t.hi), t.lo), y + ); + + __m256dx2 r; + r.hi = _mm256_add_pd(u.hi, l); + r.lo = _mm256_add_pd(_mm256_sub_pd(u.hi, r.hi), l); + + __m256d cmp_zero = _mm256_cmp_pd(y, _mm256_setzero_pd(), _CMP_EQ_OS); + r.hi = _mm256_andnot_pd(r.hi, cmp_zero); + r.lo = _mm256_andnot_pd(r.lo, cmp_zero); + return r; +} + +/** + * @brief Squares a __m256d value with the result stored as a __m256dx2 + */ +inline __m256dx2 _mm256x2_square_pd(__m256d x) { + return _mm256x2_dekker_square12_pd(x); +} + //------------------------------------------------------------------------------ // __m256dx2 bitwise operations //------------------------------------------------------------------------------ @@ -595,34 +938,6 @@ inline __m256d _mm256_cmpge_pdx2(__m256dx2 x, __m256dx2 y) { return _mm256_blendv_pd(cmp_hi, cmp_lo, cmp_eq); } -// /** -// * @brief _CMP_LT_OQ -// */ -// inline __m256d _mm256_cmplt_pdx2(__m256dx2 x, __m256dx2 y) { -// return _mm256_cmp_pd(x.hi, y.hi, _CMP_LT_OQ); -// } - -// /** -// * @brief _CMP_LE_OQ -// */ -// inline __m256d _mm256_cmple_pdx2(__m256dx2 x, __m256dx2 y) { -// return _mm256_cmp_pd(x.hi, y.hi, _CMP_LE_OQ); -// } - -// /** -// * @brief _CMP_GT_OQ -// */ -// inline __m256d _mm256_cmpgt_pdx2(__m256dx2 x, __m256dx2 y) { -// return _mm256_cmp_pd(x.hi, y.hi, _CMP_GT_OQ); -// } - -// /** -// * @brief _CMP_GE_OQ -// */ -// inline __m256d _mm256_cmpge_pdx2(__m256dx2 x, __m256dx2 y) { -// return _mm256_cmp_pd(x.hi, y.hi, _CMP_GE_OQ); -// } - //------------------------------------------------------------------------------ // _mm256_cmp_pdx2 AVX compare function //------------------------------------------------------------------------------ @@ -673,6 +988,52 @@ inline __m256dx2 _mm256x2_min_pdx2(__m256dx2 x, __m256dx2 y) { return ret; } +//------------------------------------------------------------------------------ +// __m256dx2 rounding functions +//------------------------------------------------------------------------------ + +inline __m256dx2 _mm256x2_trunc_pdx2(__m256dx2 x) { + __m256d int_hi = _mm256_trunc_pd(x.hi); + __m256d int_lo = _mm256_trunc_pd(x.lo); + __m256d frac_hi = _mm256_sub_pd(x.hi, int_hi); + __m256d frac_lo = _mm256_sub_pd(x.lo, int_lo); + + __m256d frac_ge_1 = _mm256_cmp_pd(_mm256_add_pd(frac_hi, frac_lo), _mm256_set1_pd(1.0), _CMP_GE_OQ); + + __m256dx2 trunc_all = _mm256x2_add_pdx2_pd(_mm256x2_add_pd_pd( + _mm256_blendv_pd(_mm256_setzero_pd(), _mm256_set1_pd(1.0), frac_ge_1), + int_lo + ), int_hi); + return trunc_all; +} + +inline __m256dx2 _mm256x2_floor_pdx2(__m256dx2 x) { + __m256dx2 int_part = _mm256x2_trunc_pdx2(x); + + __m256d cmp_floor; + cmp_floor = _mm256_and_pd( + _mm256_cmplt_pdx2(x, _mm256x2_setzero_pdx2()), + _mm256_cmpneq_pdx2(x, int_part) + ); + __m256d floor_part = _mm256_blendv_pd( + _mm256_setzero_pd(), _mm256_set1_pd(1.0), cmp_floor + ); + return _mm256x2_sub_pdx2_pd(int_part, floor_part); +} + +inline __m256dx2 _mm256x2_ceil_pdx2(__m256dx2 x) { + __m256dx2 int_part = _mm256x2_trunc_pdx2(x); + + __m256d cmp_ceil = _mm256_and_pd( + _mm256_cmpgt_pdx2(x, _mm256x2_setzero_pdx2()), + _mm256_cmpneq_pdx2(x, int_part) + ); + __m256d ceil_part = _mm256_blendv_pd( + _mm256_setzero_pd(), _mm256_set1_pd(1.0), cmp_ceil + ); + return _mm256x2_add_pdx2_pd(int_part, ceil_part); +} + //------------------------------------------------------------------------------ // __m256dx2 math.h functions //------------------------------------------------------------------------------ @@ -697,6 +1058,51 @@ inline __m256dx2 _mm256x2_fdim_pdx2(__m256dx2 x, __m256dx2 y) { return ret; } +inline __m256dx2 _mm256x2_copysign_pdx2(__m256dx2 x, __m256dx2 y) { + __m256d negate_mask = _mm256_xor_pd( + _mm256_cmp_pd(x.hi, _mm256_setzero_pd(), _CMP_LT_OQ), + _mm256_cmp_pd(y.hi, _mm256_setzero_pd(), _CMP_LT_OQ) + ); + __m256d negate_mul = _mm256_blendv_pd( + _mm256_set1_pd(1.0), _mm256_set1_pd(-1.0), negate_mask + ); + x.hi = _mm256_mul_pd(x.hi, negate_mul); + x.lo = _mm256_mul_pd(x.lo, negate_mul); + return x; +} + +inline __m256dx2 _mm256x2_sqrt_pdx2(__m256dx2 x) { + __m256d guess = _mm256_sqrt_pd(x.hi); + return _mm256x2_mul_pdx2_pd( + _mm256x2_add_pd_pdx2(guess, _mm256x2_div_zero_pdx2_pd(x, guess)), + _mm256_set1_pd(0.5) + ); +} + +inline __m256dx2 _mm256x2_cbrt_pdx2(__m256dx2 x) { + __m256d guess = _mm256_cbrt_pd(x.hi); + return _mm256x2_div_pdx2_pd( + _mm256x2_add_pdx2( + _mm256x2_mul_pd_pd( + _mm256_set1_pd(2.0), guess + ), + _mm256x2_div_zero_pdx2(x, _mm256x2_square_pd(guess)) + ), + _mm256_set1_pd(3.0) + ); +} + +/** + * @brief returns the fraction part of a __m256dx2 value. int_part may be NULL + */ +inline __m256dx2 _mm256x2_modf_pdx2(__m256dx2 x, __m256dx2* int_part) { + __m256dx2 trunc_part = _mm256x2_trunc_pdx2(x); + if (int_part != NULL) { + *int_part = trunc_part; + } + return _mm256x2_sub_pdx2(x, trunc_part); +} + #ifdef __cplusplus } #endif diff --git a/src/floats/double_Float64_SSE2.h b/src/floats/double_Float64_SSE2.h deleted file mode 100644 index b1091c3..0000000 --- a/src/floats/double_Float64_SSE2.h +++ /dev/null @@ -1,708 +0,0 @@ -/* -** Author: zerico2005 (2023 - 2024) -** Project: ABS-Fractal-Explorer -** License: MIT License -** A copy of the MIT License should be included with -** this project. If not, see https://opensource.org/license/MIT -*/ -#ifndef DOUBLE_FLOAT64_SSE_H -#define DOUBLE_FLOAT64_SSE_H - -/** - * @brief Double-Float64 SSE2 Dekker Float implementation. - * Source: Creel "Double it Like Dekker" on YouTube. - * - * @note Requires SSE2 or later. - * SSE3 allows for some optimizations if enabled. - * SSE4.1 allows for faster floor/ceil/round functions - */ - -#include -#include -#include - -#ifndef __SSE2__ - #error "SSE2 is not enabled in your compiler. Try -msse2 -#endif - -#include -#ifdef __SSE3__ - #include -#endif - -/** - * @brief Holds two Double-Float64 dekker floats - */ -typedef struct __m128dx2 { - __m128d hi; - __m128d lo; -} __m128dx2; - -typedef float fp32; -typedef double fp64; - -#ifndef _mm_fabs_pd -/** - * @remarks _mm_andnot_pd cannot be used because 0x8000000000000000 gets - * converted from -0.0 to 0.0 on -Ofast - */ -inline __m128d _mm_fabs_pd(const __m128d x) { - return - _mm_and_pd( - x, - _mm_castsi128_pd(_mm_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF)) - ); -} -#endif - -inline __m128dx2 _mm_add_pdx2(__m128dx2 x, __m128dx2 y) { - __m128d r_hi = _mm_add_pd(x.hi, y.hi); - __m128d r_lo = _mm_setzero_pd(); - __m128d cmp_result = __m128d_mm_cmpgt_pd(_mm_fabs_pd(x.hi), _mm_fabs_pd(y.hi)); - x.hi = _mm_shuffle_pd(x.hi, y.lo); - y = - r_lo = _mm_add_pd(_mm_sub_pd(x.hi, r_hi)); - r_lo = x.hi - r_hi + y.hi + y.lo + x.lo; - r_lo = y.hi - r_hi + x.hi + x.lo + y.lo; - - __m128dx2 c; - c.hi = r_hi + r_lo; - c.lo = r_hi - c.hi + r_lo; - return c; -} - -class Float64x2 { -public: - fp64 hi; - fp64 lo; - -private: - - /* Arithmetic */ - - inline Float64x2 Dekker_Add( - const Float64x2& x, const Float64x2& y - ) const { - fp64 r_hi = x.hi + y.hi; - fp64 r_lo = 0.0; - if (fabs(x.hi) > fabs(y.hi)) { - r_lo = x.hi - r_hi + y.hi + y.lo + x.lo; - } else { - r_lo = y.hi - r_hi + x.hi + x.lo + y.lo; - } - - Float64x2 c; - c.hi = r_hi + r_lo; - c.lo = r_hi - c.hi + r_lo; - return c; - } - - inline Float64x2 Dekker_Sub( - const Float64x2& x, const Float64x2& y - ) const { - fp64 r_hi = x.hi - y.hi; - fp64 r_lo = 0.0; - if (fabs(x.hi) > fabs(y.hi)) { - r_lo = x.hi - r_hi - y.hi - y.lo + x.lo; - } else { - r_lo = -y.hi - r_hi + x.hi + x.lo - y.lo; - } - - Float64x2 c; - c.hi = r_hi + r_lo; - c.lo = r_hi - c.hi + r_lo; - return c; - } - - static constexpr fp64 Dekker_Scale = 134217729.0; // (2^ceil(53 / 2) + 1) - - inline Float64x2 Dekker_Split(const fp64& x) const { - fp64 p = x * Dekker_Scale; - Float64x2 r; - r.hi = x - p + p; - r.lo = x - r.hi; - return r; - } - - // static constexpr uint64_t Dekker_Split_Mask = ~((uint64_t)0x3FFFFFF); - - // inline Float64x2 Dekker_Split(const fp64& x) const { - // Float64x2 r; - // uint64_t temp = (*(uint64_t*)((void*)&x)) & Dekker_Split_Mask; - // r.hi = (*(fp64*)((void*)&temp)); - // r.lo = x - r.hi; - // return r; - // } - - inline Float64x2 Dekker_Mul12( - const fp64& x, const fp64& y - ) const { - Float64x2 a = Dekker_Split(x); - Float64x2 b = Dekker_Split(y); - fp64 p = a.hi * b.hi; - fp64 q = a.hi * b.lo + a.lo * b.hi; - - Float64x2 r; - r.hi = p + q; - r.lo = p - r.hi + q + a.lo * b.lo; - return r; - } - - inline Float64x2 Dekker_Mul( - const Float64x2& x, const Float64x2& y - ) const { - Float64x2 t = Dekker_Mul12(x.hi, y.hi); - fp64 c = x.hi * y.lo + x.lo * y.hi + t.lo; - - Float64x2 r; - r.hi = t.hi + c; - r.lo = t.hi - r.hi + c; - return r; - } - - inline Float64x2 Dekker_Div( - const Float64x2& x, const Float64x2& y - ) const { - Float64x2 u; - u.hi = x.hi / y.hi; - Float64x2 t = Dekker_Mul12(u.hi, y.hi); - fp64 l = (x.hi - t.hi - t.lo + x.lo - u.hi * y.lo) / y.hi; - - Float64x2 r; - r.hi = u.hi + l; - r.lo = u.hi - r.hi + l; - return r; - } - - /* Double-Single Arithmetic */ - - inline Float64x2 Dekker_Add_Float64( - const Float64x2& x, const fp64& y - ) const { - fp64 r_hi = x.hi + y; - fp64 r_lo = 0.0f; - if (fabs(x.hi) > fabs(y)) { - r_lo = x.hi - r_hi + y + x.lo; - } else { - r_lo = y - r_hi + x.hi + x.lo; - } - - Float64x2 c; - c.hi = r_hi + r_lo; - c.lo = r_hi - c.hi + r_lo; - return c; - } - - inline Float64x2 Dekker_Sub_Float64( - const Float64x2& x, const fp64& y - ) const { - fp64 r_hi = x.hi - y; - fp64 r_lo = 0.0f; - if (fabs(x.hi) > fabs(y)) { - r_lo = x.hi - r_hi - y + x.lo; - } else { - r_lo = -y - r_hi + x.hi + x.lo; - } - - Float64x2 c; - c.hi = r_hi + r_lo; - c.lo = r_hi - c.hi + r_lo; - return c; - } - - inline Float64x2 Dekker_Mul_Float64( - const Float64x2& x, const fp64& y - ) const { - Float64x2 t = Dekker_Mul12(x.hi, y); - fp64 c = x.lo * y + t.lo; - - Float64x2 r; - r.hi = t.hi + c; - r.lo = t.hi - r.hi + c; - return r; - } - - inline Float64x2 Dekker_Div_Float64( - const Float64x2& x, const fp64& y - ) const { - Float64x2 u; - u.hi = x.hi / y; - Float64x2 t = Dekker_Mul12(u.hi, y); - fp64 l = (x.hi - t.hi - t.lo + x.lo) / y; - - Float64x2 r; - r.hi = u.hi + l; - r.lo = u.hi - r.hi + l; - return r; - } - - inline Float64x2 Float64_Div_Dekker( - const fp64& x, const Float64x2& y - ) const { - Float64x2 u; - u.hi = x / y.hi; - Float64x2 t = Dekker_Mul12(u.hi, y.hi); - fp64 l = (x - t.hi - t.lo - u.hi * y.lo) / y.hi; - - Float64x2 r; - r.hi = u.hi + l; - r.lo = u.hi - r.hi + l; - return r; - } - -public: - -/* Arithmetic */ - - inline Float64x2 operator+(const Float64x2 &value) const { - return Dekker_Add(*this, value); - } - - inline Float64x2 operator-(const Float64x2 &value) const { - return Dekker_Sub(*this, value); - } - - inline Float64x2 operator*(const Float64x2 &value) const { - return Dekker_Mul(*this, value); - } - - inline Float64x2 operator/(const Float64x2 &value) const { - return Dekker_Div(*this, value); - } - - inline Float64x2 operator+(const fp64 &value) const { - return Dekker_Add_Float64(*this, value); - } - - inline Float64x2 operator-(const fp64 &value) const { - return Dekker_Sub_Float64(*this, value); - } - - inline Float64x2 operator*(const fp64 &value) const { - return Dekker_Mul_Float64(*this, value); - } - - inline Float64x2 operator/(const fp64 &value) const { - return Dekker_Div_Float64(*this, value); - } - - inline Float64x2 operator-() const { - Float64x2 value = *this; - value.hi = -value.hi; - value.lo = -value.lo; - return value; - } - -/* Increment/Decrement */ - - inline Float64x2& operator++() { - *this = Dekker_Add(*this, static_cast(1.0)); - return *this; - } - - inline Float64x2& operator--() { - *this = Dekker_Sub(*this, static_cast(1.0)); - return *this; - } - - inline Float64x2 operator++(int) { - Float64x2 temp = *this; - *this = Dekker_Add(*this, static_cast(1.0)); - return temp; - } - - inline Float64x2 operator--(int) { - Float64x2 temp = *this; - *this = Dekker_Sub(*this, static_cast(1.0)); - return temp; - } - -/* Compound Assignment */ - - inline Float64x2& operator+=(const Float64x2 &value) { - *this = Dekker_Add(*this, value); - return *this; - } - - inline Float64x2& operator-=(const Float64x2 &value) { - *this = Dekker_Sub(*this, value); - return *this; - } - - inline Float64x2& operator*=(const Float64x2 &value) { - *this = Dekker_Mul(*this, value); - return *this; - } - - inline Float64x2& operator/=(const Float64x2 &value) { - *this = Dekker_Div(*this, value); - return *this; - } - -/* Double-Single Compound Assignment */ - - inline Float64x2& operator+=(const fp64 &value) { - *this = Dekker_Add_Float64(*this, value); - return *this; - } - - inline Float64x2& operator-=(const fp64 &value) { - *this = Dekker_Sub_Float64(*this, value); - return *this; - } - - inline Float64x2& operator*=(const fp64 &value) { - *this = Dekker_Mul_Float64(*this, value); - return *this; - } - - inline Float64x2& operator/=(const fp64 &value) { - *this = Dekker_Div_Float64(*this, value); - return *this; - } - -/* Comparison */ - - inline bool operator==(const Float64x2 &value) const { - return ( - this->hi == value.hi && - this->lo == value.lo - ) ? true : false; - } - inline bool operator!=(const Float64x2 &value) const { - return ( - this->hi != value.hi || - this->lo != value.lo - ) ? true : false; - } - - inline bool operator<(const Float64x2 &value) const { - if (this->hi == value.hi) { - return (this->lo < value.lo); - } - return (this->hi < value.hi); - } - - inline bool operator<=(const Float64x2 &value) const { - if (this->hi == value.hi) { - return (this->lo <= value.lo); - } - return (this->hi <= value.hi); - } - - inline bool operator>(const Float64x2 &value) const { - if (this->hi == value.hi) { - return (this->lo > value.lo); - } - return (this->hi > value.hi); - } - - inline bool operator>=(const Float64x2 &value) const { - if (this->hi == value.hi) { - return (this->lo >= value.lo); - } - return (this->hi >= value.hi); - } - -/* Constructors */ - - inline Float64x2() { - // this->hi = 0.0; - // this->lo = 0.0; - } - - inline Float64x2(const fp32& value) { - this->hi = (fp64)value; - this->lo = 0.0; - } - inline Float64x2(const fp64& value) { - this->hi = value; - this->lo = 0.0; - } - - template - inline Float64x2(const fpX& value) { - this->hi = (fp64)value; - this->lo = (fp64)(value - (fpX)this->hi); - } - -/* Casts */ - - constexpr inline operator fp32() const { - return (fp32)this->hi; - } - constexpr inline operator fp64() const { - return this->hi; - } - - template - constexpr inline operator fpX() const { - return (fpX)this->hi + (fpX)this->lo; - } - -}; - -typedef Float64x2 fp64x2; - -/* Math functions (Natively implemented) */ - - /* Arithmetic */ - inline fp64x2 fmax(fp64x2 x, fp64x2 y) { - return (x > y) ? x : y; - } - inline fp64x2 fmin(fp64x2 x, fp64x2 y) { - return (x < y) ? x : y; - } - inline fp64x2 fabs(fp64x2 x) { - return (x < static_cast(0.0)) ? -x : x; - } - inline fp64x2 fdim(fp64x2 x, fp64x2 y) { - return (x > y) ? (x - y) : static_cast(0.0); - } - inline fp64x2 fma(fp64x2 x, fp64x2 y, fp64x2 z) { - return (x * y) + z; - } - inline fp64x2 copysign(fp64x2 x, fp64x2 y) { - return ( - (x < static_cast(0.0)) != (y < (static_cast(0.0))) - ) ? -x : x; - } - - /* Tests */ - inline bool signbit(fp64x2 x) { - return (x < static_cast(0.0)) ? true : false; - } - /** Returns true if both x.hi and x.lo are finite */ - inline bool isfinite(fp64x2 x) { - return (isfinite(x.hi) && isfinite(x.lo)); - } - /** Returns true if either x.hi or x.lo are infinite */ - inline bool isinf(fp64x2 x) { - return (isinf(x.hi) || isinf(x.lo)); - } - /** Returns true if either x.hi or x.lo are nan */ - inline bool isnan(fp64x2 x) { - return (isnan(x.hi) || isnan(x.lo)); - } - /** Returns true if both x.hi and x.lo are normal */ - inline bool isnormal(fp64x2 x) { - return (isnormal(x.hi) && isnormal(x.lo)); - } - /** Returns true if either {x.hi, y.hi} or {x.lo, y.lo} are unordered */ - inline bool isunordered(fp64x2 x, fp64x2 y) { - return (isunordered(x.hi, y.hi) || isunordered(x.lo, y.lo)); - } - inline int fpclassify(fp64x2 x) { - if (isfinite(x)) { return FP_INFINITE; } - if (isnan(x)) { return FP_NAN; } - if (x == static_cast(0.0)) { return FP_ZERO; } - if (isnormal(x)) { return FP_NORMAL; } - return FP_SUBNORMAL; - } - - /* Comparison */ - inline bool isgreater(fp64x2 x, fp64x2 y) { - return (x > y); - } - inline bool isgreaterequal(fp64x2 x, fp64x2 y) { - return (x >= y); - } - inline bool isless(fp64x2 x, fp64x2 y) { - return (x < y); - } - inline bool islessequal(fp64x2 x, fp64x2 y) { - return (x <= y); - } - inline bool islessgreater(fp64x2 x, fp64x2 y) { - return (x < y) || (x > y); - } - - /* Rounding */ - inline fp64x2 trunc(fp64x2 x) { - fp64 frac_hi = x.hi - trunc(x.hi); - fp64 frac_lo = x.lo - trunc(x.lo); - fp64x2 int_hi = trunc(x.hi); - fp64x2 int_lo = trunc(x.lo); - // Sum in increasing order - fp64x2 trunc_all = static_cast(0.0); - trunc_all += ( - (fp64x2)frac_hi + (fp64x2)frac_lo >= static_cast(1.0) - ) ? static_cast(1.0) : static_cast(0.0); - trunc_all += int_lo; - trunc_all += int_hi; - return trunc_all; - } - inline fp64x2 floor(fp64x2 x) { - fp64x2 int_part = trunc(x); - return ( - x < static_cast(0.0) && int_part != x - ) ? int_part : int_part - static_cast(1.0); - } - inline fp64x2 ceil(fp64x2 x) { - fp64x2 int_part = trunc(x); - return ( - x > static_cast(0.0) && int_part != x - ) ? int_part + static_cast(1.0) : int_part; - } - inline fp64x2 round(fp64x2 x) { - fp64x2 int_part = trunc(x); - fp64x2 frac_part = x - int_part; - if (x >= static_cast(0.0)) { - if (frac_part >= static_cast(0.5)) { - return int_part + static_cast(1.0); - } - return int_part; - } - if (frac_part <= static_cast(-0.5)) { - return int_part - static_cast(1.0); - } - return int_part; - } - inline fp64x2 rint(fp64x2 x) { - switch (fegetround()) { - default: - case FE_TOWARDZERO: - return trunc(x); - case FE_DOWNWARD: - return floor(x); - case FE_UPWARD: - return ceil(x); - case FE_TONEAREST: - return round(x); - } - } - inline long lround(fp64x2 x) { - return (long)round(x); - } - inline long lrint(fp64x2 x) { - return (long)rint(x); - } - inline long long llround(fp64x2 x) { - return (long long)round(x); - } - inline long long llrint(fp64x2 x) { - return (long long)rint(x); - } - - /* Integer and Remainder */ - inline fp64x2 nearbyint(fp64x2 x) { - return rint(x); - } - -/* Math overloads (Casts to other types) */ - - /* Arithmetic */ - // inline fp64x2 fmax(fp64x2 x, fp64x2 y) { return (fp64x2)fmax((fp64x2_Math)x, (fp64x2_Math)y); } - // inline fp64x2 fmin(fp64x2 x, fp64x2 y) { return (fp64x2)fmin((fp64x2_Math)x, (fp64x2_Math)y); } - // inline fp64x2 fabs(fp64x2 x) { return (fp64x2)fabs((fp64x2_Math)x); } - // inline fp64x2 fdim(fp64x2 x, fp64x2 y) { return (fp64x2)fdim((fp64x2_Math)x, (fp64x2_Math)y); } - // inline fp64x2 fma(fp64x2 x, fp64x2 y, fp64x2 z) { return (fp64x2)fma((fp64x2_Math)x, (fp64x2_Math)y, (fp64x2_Math)z); } - // inline fp64x2 copysign(fp64x2 x, fp64x2 y) { return (fp64x2)copysign((fp64x2_Math)x, (fp64x2_Math)y); } - inline fp64x2 sqrt(fp64x2 x) { return (fp64x2)sqrt((fp64x2_Math)x); } - inline fp64x2 cbrt(fp64x2 x) { return (fp64x2)cbrt((fp64x2_Math)x); } - inline fp64x2 hypot(fp64x2 x, fp64x2 y) { return (fp64x2)hypot((fp64x2_Math)x, (fp64x2_Math)y); } - /* Trigonometry */ - inline fp64x2 sin (fp64x2 x) { return (fp64x2) sin ((fp64x2_Math)x); } - inline fp64x2 cos (fp64x2 x) { return (fp64x2) cos ((fp64x2_Math)x); } - inline fp64x2 tan (fp64x2 x) { return (fp64x2) tan ((fp64x2_Math)x); } - inline fp64x2 asin (fp64x2 x) { return (fp64x2)asin ((fp64x2_Math)x); } - inline fp64x2 acos (fp64x2 x) { return (fp64x2)acos ((fp64x2_Math)x); } - inline fp64x2 atan (fp64x2 x) { return (fp64x2)atan ((fp64x2_Math)x); } - inline fp64x2 sinh(fp64x2 x) { return (fp64x2) sinh((fp64x2_Math)x); } - inline fp64x2 cosh(fp64x2 x) { return (fp64x2) cosh((fp64x2_Math)x); } - inline fp64x2 tanh(fp64x2 x) { return (fp64x2) tanh((fp64x2_Math)x); } - inline fp64x2 asinh(fp64x2 x) { return (fp64x2)asinh((fp64x2_Math)x); } - inline fp64x2 acosh(fp64x2 x) { return (fp64x2)acosh((fp64x2_Math)x); } - inline fp64x2 atanh(fp64x2 x) { return (fp64x2)atanh((fp64x2_Math)x); } - inline fp64x2 atan2(fp64x2 y, fp64x2 x) { return (fp64x2)atan2((fp64x2_Math)y, (fp64x2_Math)x); } - inline void sincos(fp64x2 x, fp64x2* p_sin, fp64x2* p_cos) { - *p_sin = sin(x); - *p_cos = cos(x); - } - /* Logarithms and Exponents */ - inline fp64x2 log (fp64x2 x) { return (fp64x2)log ((fp64x2_Math)x); } - inline fp64x2 log1p(fp64x2 x) { return (fp64x2)log1p((fp64x2_Math)x); } - inline fp64x2 logb (fp64x2 x) { return (fp64x2)logb ((fp64x2_Math)x); } - inline fp64x2 log2 (fp64x2 x) { return (fp64x2)log2 ((fp64x2_Math)x); } - inline fp64x2 log10(fp64x2 x) { return (fp64x2)log10((fp64x2_Math)x); } - inline fp64x2 exp (fp64x2 x) { return (fp64x2)exp ((fp64x2_Math)x); } - inline fp64x2 expm1(fp64x2 x) { return (fp64x2)expm1((fp64x2_Math)x); } - inline fp64x2 exp2 (fp64x2 x) { return (fp64x2)exp2 ((fp64x2_Math)x); } - inline fp64x2 pow(fp64x2 x, fp64x2 y) { return (fp64x2)pow((fp64x2_Math)x, (fp64x2_Math)y); } - /* Rounding */ - // inline fp64x2 trunc(fp64x2 x) { return (fp64x2)trunc((fp64x2_Math)x); } - // inline fp64x2 floor(fp64x2 x) { return (fp64x2)floor((fp64x2_Math)x); } - // inline fp64x2 ceil (fp64x2 x) { return (fp64x2)ceil ((fp64x2_Math)x); } - // inline fp64x2 rint (fp64x2 x) { return (fp64x2)rint ((fp64x2_Math)x); } - // inline fp64x2 round(fp64x2 x) { return (fp64x2)round((fp64x2_Math)x); } - // inline long lrint (fp64x2 x) { return lrint ((fp64x2_Math)x); } - // inline long lround(fp64x2 x) { return lround((fp64x2_Math)x); } - // inline long long llrint (fp64x2 x) { return llrint ((fp64x2_Math)x); } - // inline long long llround(fp64x2 x) { return llround((fp64x2_Math)x); } - /* Integer and Remainder */ - inline fp64x2 fmod(fp64x2 x, fp64x2 y) { return (fp64x2)fmod((fp64x2_Math)x, (fp64x2_Math)y); } - inline fp64x2 modf(fp64x2 x, fp64x2* y) { - fp64x2_Math y_temp; - fp64x2 result = modf((fp64x2_Math)x, &y_temp); - *y = (fp64x2)y_temp; - return result; - } - // inline fp64x2 nearbyint(fp64x2 x) { return (fp64x2)nearbyint((fp64x2_Math)x); } - // Incorrect Function // inline fp64x2 nextafter(fp64x2 x, fp64x2 y) { return (fp64x2)nextafter((fp64x2_Math)x, (fp64x2_Math)y); } - inline fp64x2 remainder(fp64x2 x, fp64x2 y) { return (fp64x2)remainder((fp64x2_Math)x, (fp64x2_Math)y); } - inline fp64x2 remquo(fp64x2 x, fp64x2 y, int* quo) { return (fp64x2)remquo((fp64x2_Math)x, (fp64x2_Math)y, quo); } - /* Float Exponents */ - inline int ilogb(fp64x2 x) { return ilogb((fp64x2_Math)x); } - inline fp64x2 frexp (fp64x2 x, int* exp) { return (fp64x2)frexp ((fp64x2_Math)x, exp); } - inline fp64x2 ldexp (fp64x2 x, int exp) { return (fp64x2)ldexp ((fp64x2_Math)x, exp); } - inline fp64x2 scalbn (fp64x2 x, int exp) { return (fp64x2)scalbn ((fp64x2_Math)x, exp); } - inline fp64x2 scalbln(fp64x2 x, long exp) { return (fp64x2)scalbln((fp64x2_Math)x, exp); } - /* Tests */ - // inline bool signbit(fp64x2 x) { return (signbit((fp64x2_Math)x) != 0) ? true : false; } - // inline bool isfinite(fp64x2 x) { return (isfinite((fp64x2_Math)x) != 0) ? true : false; } - // inline bool isinf(fp64x2 x) { return (isinf((fp64x2_Math)x) != 0) ? true : false; } - // inline bool isnan(fp64x2 x) { return (isnan((fp64x2_Math)x) != 0) ? true : false; } - /* Transcendental Functions */ - inline fp64x2 erf (fp64x2 x) { return (fp64x2)erf ((fp64x2_Math)x); } - inline fp64x2 erfc(fp64x2 x) { return (fp64x2)erfc((fp64x2_Math)x); } - inline fp64x2 lgamma(fp64x2 x) { return (fp64x2)lgamma((fp64x2_Math)x); } - inline fp64x2 tgamma(fp64x2 x) { return (fp64x2)tgamma((fp64x2_Math)x); } - - /* Strings */ - - #include "double_FloatN_stringTo.hpp" - - inline Float64x2 stringTo_Float64x2(const char* nPtr, char** endPtr = nullptr) { - internal_double_FloatN_stringTo stringTo_func; - return stringTo_func.stringTo_FloatNx2(nPtr, endPtr); - } - - #include "double_FloatN_snprintf.hpp" - - #define PRIFloat64x2 "D" - #define PRIfp64x2 "D" - - /** - * @brief snprintf a singular Float64x2/fp64x2. - * Similar in functionallity to quadmath_snprintf. - * @note $ not supported. This function ignore additional - * format specifiers. - * @return -1 on encoding failure. Otherwise the total length of the - * string excluding the \0 terminator and ignoring the buffer size. - */ - inline int Float64x2_snprintf( - char* buf, size_t len, - const char* format, ... - ) { - va_list args; - va_start(args, format); - internal_double_FloatN_snprintf func_snprintf; - int ret_val = func_snprintf.FloatNx2_snprintf( - PRIFloat64x2, buf, len, - format, args - ); - va_end(args); - return ret_val; - } - -#endif /* DOUBLE_FLOAT64_SSE_H */ \ No newline at end of file diff --git a/src/floats/double_Float80.hpp b/src/floats/double_Float80.hpp index fb35dfc..3d70c14 100644 --- a/src/floats/double_Float80.hpp +++ b/src/floats/double_Float80.hpp @@ -27,25 +27,22 @@ typedef double fp64; typedef fp80 fp80x2_Math; #endif - /** * @brief Double-Float80 Dekker Float implementation. * Source: Creel "Double it Like Dekker" on YouTube. */ -class Float80x2 { -public: +struct Float80x2 { + fp80 hi; fp80 lo; - -private: /* Arithmetic */ - inline Float80x2 Dekker_Add( + static inline Float80x2 Dekker_Add( const Float80x2& x, const Float80x2& y - ) const { + ) { fp80 r_hi = x.hi + y.hi; - fp80 r_lo = 0.0; + fp80 r_lo = static_cast(0.0); if (fabsl(x.hi) > fabsl(y.hi)) { r_lo = x.hi - r_hi + y.hi + y.lo + x.lo; } else { @@ -58,11 +55,11 @@ class Float80x2 { return c; } - inline Float80x2 Dekker_Sub( + static inline Float80x2 Dekker_Sub( const Float80x2& x, const Float80x2& y - ) const { + ) { fp80 r_hi = x.hi - y.hi; - fp80 r_lo = 0.0; + fp80 r_lo = static_cast(0.0); if (fabsl(x.hi) > fabsl(y.hi)) { r_lo = x.hi - r_hi - y.hi - y.lo + x.lo; } else { @@ -77,7 +74,7 @@ class Float80x2 { static constexpr fp80 Dekker_Scale = 4294967297.0; // (2^ceil(64 / 2) + 1) - inline Float80x2 Dekker_Split(const fp80& x) const { + static inline Float80x2 Dekker_Split(const fp80& x) { fp80 p = x * Dekker_Scale; Float80x2 r; r.hi = (x - p) + p; @@ -96,9 +93,9 @@ class Float80x2 { // return r; // } - inline Float80x2 Dekker_Mul12( + static inline Float80x2 Dekker_Mul12( const fp80& x, const fp80& y - ) const { + ) { Float80x2 a = Dekker_Split(x); Float80x2 b = Dekker_Split(y); fp80 p = a.hi * b.hi; @@ -110,9 +107,9 @@ class Float80x2 { return r; } - inline Float80x2 Dekker_Mul( + static inline Float80x2 Dekker_Mul( const Float80x2& x, const Float80x2& y - ) const { + ) { Float80x2 t = Dekker_Mul12(x.hi, y.hi); fp80 c = x.hi * y.lo + x.lo * y.hi + t.lo; @@ -122,9 +119,9 @@ class Float80x2 { return r; } - inline Float80x2 Dekker_Div( + static inline Float80x2 Dekker_Div( const Float80x2& x, const Float80x2& y - ) const { + ) { Float80x2 u; u.hi = x.hi / y.hi; Float80x2 t = Dekker_Mul12(u.hi, y.hi); @@ -135,13 +132,39 @@ class Float80x2 { r.lo = u.hi - r.hi + l; return r; } + + static inline Float80x2 Dekker_Sqr12( + const fp80& x + ) { + Float80x2 a = Dekker_Split(x); + fp80 p = a.hi * a.hi; + fp80 q = static_cast(2.0) * a.hi * a.lo; + + Float80x2 r; + r.hi = p + q; + r.lo = p - r.hi + q + a.lo * a.lo; + return r; + } + + static inline Float80x2 Dekker_Sqr( + const Float80x2& x + ) { + Float80x2 t = Dekker_Sqr12(x.hi); + fp80 c = static_cast(2.0) * x.hi * x.lo + t.lo; + + Float80x2 r; + r.hi = t.hi + c; + r.lo = t.hi - r.hi + c; + return r; + } + /* Double-Single Arithmetic */ - inline Float80x2 Dekker_Add_Float80( + static inline Float80x2 Dekker_Add_Float80( const Float80x2& x, const fp80& y - ) const { + ) { fp80 r_hi = x.hi + y; - fp80 r_lo = 0.0f; + fp80 r_lo = static_cast(0.0); if (fabsl(x.hi) > fabsl(y)) { r_lo = x.hi - r_hi + y + x.lo; } else { @@ -154,11 +177,11 @@ class Float80x2 { return c; } - inline Float80x2 Dekker_Sub_Float80( + static inline Float80x2 Dekker_Sub_Float80( const Float80x2& x, const fp80& y - ) const { + ) { fp80 r_hi = x.hi - y; - fp80 r_lo = 0.0f; + fp80 r_lo = static_cast(0.0); if (fabsl(x.hi) > fabsl(y)) { r_lo = x.hi - r_hi - y + x.lo; } else { @@ -171,9 +194,9 @@ class Float80x2 { return c; } - inline Float80x2 Dekker_Mul_Float80( + static inline Float80x2 Dekker_Mul_Float80( const Float80x2& x, const fp80& y - ) const { + ) { Float80x2 t = Dekker_Mul12(x.hi, y); fp80 c = x.lo * y + t.lo; @@ -183,9 +206,9 @@ class Float80x2 { return r; } - inline Float80x2 Dekker_Div_Float80( + static inline Float80x2 Dekker_Div_Float80( const Float80x2& x, const fp80& y - ) const { + ) { Float80x2 u; u.hi = x.hi / y; Float80x2 t = Dekker_Mul12(u.hi, y); @@ -197,9 +220,9 @@ class Float80x2 { return r; } - inline Float80x2 Float80_Div_Dekker( + static inline Float80x2 Float80_Div_Dekker( const fp80& x, const Float80x2& y - ) const { + ) { Float80x2 u; u.hi = x / y.hi; Float80x2 t = Dekker_Mul12(u.hi, y.hi); @@ -211,8 +234,6 @@ class Float80x2 { return r; } -public: - /* Arithmetic */ inline Float80x2 operator+(const Float80x2 &value) const { @@ -415,12 +436,23 @@ typedef Float80x2 fp80x2; /* Math functions (Natively implemented) */ /* Arithmetic */ + inline fp80x2 fmax(fp80x2 x, fp80x2 y) { return (x > y) ? x : y; } + // inline fp80x2 fmax(fp80x2 x, fp80x2 y, fp80x2 z) { + // return (x > y) ? + // ((x > z) ? x : z) : + // ((y > z) ? y : z); + // } inline fp80x2 fmin(fp80x2 x, fp80x2 y) { return (x < y) ? x : y; } + // inline fp80x2 fmin(fp80x2 x, fp80x2 y, fp80x2 z) { + // return (x < y) ? + // ((x < z) ? x : z) : + // ((y < z) ? y : z); + // } inline fp80x2 fabs(fp80x2 x) { return (x < static_cast(0.0)) ? -x : x; } @@ -432,13 +464,44 @@ typedef Float80x2 fp80x2; } inline fp80x2 copysign(fp80x2 x, fp80x2 y) { return ( - (x < static_cast(0.0)) != (y < (static_cast(0.0))) + (x.hi < static_cast(0.0)) != (y.hi < (static_cast(0.0))) ) ? -x : x; } + /** @note This function name may change to square() or etc */ + inline fp80x2 sqr(fp80x2 x) { + return Float80x2::Dekker_Sqr(x); + } + inline fp80x2 sqrt(fp80x2 x) { + if (x == static_cast(0.0)) { + return x; + } + fp80x2 guess = (fp80x2)sqrt(x.hi); + return (guess + x / guess) * static_cast(0.5); + } + inline fp80x2 cbrt(fp80x2 x) { + if (x == static_cast(0.0)) { + return x; + } + fp80x2 guess = (fp80x2)cbrt(x.hi); + return ( + guess * static_cast(2.0) + (x) / Float80x2::Dekker_Sqr(guess) + ) / static_cast(3.0); + } + inline fp80x2 hypot(fp80x2 x, fp80x2 y) { + return sqrt( + Float80x2::Dekker_Sqr(x) + Float80x2::Dekker_Sqr(y) + ); + } + // inline fp80x2 hypot(fp80x2 x, fp80x2 y, fp80x2 z) { + // return sqrt( + // Float80x2::Dekker_Sqr(x) + Float80x2::Dekker_Sqr(y) + Float80x2::Dekker_Sqr(z) + // ); + // } /* Tests */ + inline bool signbit(fp80x2 x) { - return (x < static_cast(0.0)) ? true : false; + return (x.hi < static_cast(0.0)) ? true : false; } /** Returns true if both x.hi and x.lo are finite */ inline bool isfinite(fp80x2 x) { @@ -469,6 +532,7 @@ typedef Float80x2 fp80x2; } /* Comparison */ + inline bool isgreater(fp80x2 x, fp80x2 y) { return (x > y); } @@ -486,25 +550,26 @@ typedef Float80x2 fp80x2; } /* Rounding */ - inline fp80x2 trunc(fp80x2 x) { - fp80 frac_hi = x.hi - trunc(x.hi); - fp80 frac_lo = x.lo - trunc(x.lo); - fp80x2 int_hi = trunc(x.hi); - fp80x2 int_lo = trunc(x.lo); - // Sum in increasing order - fp80x2 trunc_all = static_cast(0.0); - trunc_all += ( + + inline fp80x2 trunc(fp80x2 x) { + fp80x2 int_hi = trunc(x.hi); + fp80x2 int_lo = trunc(x.lo); + fp80 frac_hi = x.hi - int_hi.hi; + fp80 frac_lo = x.lo - int_lo.hi; + // Sum in increasing order + fp80x2 trunc_all = static_cast(0.0); + trunc_all += ( (fp80x2)frac_hi + (fp80x2)frac_lo >= static_cast(1.0) ) ? static_cast(1.0) : static_cast(0.0); - trunc_all += int_lo; - trunc_all += int_hi; - return trunc_all; - } + trunc_all += int_lo; + trunc_all += int_hi; + return trunc_all; + } inline fp80x2 floor(fp80x2 x) { fp80x2 int_part = trunc(x); return ( x < static_cast(0.0) && int_part != x - ) ? int_part : int_part - static_cast(1.0); + ) ? int_part - static_cast(1.0) : int_part; } inline fp80x2 ceil(fp80x2 x) { fp80x2 int_part = trunc(x); @@ -553,10 +618,36 @@ typedef Float80x2 fp80x2; } /* Integer and Remainder */ + + inline fp80x2 modf(fp80x2 x, fp80x2* int_part) { + fp80x2 trunc_part = trunc(x); + if (int_part != nullptr) { + *int_part = trunc_part; + } + return x - trunc_part; + } inline fp80x2 nearbyint(fp80x2 x) { return rint(x); } + /* Float Exponents */ + + inline fp80x2 ldexp(fp80x2 x, int exp) { + x.hi = ldexp(x.hi, exp); + x.lo = isfinite(x.hi) ? ldexp(x.lo, exp) : x.hi; + return x; + } + inline fp80x2 scalbn(fp80x2 x, int exp) { + x.hi = scalbn(x.hi, exp); + x.lo = isfinite(x.hi) ? scalbn(x.lo, exp) : x.hi; + return x; + } + inline fp80x2 scalbln(fp80x2 x, long exp) { + x.hi = scalbln(x.hi, exp); + x.lo = isfinite(x.hi) ? scalbln(x.lo, exp) : x.hi; + return x; + } + /* Math overloads (Casts to other types) */ /* Arithmetic */ @@ -566,9 +657,9 @@ typedef Float80x2 fp80x2; // inline fp80x2 fdim(fp80x2 x, fp80x2 y) { return (fp80x2)fdim((fp80x2_Math)x, (fp80x2_Math)y); } // inline fp80x2 fma(fp80x2 x, fp80x2 y, fp80x2 z) { return (fp80x2)fma((fp80x2_Math)x, (fp80x2_Math)y, (fp80x2_Math)z); } // inline fp80x2 copysign(fp80x2 x, fp80x2 y) { return (fp80x2)copysign((fp80x2_Math)x, (fp80x2_Math)y); } - inline fp80x2 sqrt(fp80x2 x) { return (fp80x2)sqrt((fp80x2_Math)x); } - inline fp80x2 cbrt(fp80x2 x) { return (fp80x2)cbrt((fp80x2_Math)x); } - inline fp80x2 hypot(fp80x2 x, fp80x2 y) { return (fp80x2)hypot((fp80x2_Math)x, (fp80x2_Math)y); } + // inline fp80x2 sqrt(fp80x2 x) { return (fp80x2)sqrt((fp80x2_Math)x); } + // inline fp80x2 cbrt(fp80x2 x) { return (fp80x2)cbrt((fp80x2_Math)x); } + // inline fp80x2 hypot(fp80x2 x, fp80x2 y) { return (fp80x2)hypot((fp80x2_Math)x, (fp80x2_Math)y); } /* Trigonometry */ inline fp80x2 sin (fp80x2 x) { return (fp80x2) sin ((fp80x2_Math)x); } inline fp80x2 cos (fp80x2 x) { return (fp80x2) cos ((fp80x2_Math)x); } @@ -609,12 +700,12 @@ typedef Float80x2 fp80x2; // inline long long llround(fp80x2 x) { return llround((fp80x2_Math)x); } /* Integer and Remainder */ inline fp80x2 fmod(fp80x2 x, fp80x2 y) { return (fp80x2)fmod((fp80x2_Math)x, (fp80x2_Math)y); } - inline fp80x2 modf(fp80x2 x, fp80x2* y) { - fp80x2_Math y_temp; - fp80x2 result = modf((fp80x2_Math)x, &y_temp); - *y = (fp80x2)y_temp; - return result; - } + // inline fp80x2 modf(fp80x2 x, fp80x2* y) { + // fp80x2_Math y_temp; + // fp80x2 result = modf((fp80x2_Math)x, &y_temp); + // *y = (fp80x2)y_temp; + // return result; + // } // inline fp80x2 nearbyint(fp80x2 x) { return (fp80x2)nearbyint((fp80x2_Math)x); } // Incorrect Function // inline fp80x2 nextafter(fp80x2 x, fp80x2 y) { return (fp80x2)nextafter((fp80x2_Math)x, (fp80x2_Math)y); } inline fp80x2 remainder(fp80x2 x, fp80x2 y) { return (fp80x2)remainder((fp80x2_Math)x, (fp80x2_Math)y); } @@ -622,9 +713,9 @@ typedef Float80x2 fp80x2; /* Float Exponents */ inline int ilogb(fp80x2 x) { return ilogb((fp80x2_Math)x); } inline fp80x2 frexp (fp80x2 x, int* exp) { return (fp80x2)frexp ((fp80x2_Math)x, exp); } - inline fp80x2 ldexp (fp80x2 x, int exp) { return (fp80x2)ldexp ((fp80x2_Math)x, exp); } - inline fp80x2 scalbn (fp80x2 x, int exp) { return (fp80x2)scalbn ((fp80x2_Math)x, exp); } - inline fp80x2 scalbln(fp80x2 x, long exp) { return (fp80x2)scalbln((fp80x2_Math)x, exp); } + // inline fp80x2 ldexp (fp80x2 x, int exp) { return (fp80x2)ldexp ((fp80x2_Math)x, exp); } + // inline fp80x2 scalbn (fp80x2 x, int exp) { return (fp80x2)scalbn ((fp80x2_Math)x, exp); } + // inline fp80x2 scalbln(fp80x2 x, long exp) { return (fp80x2)scalbln((fp80x2_Math)x, exp); } /* Tests */ // inline bool signbit(fp80x2 x) { return (signbit((fp80x2_Math)x) != 0) ? true : false; } // inline bool isfinite(fp80x2 x) { return (isfinite((fp80x2_Math)x) != 0) ? true : false; } @@ -675,4 +766,4 @@ typedef Float80x2 fp80x2; #endif -#endif /* DOUBLE_FLOAT64_HPP */ \ No newline at end of file +#endif /* DOUBLE_FLOAT80_HPP */ \ No newline at end of file diff --git a/src/floats/double_FloatN_snprintf.hpp b/src/floats/double_FloatN_snprintf.hpp index 71f60a3..b73aae3 100644 --- a/src/floats/double_FloatN_snprintf.hpp +++ b/src/floats/double_FloatN_snprintf.hpp @@ -1,6 +1,6 @@ /* -** Author: zerico2005 (2023 - 2024) -** Project: ABS-Fractal-Explorer +** Author: zerico2005 (2024) +** Project: LIB-Dekker-Float ** License: MIT License ** A copy of the MIT License should be included with ** this project. If not, see https://opensource.org/license/MIT @@ -195,6 +195,53 @@ class internal_double_FloatN_snprintf { return '?'; } + // /** + // * @brief Increments decimial/hexadecimal digit + // * @note Assumes base is either 10 or 16 for %f, %F, %a, and %A + // * @returns true if an overflow occured. + // */ + // static bool inc_digit_from_base( + // char* digit, int base, bool upperCase + // ) { + // printf("L: %d %02X %p\n", __LINE__, *digit, digit); + // if (digit == nullptr) { return false; } + // assert(base == 10 || base == 16); + // if (base == 16) { + // printf("L: %d\n", __LINE__); + // if ( + // (upperCase == true && *digit == 'F') || + // (upperCase == false && *digit == 'f') + // ) { + // *digit = '0'; + // return true; // Overflow + // } + // if (*digit == '9') { + // *digit = (upperCase ? 'A' : 'a'); + // return false; + // } + // if ( + // (upperCase == true && *digit >= 'A' && *digit <= 'E') || + // (upperCase == false && *digit >= 'a' && *digit <= 'e') + // ) { + // (*digit)++; + // return false; + // } + // } else { + // printf("L: %d\n", __LINE__); + // if (*digit == '9') { + // *digit = '0'; + // return true; // Overflow + // } + // } + // printf("L: %d\n", __LINE__); + // if (*digit >= '0' && *digit <= '8') { + // (*digit)++; + // return false; + // } + // printf("L: %d\n", __LINE__); + // return false; + // } + static std::string FloatNx2_write( const FloatNx2_format_param& param, FloatNx2 value, const int base, const bool upperCase @@ -255,8 +302,25 @@ class internal_double_FloatN_snprintf { { /* Rounds the last digit */ value_frac *= (FloatNx2)base; FloatNx2 digit_temp = round(value_frac); - char digit = get_digit_from_base((int)digit_temp, base, upperCase); - frac_digits += digit; + if ((int)digit_temp == base) { + /* round() + char* round_ptr = (char*)frac_digits.c_str(); + frac_digits += '0'; + while (inc_digit_from_base(round_ptr, base, upperCase)) { + round_ptr--; + } + */ + // /* trunc() + if (base == 16) { + frac_digits += upperCase ? 'F' : 'f'; + } else { + frac_digits += '9'; + } + // */ + } else { + char digit = get_digit_from_base((int)digit_temp, base, upperCase); + frac_digits += digit; + } } str += frac_digits; @@ -297,7 +361,6 @@ class internal_double_FloatN_snprintf { if (fm_ptr == nullptr || *fm_ptr == '\0') { return -1; // Invalid format } - const char* const fm_start = fm_ptr; FloatNx2_format_param param; @@ -384,7 +447,7 @@ class internal_double_FloatN_snprintf { strncpy(buf, output_str.c_str(), len); // va_end(args); - return (int)strlen(format); + return (int)strlen(output_str.c_str()); } }; diff --git a/src/floats/double_FloatN_stringTo.hpp b/src/floats/double_FloatN_stringTo.hpp index 27b51b8..b4fc628 100644 --- a/src/floats/double_FloatN_stringTo.hpp +++ b/src/floats/double_FloatN_stringTo.hpp @@ -1,6 +1,6 @@ /* -** Author: zerico2005 (2023 - 2024) -** Project: ABS-Fractal-Explorer +** Author: zerico2005 (2024) +** Project: LIB-Dekker-Float ** License: MIT License ** A copy of the MIT License should be included with ** this project. If not, see https://opensource.org/license/MIT @@ -43,9 +43,12 @@ class internal_double_FloatN_stringTo { static const char* get_exponent(const char* ptr, int& exponent) { exponent = 0; - bool exp_sign = (*ptr == '-') ? true : false; + bool exp_sign = false; if (*ptr == 'E' || *ptr == 'e') { ptr++; + if (*ptr == '-') { + exp_sign = true; + } if (*ptr == '+' || *ptr == '-') { ptr++; } diff --git a/src/render_CPU/frac_Multi_Float64x2_AVX.cpp b/src/render_CPU/frac_Multi_Float64x2_AVX.cpp index ca43ede..ceb8508 100644 --- a/src/render_CPU/frac_Multi_Float64x2_AVX.cpp +++ b/src/render_CPU/frac_Multi_Float64x2_AVX.cpp @@ -253,9 +253,9 @@ } // Bits 0-2 will flip signage - const __m256dx2 s1 = (f[0]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s2 = (f[1]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s3 = (f[2]) ? _mm256x2_set1_pd(-2.0) : _mm256x2_set1_pd(2.0); + const __m256d s1 = (f[0]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s2 = (f[1]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s3 = (f[2]) ? _mm256_set1_pd(-2.0) : _mm256_set1_pd(2.0); // Bits 3-7 will apply fabs() via a mask const __m256dx2 zr1_mask = (f[3]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); @@ -271,15 +271,15 @@ zi2 = _mm256x2_and_pdx2(zi2_mask, zi); // Calculates the new zr and zi - zr = _mm256x2_add_pdx2(_mm256x2_and_pdx2(zr_mask, _mm256x2_mul_pdx2(s1, + zr = _mm256x2_add_pdx2(_mm256x2_and_pdx2(zr_mask, _mm256x2_mul_pd_pdx2(s1, _mm256x2_sub_pdx2( _mm256x2_mul_pdx2(zr1, zr), - _mm256x2_mul_pdx2(s2, + _mm256x2_mul_pd_pdx2(s2, _mm256x2_mul_pdx2(zi1, zi) ) ) )), cr); - zi = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr2, zi2), s3), ci); + zi = _mm256x2_add_pdx2(_mm256x2_mul_pdx2_pd(_mm256x2_mul_pdx2(zr2, zi2), s3), ci); zs = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(zr, zr), _mm256x2_mul_pdx2(zi, zi)); @@ -295,12 +295,12 @@ for (uint8_t q = 0; q < 14; q++) { f[q] = ((param.formula >> q) & 1) ? true : false; } - const __m256dx2 s1 = (f[0]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s2 = (f[1]) ? _mm256x2_set1_pd(-3.0) : _mm256x2_set1_pd(3.0); - const __m256dx2 s3 = (f[2]) ? _mm256x2_set1_pd(-3.0) : _mm256x2_set1_pd(3.0); - const __m256dx2 s4 = (f[3]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s5 = (f[4]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s6 = (f[5]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); + const __m256d s1 = (f[0]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s2 = (f[1]) ? _mm256_set1_pd(-3.0) : _mm256_set1_pd(3.0); + const __m256d s3 = (f[2]) ? _mm256_set1_pd(-3.0) : _mm256_set1_pd(3.0); + const __m256d s4 = (f[3]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s5 = (f[4]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s6 = (f[5]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); const __m256dx2 zr1_mask = (f[ 6]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); const __m256dx2 zi1_mask = (f[ 7]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); @@ -320,16 +320,16 @@ zr3 = _mm256x2_and_pdx2(zr3_mask, zr); zi3 = _mm256x2_and_pdx2(zi3_mask, zi); - temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s5, _mm256x2_and_pdx2(zr_mask, + temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s5, _mm256x2_and_pdx2(zr_mask, _mm256x2_sub_pdx2( - _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(s1, zr1), _mm256x2_mul_pdx2(zr , zr)), - _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(s2, zr2), _mm256x2_mul_pdx2(zi1, zi)) + _mm256x2_mul_pdx2(_mm256x2_mul_pd_pdx2(s1, zr1), _mm256x2_mul_pdx2(zr , zr)), + _mm256x2_mul_pdx2(_mm256x2_mul_pd_pdx2(s2, zr2), _mm256x2_mul_pdx2(zi1, zi)) ) )), cr); - zi = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s6, _mm256x2_and_pdx2(zi_mask, + zi = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s6, _mm256x2_and_pdx2(zi_mask, _mm256x2_sub_pdx2( - _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(s3, zr3), _mm256x2_mul_pdx2(zr, zi2)), - _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(s4, zi3), _mm256x2_mul_pdx2(zi, zi )) + _mm256x2_mul_pdx2(_mm256x2_mul_pd_pdx2(s3, zr3), _mm256x2_mul_pdx2(zr, zi2)), + _mm256x2_mul_pdx2(_mm256x2_mul_pd_pdx2(s4, zi3), _mm256x2_mul_pdx2(zi, zi )) ) )), ci); zr = temp_zr; @@ -348,13 +348,13 @@ for (uint8_t q = 0; q < 17; q++) { f[q] = ((param.formula >> q) & 1) ? true : false; } - const __m256dx2 s1 = (f[0]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s2 = (f[1]) ? _mm256x2_set1_pd(-6.0) : _mm256x2_set1_pd(6.0); - const __m256dx2 s3 = (f[2]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s4 = (f[3]) ? _mm256x2_set1_pd(-4.0) : _mm256x2_set1_pd(4.0); - const __m256dx2 s5 = (f[4]) ? _mm256x2_set1_pd(-4.0) : _mm256x2_set1_pd(4.0); - const __m256dx2 s6 = (f[5]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); - const __m256dx2 s7 = (f[6]) ? _mm256x2_set1_pd(-1.0) : _mm256x2_set1_pd(1.0); + const __m256d s1 = (f[0]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s2 = (f[1]) ? _mm256_set1_pd(-6.0) : _mm256_set1_pd(6.0); + const __m256d s3 = (f[2]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s4 = (f[3]) ? _mm256_set1_pd(-4.0) : _mm256_set1_pd(4.0); + const __m256d s5 = (f[4]) ? _mm256_set1_pd(-4.0) : _mm256_set1_pd(4.0); + const __m256d s6 = (f[5]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); + const __m256d s7 = (f[6]) ? _mm256_set1_pd(-1.0) : _mm256_set1_pd(1.0); const __m256dx2 zr1_mask = (f[ 7]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); const __m256dx2 zi1_mask = (f[ 8]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); @@ -378,27 +378,27 @@ zr4 = _mm256x2_and_pdx2(zr4_mask, zr); zi4 = _mm256x2_and_pdx2(zi4_mask, zi); - temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s6, _mm256x2_and_pdx2(zr_mask, + temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s6, _mm256x2_and_pdx2(zr_mask, _mm256x2_add_pdx2( _mm256x2_sub_pdx2( - _mm256x2_mul_pdx2(s1, + _mm256x2_mul_pd_pdx2(s1, _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr1, zr ), _mm256x2_mul_pdx2(zr, zr)) ), - _mm256x2_mul_pdx2(s2, + _mm256x2_mul_pd_pdx2(s2, _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr2, zi1), _mm256x2_mul_pdx2(zr, zi)) ) ), - _mm256x2_mul_pdx2(s3, + _mm256x2_mul_pd_pdx2(s3, _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zi2, zi), _mm256x2_mul_pdx2(zi, zi)) ) ) )), cr); - zi = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s7, _mm256x2_and_pdx2(zi_mask, + zi = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s7, _mm256x2_and_pdx2(zi_mask, _mm256x2_sub_pdx2( - _mm256x2_mul_pdx2(s4, + _mm256x2_mul_pd_pdx2(s4, _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr3, zi3), _mm256x2_mul_pdx2(zr, zr)) ), - _mm256x2_mul_pdx2(s5, + _mm256x2_mul_pd_pdx2(s5, _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr4, zi4), _mm256x2_mul_pdx2(zi, zi)) ) ) @@ -430,14 +430,14 @@ for (uint8_t q = 18; q < 20; q++) { /* 18-19 */ fO[q - 16] = ((param.formula >> q) & 1) ? true : false; } - const __m256dx2 s1 = (fS[0]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); - const __m256dx2 s2 = (fS[1]) ? _mm256x2_set1_pd(-10.0) : _mm256x2_set1_pd(10.0); - const __m256dx2 s3 = (fS[2]) ? _mm256x2_set1_pd(-5.0 ) : _mm256x2_set1_pd(5.0 ); - const __m256dx2 s4 = (fS[3]) ? _mm256x2_set1_pd(-5.0 ) : _mm256x2_set1_pd(5.0 ); - const __m256dx2 s5 = (fS[4]) ? _mm256x2_set1_pd(-10.0) : _mm256x2_set1_pd(10.0); - const __m256dx2 s6 = (fS[5]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); - const __m256dx2 s7 = (fO[0]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); - const __m256dx2 s8 = (fO[1]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); + const __m256d s1 = (fS[0]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); + const __m256d s2 = (fS[1]) ? _mm256_set1_pd(-10.0) : _mm256_set1_pd(10.0); + const __m256d s3 = (fS[2]) ? _mm256_set1_pd(-5.0 ) : _mm256_set1_pd(5.0 ); + const __m256d s4 = (fS[3]) ? _mm256_set1_pd(-5.0 ) : _mm256_set1_pd(5.0 ); + const __m256d s5 = (fS[4]) ? _mm256_set1_pd(-10.0) : _mm256_set1_pd(10.0); + const __m256d s6 = (fS[5]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); + const __m256d s7 = (fO[0]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); + const __m256d s8 = (fO[1]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); const __m256dx2 zr1_mask = (fA[0]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); const __m256dx2 zi1_mask = (fA[1]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); @@ -465,38 +465,38 @@ zr5 = _mm256x2_and_pdx2(zr5_mask, zr); zi5 = _mm256x2_and_pdx2(zi5_mask, zi); - temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s7, _mm256x2_and_pdx2(zr_mask , + temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s7, _mm256x2_and_pdx2(zr_mask , _mm256x2_add_pdx2( _mm256x2_sub_pdx2( _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s1, zr1), + _mm256x2_mul_pd_pdx2(s1, zr1), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr, zr), _mm256x2_mul_pdx2(zr , zr)) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s2, _mm256x2_mul_pdx2(zr2, zi1)), + _mm256x2_mul_pd_pdx2(s2, _mm256x2_mul_pdx2(zr2, zi1)), _mm256x2_mul_pdx2(zi, _mm256x2_mul_pdx2(zr , zr )) ) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s3, _mm256x2_mul_pdx2(zr3, zi2)), + _mm256x2_mul_pd_pdx2(s3, _mm256x2_mul_pdx2(zr3, zi2)), _mm256x2_mul_pdx2(zi, _mm256x2_mul_pdx2(zi , zi )) ) ) )), cr); - zi = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s8, _mm256x2_and_pdx2(zi_mask, + zi = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s8, _mm256x2_and_pdx2(zi_mask, _mm256x2_add_pdx2( _mm256x2_sub_pdx2( _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s4 , _mm256x2_mul_pdx2(zr4, zr)), + _mm256x2_mul_pd_pdx2(s4 , _mm256x2_mul_pdx2(zr4, zr)), _mm256x2_mul_pdx2(zi3, _mm256x2_mul_pdx2(zr , zr)) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s5 , _mm256x2_mul_pdx2(zr5, zr)), + _mm256x2_mul_pd_pdx2(s5 , _mm256x2_mul_pdx2(zr5, zr)), _mm256x2_mul_pdx2(zi4, _mm256x2_mul_pdx2(zi , zi)) ) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s6 , zi5), + _mm256x2_mul_pd_pdx2(s6 , zi5), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zi, zi), _mm256x2_mul_pdx2(zi , zi)) ) ) @@ -529,15 +529,15 @@ fO[q - 19] = ((param.formula >> q) & 1) ? true : false; } - const __m256dx2 s1 = (fS[0]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); - const __m256dx2 s2 = (fS[1]) ? _mm256x2_set1_pd(-15.0) : _mm256x2_set1_pd(15.0); - const __m256dx2 s3 = (fS[2]) ? _mm256x2_set1_pd(-15.0) : _mm256x2_set1_pd(15.0); - const __m256dx2 s4 = (fS[3]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); - const __m256dx2 s5 = (fS[4]) ? _mm256x2_set1_pd(-6.0 ) : _mm256x2_set1_pd(6.0 ); - const __m256dx2 s6 = (fS[5]) ? _mm256x2_set1_pd(-20.0) : _mm256x2_set1_pd(20.0); - const __m256dx2 s7 = (fS[6]) ? _mm256x2_set1_pd(-6.0 ) : _mm256x2_set1_pd(6.0 ); - const __m256dx2 s8 = (fO[0]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); - const __m256dx2 s9 = (fO[1]) ? _mm256x2_set1_pd(-1.0 ) : _mm256x2_set1_pd(1.0 ); + const __m256d s1 = (fS[0]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); + const __m256d s2 = (fS[1]) ? _mm256_set1_pd(-15.0) : _mm256_set1_pd(15.0); + const __m256d s3 = (fS[2]) ? _mm256_set1_pd(-15.0) : _mm256_set1_pd(15.0); + const __m256d s4 = (fS[3]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); + const __m256d s5 = (fS[4]) ? _mm256_set1_pd(-6.0 ) : _mm256_set1_pd(6.0 ); + const __m256d s6 = (fS[5]) ? _mm256_set1_pd(-20.0) : _mm256_set1_pd(20.0); + const __m256d s7 = (fS[6]) ? _mm256_set1_pd(-6.0 ) : _mm256_set1_pd(6.0 ); + const __m256d s8 = (fO[0]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); + const __m256d s9 = (fO[1]) ? _mm256_set1_pd(-1.0 ) : _mm256_set1_pd(1.0 ); const __m256dx2 zr1_mask = (fA[ 0]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); const __m256dx2 zi1_mask = (fA[ 1]) ? _mm256x2_set1_epi64x((int64_t)0x7FFFFFFFFFFFFFFF) : _mm256x2_set1_epi64x((int64_t)0xFFFFFFFFFFFFFFFF); @@ -569,44 +569,44 @@ zr6 = _mm256x2_and_pdx2(zr6_mask, zr); zi6 = _mm256x2_and_pdx2(zi6_mask, zi); - temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s8, _mm256x2_and_pdx2(zr_mask , + temp_zr = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s8, _mm256x2_and_pdx2(zr_mask , _mm256x2_sub_pdx2( _mm256x2_add_pdx2( _mm256x2_sub_pdx2( _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s1, _mm256x2_mul_pdx2(zr1, zr)), + _mm256x2_mul_pd_pdx2(s1, _mm256x2_mul_pdx2(zr1, zr)), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr, zr), _mm256x2_mul_pdx2(zr, zr)) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s2, _mm256x2_mul_pdx2(zr2, zi1)), + _mm256x2_mul_pd_pdx2(s2, _mm256x2_mul_pdx2(zr2, zi1)), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr, zr), _mm256x2_mul_pdx2(zr, zi)) ) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s3, _mm256x2_mul_pdx2(zr3, zi2)), + _mm256x2_mul_pd_pdx2(s3, _mm256x2_mul_pdx2(zr3, zi2)), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr, zi), _mm256x2_mul_pdx2(zi, zi)) ) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s4, _mm256x2_mul_pdx2(zi3, zi)), + _mm256x2_mul_pd_pdx2(s4, _mm256x2_mul_pdx2(zi3, zi)), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zi, zi), _mm256x2_mul_pdx2(zi, zi)) ) ) )), cr); - zi = _mm256x2_add_pdx2(_mm256x2_mul_pdx2(s9, _mm256x2_and_pdx2(zi_mask, + zi = _mm256x2_add_pdx2(_mm256x2_mul_pd_pdx2(s9, _mm256x2_and_pdx2(zi_mask, _mm256x2_add_pdx2( _mm256x2_sub_pdx2( _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s5, _mm256x2_mul_pdx2(zr4, zi4)), + _mm256x2_mul_pd_pdx2(s5, _mm256x2_mul_pdx2(zr4, zi4)), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr, zr), _mm256x2_mul_pdx2(zr, zr)) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s6, _mm256x2_mul_pdx2(zr5, zi5)), + _mm256x2_mul_pd_pdx2(s6, _mm256x2_mul_pdx2(zr5, zi5)), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zr, zr), _mm256x2_mul_pdx2(zi, zi)) ) ), _mm256x2_mul_pdx2( - _mm256x2_mul_pdx2(s7, _mm256x2_mul_pdx2(zr6, zi6)), + _mm256x2_mul_pd_pdx2(s7, _mm256x2_mul_pdx2(zr6, zi6)), _mm256x2_mul_pdx2(_mm256x2_mul_pdx2(zi, zi), _mm256x2_mul_pdx2(zi, zi)) ) )