diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,19 @@ # Changelog for fp-ieee +## Version 0.1.0.7 (2026-09-28) + +* Fix bugs with `Float128`. +* Fix `nextTowardZeroHalf`. +* Fix correctness issue with the generic FMA. +* Fix `remainder` when the result is zero. +* Fix `doubleToHalf` on the F16C configuration. +* Make use of RISC-V instructions. Some features need `rva22u64` or `rva23u64` package flags. +* Support GHC 10.0. + +## Version 0.1.0.6 (2025-12-29) + +* Support GHC 9.14. + ## Version 0.1.0.5 (2024-12-15) * Support GHC 9.10/9.12. diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright ARATA Mizuki (c) 2020-2024 +Copyright ARATA Mizuki (c) 2020-2026 All rights reserved. diff --git a/cbits/canonicalize.c b/cbits/canonicalize.c --- a/cbits/canonicalize.c +++ b/cbits/canonicalize.c @@ -8,7 +8,7 @@ float hs_canonicalizeFloat(float x) { - asm volatile("mulss %1, %0" : "+x"(x) : "x"(1.0f)); + __asm__ __volatile__("mulss %1, %0" : "+x"(x) : "x"(1.0f)); return x; /* Clang optimizes away this: @@ -22,7 +22,7 @@ } double hs_canonicalizeDouble(double x) { - asm volatile("mulsd %1, %0" : "+x"(x) : "x"(1.0)); + __asm__ __volatile__("mulsd %1, %0" : "+x"(x) : "x"(1.0)); return x; /* Clang optimizes away this: @@ -39,12 +39,25 @@ float hs_canonicalizeFloat(float x) { - asm volatile("fmul %s0, %s0, %s1" : "+w"(x) : "w"(1.0f)); + __asm__ __volatile__("fmul %s0, %s0, %s1" : "+w"(x) : "w"(1.0f)); return x; } double hs_canonicalizeDouble(double x) { - asm volatile("fmul %d0, %d0, %d1" : "+w"(x) : "w"(1.0)); + __asm__ __volatile__("fmul %d0, %d0, %d1" : "+w"(x) : "w"(1.0)); + return x; +} + +#elif defined(__riscv) + +float hs_canonicalizeFloat(float x) +{ + __asm__ __volatile__("fmul.s %0, %0, %1" : "+f"(x) : "f"(1.0f)); + return x; +} +double hs_canonicalizeDouble(double x) +{ + __asm__ __volatile__("fmul.d %0, %0, %1" : "+f"(x) : "f"(1.0)); return x; } diff --git a/cbits/half.c b/cbits/half.c --- a/cbits/half.c +++ b/cbits/half.c @@ -1,11 +1,50 @@ #include // uint16_t +#include #include +/* + * binary16 + * p = 11 + * emin = -14 + * emax = 15 + * floatRange _ = (-13,16) + * minPositive = 0x1p-24 + * minPositiveNormal = 0x1p-14 + * maxFinite = 0x1.ffcp15 + * exponentBias = 15 + */ + +/* + * binary32 + * p = 24 + * emin = -126 + * emax = 127 + * floatRange _ = (-125,128) + * minPositive = 0x1p-149 + * minPositiveNormal = 0x1p-126 + * maxFinite = 0x1.fffffep127 + * exponentBias = 127 + */ + +/* + * binary64 + * p = 53 + * emin = -1022 + * emax = 1023 + * floatRange _ = (-1021,1024) + * minPositive = 0x1p-1074 + * minPositiveNormal = 0x1p-1022 + * maxFinite = 0x1.ffff_ffff_ffff_fp1023 + * exponentBias = 1023 + */ + +// TODO: Use AVX-512 FP16 if available + #if defined(__F16C__) // x86 F16C #include -uint16_t hs_fastFloatToHalf(float f) +uint16_t hs_fp_ieee_floatToHalf(float f) { __m128 x = _mm_set_ss(f); union { @@ -17,7 +56,7 @@ return u.c; } -float hs_fastHalfToFloat(uint16_t c) +float hs_fp_ieee_halfToFloat(uint16_t c) { union { __m128i v; @@ -30,37 +69,71 @@ return d; } -// Is this really faster than bit manipulation? -uint16_t hs_fastDoubleToHalf(double d) +uint16_t hs_fp_ieee_doubleToHalf(double d) { - float f = (float)d; - if ((double)f != d && isfinite(f)) { - // The conversion was inexact. - // Use "round-to-odd" trick. - union { - float x; - struct { - // little-endian - unsigned mant: 23; - unsigned exp: 8; - unsigned sign: 1; - }; - } w; - w.x = f; - w.mant |= 1; - f = w.x; + uint64_t q; + memcpy(&q, &d, sizeof(d)); + const uint64_t F64_SIGN_MASK = UINT64_C(0x8000000000000000); + const uint64_t F64_EXP_MASK = UINT64_C(0x7ff0000000000000); + const uint64_t F64_MANT_MASK = UINT64_C(0x000fffffffffffff); + const uint64_t UPPER_MANT_MASK = UINT64_C(0x000ffc0000000000); + const uint64_t LOWER_MANT_MASK = UINT64_C(0x000003ffffffffff); + const uint64_t TIE = UINT64_C(0x0000020000000000); + uint16_t f16sign = (uint16_t)((q & F64_SIGN_MASK) >> (63 - 15)); + uint64_t biasedExp = (q & F64_EXP_MASK) >> 52; + // biasedExp - 1023 <= -26: round to zero + // biasedExp - 1023 == -25: if abs d == 0x1p-25 then +-0 else +-0x1p-24 + // -24 <= biasedExp - 1023 <= -16: subnormal + // biasedExp - 1023 == -15: if abs d >= 0x1.ffep-15 then +-0x1p-14 else subnormal + // -14 <= biasedExp - 1023 <= 14: normal + // biasedExp - 1023 == 15: if abs d >= 0x1.ffep15 then +-infinity else normal + // 16 <= biasedExp - 1023 : infinity + if (biasedExp < 1023 - 25) { + // biasedExp - 1023 < -25: zero + return f16sign | 0; + } else if (biasedExp < 1023 - 14) { + // biasedExp - 1023 < -14: zero / subnormal / minPositiveNormal + // 2^(biasedExp - 1023) <= abs d < 2^(biasedExp - 1022) + // ulp = 0x1p-24 + // 0b1XX..X * 0x1p-24 + // ^^ ^ + // || +- 2^(-24) + // | \ : (biasedExp - 1024) - (-24) + 1 = biasedExp - (1024 - 24 - 1) + // \ ---- 2^(biasedExp - 1024) + // ----- 2^(biasedExp - 1023) + int bitWidth = biasedExp - (1024 - 24 - 1); + uint64_t mant = (F64_MANT_MASK + 1) | (q & F64_MANT_MASK); + uint16_t f16mant = (uint16_t)(mant >> (52 - bitWidth)); + uint64_t lower = mant & (((LOWER_MANT_MASK + 1) << (10 - bitWidth)) - 1); + uint64_t tie = TIE << (10 - bitWidth); + if (lower < tie || (lower == tie && (f16mant & 1) == 0)) { + return f16sign | f16mant; + } else { + return (f16sign | f16mant) + 1; // may be +-minPositiveNormal + } + } else if (biasedExp < 1023 + 16) { + // biasedExp - 1023 < 16: normal / infinity + uint16_t f16mant = (uint16_t)((q & UPPER_MANT_MASK) >> (53 - 11)); + uint16_t f16exp = (uint16_t)((biasedExp - (1023 - 15)) << 10); + uint64_t lower = q & LOWER_MANT_MASK; + if (lower < TIE || (lower == TIE && (f16mant & 1) == 0)) { + return f16sign | f16exp | f16mant; + } else { + return (f16sign | f16exp | f16mant) + 1; // may be +-infinity + } + } else { + // infinity or NaN + if (isnan(d)) { + // keep the sign bit + // discard the payload + return f16sign | 0x7e00; + } else { + return f16sign | 0x7c00; + } } - __m128 x = _mm_set_ss(f); - union { - __m128i v; - uint16_t c; - } u; - // A floating-point exception can be raised - u.v = _mm_cvtps_ph(x, _MM_FROUND_TO_NEAREST_INT); // VCVTPS2PH - return u.c; } -double hs_fastHalfToDouble(uint16_t c) +double hs_fp_ieee_halfToDouble(uint16_t c) { union { __m128i v; @@ -77,7 +150,7 @@ // Let's hope _Float16 is available -uint16_t hs_fastFloatToHalf(float x) +uint16_t hs_fp_ieee_floatToHalf(float x) { union { _Float16 f; @@ -87,7 +160,7 @@ return u.u; } -float hs_fastHalfToFloat(uint16_t x) +float hs_fp_ieee_halfToFloat(uint16_t x) { union { _Float16 f; @@ -97,7 +170,7 @@ return (float)u.f; } -uint16_t hs_fastDoubleToHalf(double x) +uint16_t hs_fp_ieee_doubleToHalf(double x) { union { _Float16 f; @@ -107,7 +180,7 @@ return u.u; } -double hs_fastHalfToDouble(uint16_t x) +double hs_fp_ieee_halfToDouble(uint16_t x) { union { _Float16 f; diff --git a/cbits/minmax.c b/cbits/minmax.c --- a/cbits/minmax.c +++ b/cbits/minmax.c @@ -13,28 +13,28 @@ float hs_minimumFloat(float x, float y) { float result; - asm("fmin %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmin %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); return result; } float hs_maximumFloat(float x, float y) { float result; - asm("fmax %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmax %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); return result; } double hs_minimumDouble(double x, double y) { double result; - asm("fmin %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmin %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); return result; } double hs_maximumDouble(double x, double y) { double result; - asm("fmax %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmax %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); return result; } @@ -50,9 +50,9 @@ // Therefore, we convert signaling NaNs to quiet ones before applying FMINNM. // x *= 1.0f; // y *= 1.0f; - asm("fmul %s0, %s0, %s1" : "+w"(x) : "w"(1.0f)); - asm("fmul %s0, %s0, %s1" : "+w"(y) : "w"(1.0f)); - asm("fminnm %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmul %s0, %s0, %s1" : "+w"(x) : "w"(1.0f)); + __asm__("fmul %s0, %s0, %s1" : "+w"(y) : "w"(1.0f)); + __asm__("fminnm %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); return result; } @@ -63,9 +63,9 @@ // Therefore, we convert signaling NaNs to quiet ones before applying FMAXNM. // x *= 1.0f; // y *= 1.0f; - asm("fmul %s0, %s0, %s1" : "+w"(x) : "w"(1.0f)); - asm("fmul %s0, %s0, %s1" : "+w"(y) : "w"(1.0f)); - asm("fmaxnm %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmul %s0, %s0, %s1" : "+w"(x) : "w"(1.0f)); + __asm__("fmul %s0, %s0, %s1" : "+w"(y) : "w"(1.0f)); + __asm__("fmaxnm %s0, %s1, %s2" : "=w"(result) : "w"(x), "w"(y)); return result; } @@ -76,9 +76,9 @@ // Therefore, we convert signaling NaNs to quiet ones before applying FMINNM. // x *= 1.0; // y *= 1.0; - asm("fmul %d0, %d0, %d1" : "+w"(x) : "w"(1.0)); - asm("fmul %d0, %d0, %d1" : "+w"(y) : "w"(1.0)); - asm("fminnm %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmul %d0, %d0, %d1" : "+w"(x) : "w"(1.0)); + __asm__("fmul %d0, %d0, %d1" : "+w"(y) : "w"(1.0)); + __asm__("fminnm %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); return result; } @@ -89,9 +89,69 @@ // Therefore, we convert signaling NaNs to quiet ones before applying FMAXNM. // x *= 1.0; // y *= 1.0; - asm("fmul %d0, %d0, %d1" : "+w"(x) : "w"(1.0)); - asm("fmul %d0, %d0, %d1" : "+w"(y) : "w"(1.0)); - asm("fmaxnm %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); + __asm__("fmul %d0, %d0, %d1" : "+w"(x) : "w"(1.0)); + __asm__("fmul %d0, %d0, %d1" : "+w"(y) : "w"(1.0)); + __asm__("fmaxnm %d0, %d1, %d2" : "=w"(result) : "w"(x), "w"(y)); + return result; +} + +#elif defined(__riscv) + +#if defined(__riscv_zfa) +float hs_minimumFloat(float x, float y) +{ + float result; + __asm__("fminm.s %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); + return result; +} + +float hs_maximumFloat(float x, float y) +{ + float result; + __asm__("fmaxm.s %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); + return result; +} + +double hs_minimumDouble(double x, double y) +{ + double result; + __asm__("fminm.d %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); + return result; +} + +double hs_maximumDouble(double x, double y) +{ + double result; + __asm__("fmaxm.d %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); + return result; +} +#endif + +float hs_minimumNumberFloat(float x, float y) +{ + float result; + __asm__("fmin.s %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); + return result; +} + +float hs_maximumNumberFloat(float x, float y) +{ + float result; + __asm__("fmax.s %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); + return result; +} + +double hs_minimumNumberDouble(double x, double y) +{ + double result; + __asm__("fmin.d %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); + return result; +} + +double hs_maximumNumberDouble(double x, double y) +{ + double result; + __asm__("fmax.d %0, %1, %2" : "=f"(result) : "f"(x), "f"(y)); return result; } diff --git a/cbits/roundeven.c b/cbits/roundeven.c --- a/cbits/roundeven.c +++ b/cbits/roundeven.c @@ -29,7 +29,7 @@ { float result; // a floating-exception can be generated - asm("frintn %s0, %s1" : "=w"(result) : "w"(x)); + __asm__("frintn %s0, %s1" : "=w"(result) : "w"(x)); return result; } @@ -37,7 +37,25 @@ { double result; // a floating-exception can be generated - asm("frintn %d0, %d1" : "=w"(result) : "w"(x)); + __asm__("frintn %d0, %d1" : "=w"(result) : "w"(x)); + return result; +} + +#elif defined(__riscv_zfa) + +float hs_roundevenFloat(float x) +{ + float result; + // a floating-exception can be generated for sNaN + __asm__("fround.s %0, %1, rne" : "=f"(result) : "f"(x)); + return result; +} + +double hs_roundevenDouble(double x) +{ + double result; + // a floating-exception can be generated for sNaN + __asm__("fround.d %0, %1, rne" : "=f"(result) : "f"(x)); return result; } diff --git a/fp-ieee.cabal b/fp-ieee.cabal --- a/fp-ieee.cabal +++ b/fp-ieee.cabal @@ -1,7 +1,7 @@ cabal-version: 2.2 name: fp-ieee -version: 0.1.0.6 +version: 0.1.0.7 synopsis: IEEE 754-2019 compliant operations description: Please see the README on GitHub at category: Numeric, Math @@ -9,7 +9,7 @@ bug-reports: https://github.com/minoki/haskell-floating-point/issues author: ARATA Mizuki maintainer: minorinoki@gmail.com -copyright: 2020-2025 ARATA Mizuki +copyright: 2020-2026 ARATA Mizuki license: BSD-3-Clause license-file: LICENSE build-type: Simple @@ -42,6 +42,18 @@ manual: True default: False +-- Zfhmin +flag rva22u64 + description: Enable RVA22U64 profile on RISC-V + manual: True + default: False + +-- Zfhmin, Zfa +flag rva23u64 + description: Enable RVA23U64 profile on RISC-V + manual: True + default: False + -- flag x87-long-double -- description: Support x87 "long double" via long-double package -- manual: True @@ -72,7 +84,7 @@ -- We use a post-GHC 8.6 language extension: NumericUnderscores -- cast{Word32,Word64}To{Float,Double}, cast{Float,Double}To{Word32,Word64} are since base-4.10.0.0 (GHC 8.2) -- Semigroup((<>)) is exported from Prelude since base-4.11.0.0 (GHC 8.4) - base >=4.12 && <4.23 + base >=4.12 && <4.24 if !flag(pure-hs) cpp-options: -DUSE_FFI if !flag(pure-hs) && os(windows) @@ -155,7 +167,7 @@ integer-gmp ==1.0.* if flag(ghc-bignum) && impl(ghc >= 9.0.0) build-depends: - ghc-bignum >=1.0 && <1.4 + ghc-bignum >=1.0 && <1.5 -- Fast roundeven: needs SSE4.1 on x86 if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1) cpp-options: -DHAS_FAST_ROUNDEVEN @@ -167,6 +179,11 @@ cpp-options: -DHAS_FAST_ROUNDEVEN c-sources: cbits/roundeven.c + if !flag(pure-hs) && arch(riscv64) && flag(rva23u64) + -- RISC-V Zfa + cpp-options: -DHAS_FAST_ROUNDEVEN + c-sources: + cbits/roundeven.c -- Fast FMA: needs FMA3 on x86 (FMA4 is not supported by this package) if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3) && impl(ghc < 9.8) cpp-options: -DHAS_FAST_FMA @@ -176,12 +193,12 @@ if (arch(i386) || arch(x86_64)) && flag(fma3) && impl(ghc >= 9.8) cpp-options: -DHAS_FMA_PRIM ghc-options: -mfma - -- Fast FMA: always available on AArch64 - if !flag(pure-hs) && arch(aarch64) && impl(ghc < 9.8) + -- Fast FMA: always available on AArch64 and RISC-V + if !flag(pure-hs) && (arch(aarch64) || arch(riscv64)) && impl(ghc < 9.8) cpp-options: -DHAS_FAST_FMA c-sources: cbits/fma.c - if !flag(pure-hs) && arch(aarch64) && impl(ghc >= 9.8) + if !flag(pure-hs) && (arch(aarch64) || arch(riscv64)) && impl(ghc >= 9.8) cpp-options: -DHAS_FMA_PRIM -- Enable use of libm's fma unless "pure-hs" is set; but not on Windows -- (mingw-w64's fma is not reliable) @@ -191,9 +208,16 @@ cpp-options: -DDONT_INLINE_FMA_PRIM -- Fast min/max: available on AArch64 if !flag(pure-hs) && arch(aarch64) - cpp-options: -DHAS_FAST_MINMAX + cpp-options: -DHAS_FAST_MINMAX -DHAS_FAST_MINMAXNUM c-sources: cbits/minmax.c + if !flag(pure-hs) && arch(riscv64) + cpp-options: -DHAS_FAST_MINMAXNUM + if flag(rva23u64) + -- Zfa + cpp-options: -DHAS_FAST_MINMAX + c-sources: + cbits/minmax.c if !flag(pure-hs) && flag(half) && arch(x86_64) && flag(f16c) cpp-options: -DHAS_FAST_HALF_CONVERSION cc-options: -mf16c @@ -203,10 +227,22 @@ cpp-options: -DHAS_FAST_HALF_CONVERSION c-sources: cbits/half.c - if !flag(pure-hs) && (arch(aarch64) || arch(x86_64)) + if !flag(pure-hs) && flag(half) && arch(riscv64) && (flag(rva22u64) || flag(rva23u64)) + -- Zfhmin + cpp-options: -DHAS_FAST_HALF_CONVERSION + c-sources: + cbits/half.c + if !flag(pure-hs) && (arch(aarch64) || arch(x86_64) || arch(riscv64)) cpp-options: -DHAS_FAST_CANONICALIZE c-sources: cbits/canonicalize.c + if arch(riscv64) + if flag(rva23u64) + -- Zfhmin, Zfa + cc-options: -march=rva23u64 + elif flag(rva22u64) + -- Zfhmin + cc-options: -march=rva22u64 default-language: Haskell2010 test-suite fp-ieee-doctests @@ -214,8 +250,8 @@ type: exitcode-stdio-1.0 main-is: doctests.hs build-depends: - doctest >=0.22.2 && <0.25 - , QuickCheck >=2.14.3 && <2.18 + doctest >=0.22.2 && <0.26 + , QuickCheck >=2.14.3 && <2.19 default-language: Haskell2010 if impl(ghc >= 9.14) buildable: False @@ -227,11 +263,13 @@ other-modules: AugmentedArithSpec ClassificationSpec + ConversionSpec FMASpec IntegerInternalsSpec MinMaxSpec NaNSpec NextFloatSpec + RemainderSpec RoundingSpec RoundToIntegralSpec TwoSumSpec @@ -240,7 +278,7 @@ test ghc-options: -threaded -rtsopts -with-rtsopts=-N -fno-ignore-asserts build-depends: - QuickCheck >=2.14.3 && <2.18 + QuickCheck >=2.14.3 && <2.19 , fp-ieee , hspec ^>=2.11.7 , hspec-core ^>=2.11.7 @@ -277,4 +315,37 @@ CPP HexFloatLiterals NumericUnderscores + if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1) + cpp-options: -DHAS_FAST_ROUNDEVEN + if !flag(pure-hs) && arch(aarch64) + cpp-options: -DHAS_FAST_ROUNDEVEN + if !flag(pure-hs) && arch(riscv64) && flag(rva23u64) + cpp-options: -DHAS_FAST_ROUNDEVEN + if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3) && impl(ghc < 9.8) + cpp-options: -DHAS_FAST_FMA + if (arch(i386) || arch(x86_64)) && flag(fma3) && impl(ghc >= 9.8) + cpp-options: -DHAS_FMA_PRIM + ghc-options: -mfma + if !flag(pure-hs) && (arch(aarch64) || arch(riscv64)) && impl(ghc < 9.8) + cpp-options: -DHAS_FAST_FMA + if !flag(pure-hs) && (arch(aarch64) || arch(riscv64)) && impl(ghc >= 9.8) + cpp-options: -DHAS_FMA_PRIM + if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && !os(windows) + cpp-options: -DUSE_C99_FMA + if (arch(i386) || arch(x86_64)) && os(windows) + cpp-options: -DDONT_INLINE_FMA_PRIM + if !flag(pure-hs) && arch(aarch64) + cpp-options: -DHAS_FAST_MINMAX -DHAS_FAST_MINMAXNUM + if !flag(pure-hs) && arch(riscv64) + cpp-options: -DHAS_FAST_MINMAXNUM + if flag(rva23u64) + cpp-options: -DHAS_FAST_MINMAX + if !flag(pure-hs) && flag(half) && arch(x86_64) && flag(f16c) + cpp-options: -DHAS_FAST_HALF_CONVERSION + if !flag(pure-hs) && flag(half) && arch(aarch64) + cpp-options: -DHAS_FAST_HALF_CONVERSION + if !flag(pure-hs) && flag(half) && arch(riscv64) && (flag(rva22u64) || flag(rva23u64)) + cpp-options: -DHAS_FAST_HALF_CONVERSION + if !flag(pure-hs) && (arch(aarch64) || arch(x86_64) || arch(riscv64)) + cpp-options: -DHAS_FAST_CANONICALIZE default-language: Haskell2010 diff --git a/src/Numeric/Floating/IEEE/Internal/Classify.hs b/src/Numeric/Floating/IEEE/Internal/Classify.hs --- a/src/Numeric/Floating/IEEE/Internal/Classify.hs +++ b/src/Numeric/Floating/IEEE/Internal/Classify.hs @@ -46,6 +46,9 @@ -- IEEE 754 @isZero@ operation. isZero :: RealFloat a => a -> Bool isZero x = x == 0 +{-# NOINLINE [1] isZero #-} +{-# SPECIALIZE isZero :: Float -> Bool #-} +{-# SPECIALIZE isZero :: Double -> Bool #-} -- | -- Returns @True@ if the argument is negative (including negative zero). @@ -57,6 +60,9 @@ -- IEEE 754 @isSignMinus@ operation. isSignMinus :: RealFloat a => a -> Bool isSignMinus x = x < 0 || isNegativeZero x +{-# NOINLINE [1] isSignMinus #-} +{-# SPECIALIZE isSignMinus :: Float -> Bool #-} +{-# SPECIALIZE isSignMinus :: Double -> Bool #-} -- | -- Comparison with IEEE 754 @totalOrder@ predicate. diff --git a/src/Numeric/Floating/IEEE/Internal/FMA.hs b/src/Numeric/Floating/IEEE/Internal/FMA.hs --- a/src/Numeric/Floating/IEEE/Internal/FMA.hs +++ b/src/Numeric/Floating/IEEE/Internal/FMA.hs @@ -269,7 +269,9 @@ result0 = v1 + w !_ = assert (result0 == fromRational (toRational x + toRational y + toRational c'')) () result = scaleFloat e result0 - !_ = assert (result == fromRational (toRational a * toRational b + toRational c) || isDenormalized result) () + -- result is normal: fst (floatRange _) <= exponent result <= snd (floatRange _) + resultMaybeInexact = exponent result0 + e < fst (floatRange result0) + !_ = assert (result == fromRational (toRational a * toRational b + toRational c) || resultMaybeInexact) () in if result0 == 0 then -- We need to handle the sign of zero if c == 0 && a /= 0 && b /= 0 then @@ -277,7 +279,7 @@ else a * b + c -- -0 if both a * b and c are -0 else - if isDenormalized result then + if resultMaybeInexact then -- The rounding in 'scaleFloat e result0' may yield an incorrect result. -- Take the slow path. case toRational a * toRational b + toRational c of diff --git a/src/Numeric/Floating/IEEE/Internal/Float128.hs b/src/Numeric/Floating/IEEE/Internal/Float128.hs --- a/src/Numeric/Floating/IEEE/Internal/Float128.hs +++ b/src/Numeric/Floating/IEEE/Internal/Float128.hs @@ -55,7 +55,7 @@ nextUpF128 :: Float128 -> Float128 nextUpF128 x = case float128ToWord64Pair x of - (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000 + (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_0000 , (hi, lo) /= (0xffff_0000_0000_0000, 0) -> x + x -- NaN or positive infinity -> itself (0x8000_0000_0000_0000, 0x0000_0000_0000_0000) -> minPositive -- -0 -> min positive (hi, lo) | testBit hi 63 -> -- negative @@ -68,7 +68,7 @@ nextDownF128 :: Float128 -> Float128 nextDownF128 x = case float128ToWord64Pair x of - (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000 + (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_0000 , (hi, lo) /= (0x7fff_0000_0000_0000, 0) -> x + x -- NaN or negative infinity -> itself (0x0000_0000_0000_0000, 0x0000_0000_0000_0000) -> - minPositive -- +0 -> max negative (hi, lo) | testBit hi 63 -> -- negative @@ -81,7 +81,7 @@ nextTowardZeroF128 :: Float128 -> Float128 nextTowardZeroF128 x = case float128ToWord64Pair x of - (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000 + (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_0000 , (lo, hi .&. 0x0000_ffff_ffff_ffff) /= (0, 0) -> x + x -- NaN -> itself (0x8000_0000_0000_0000, 0x0000_0000_0000_0000) -> x -- -0 -> itself (0x0000_0000_0000_0000, 0x0000_0000_0000_0000) -> x -- +0 -> itself @@ -97,7 +97,7 @@ isFiniteF128 :: Float128 -> Bool isFiniteF128 x = case float128ToWord64Pair x of (hi, _) -> let hi' = hi .&. 0x7fff_0000_0000_0000 - in hi' /= 0 && hi' /= 0x7fff_0000_0000_0000 + in hi' /= 0x7fff_0000_0000_0000 classifyF128DiscardingSignalingNaNs :: Float128 -> Class classifyF128DiscardingSignalingNaNs x = diff --git a/src/Numeric/Floating/IEEE/Internal/Half.hs b/src/Numeric/Floating/IEEE/Internal/Half.hs --- a/src/Numeric/Floating/IEEE/Internal/Half.hs +++ b/src/Numeric/Floating/IEEE/Internal/Half.hs @@ -53,7 +53,7 @@ nextTowardZeroHalf x = case castHalfToWord16 x of w | w .&. 0x7c00 == 0x7c00 - , w /= 0x7fff -> x + x -- NaN -> itself + , w .&. 0x03ff /= 0 -> x + x -- NaN -> itself 0x8000 -> x -- -0 -> itself 0x0000 -> x -- +0 -> itself w -> castWord16ToHalf (w - 1) -- positive / negative @@ -165,7 +165,7 @@ {-# SPECIALIZE fromPositiveIntegerR :: RoundingStrategy f => Bool -> Integer -> f Half #-} {-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTiesToEven Half #-} {-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTiesToAway Half #-} -{-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTowardPositive Half　#-} +{-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTowardPositive Half #-} {-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTowardNegative Half #-} {-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTowardZero Half #-} @@ -198,13 +198,13 @@ #if defined(HAS_FAST_HALF_CONVERSION) -foreign import ccall unsafe "hs_fastHalfToFloat" +foreign import ccall unsafe "hs_fp_ieee_halfToFloat" c_fastHalfToFloat :: Word16 -> Float -foreign import ccall unsafe "hs_fastHalfToDouble" +foreign import ccall unsafe "hs_fp_ieee_halfToDouble" c_fastHalfToDouble :: Word16 -> Double -foreign import ccall unsafe "hs_fastFloatToHalf" +foreign import ccall unsafe "hs_fp_ieee_floatToHalf" c_fastFloatToHalf :: Float -> Word16 -foreign import ccall unsafe "hs_fastDoubleToHalf" +foreign import ccall unsafe "hs_fp_ieee_doubleToHalf" c_fastDoubleToHalf :: Double -> Word16 halfToFloat = coerce c_fastHalfToFloat diff --git a/src/Numeric/Floating/IEEE/Internal/MinMax.hs b/src/Numeric/Floating/IEEE/Internal/MinMax.hs --- a/src/Numeric/Floating/IEEE/Internal/MinMax.hs +++ b/src/Numeric/Floating/IEEE/Internal/MinMax.hs @@ -81,47 +81,59 @@ minimumFloat :: Float -> Float -> Float foreign import ccall unsafe "hs_maximumFloat" maximumFloat :: Float -> Float -> Float -foreign import ccall unsafe "hs_minimumNumberFloat" - minimumNumberFloat :: Float -> Float -> Float -foreign import ccall unsafe "hs_maximumNumberFloat" - maximumNumberFloat :: Float -> Float -> Float foreign import ccall unsafe "hs_minimumDouble" minimumDouble :: Double -> Double -> Double foreign import ccall unsafe "hs_maximumDouble" maximumDouble :: Double -> Double -> Double + +{-# RULES +"minimum'/Float" minimum' = minimumFloat +"maximum'/Float" maximum' = maximumFloat +"minimum'/Double" minimum' = minimumDouble +"maximum'/Double" maximum' = maximumDouble + #-} + +#else + +minimumFloat :: Float -> Float -> Float +maximumFloat :: Float -> Float -> Float +minimumDouble :: Double -> Double -> Double +maximumDouble :: Double -> Double -> Double + +minimumFloat = minimum' +minimumDouble = minimum' +maximumFloat = maximum' +maximumDouble = maximum' + +#endif + +#if defined(HAS_FAST_MINMAXNUM) + +foreign import ccall unsafe "hs_minimumNumberFloat" + minimumNumberFloat :: Float -> Float -> Float +foreign import ccall unsafe "hs_maximumNumberFloat" + maximumNumberFloat :: Float -> Float -> Float foreign import ccall unsafe "hs_minimumNumberDouble" minimumNumberDouble :: Double -> Double -> Double foreign import ccall unsafe "hs_maximumNumberDouble" maximumNumberDouble :: Double -> Double -> Double {-# RULES -"minimum'/Float" minimum' = minimumFloat -"maximum'/Float" maximum' = maximumFloat "minimumNumber/Float" minimumNumber = minimumNumberFloat "maximumNumber/Float" maximumNumber = maximumNumberFloat -"minimum'/Double" minimum' = minimumDouble -"maximum'/Double" maximum' = maximumDouble "minimumNumber/Double" minimumNumber = minimumNumberDouble "maximumNumber/Double" maximumNumber = maximumNumberDouble #-} #else -minimumFloat :: Float -> Float -> Float -maximumFloat :: Float -> Float -> Float minimumNumberFloat :: Float -> Float -> Float maximumNumberFloat :: Float -> Float -> Float -minimumDouble :: Double -> Double -> Double -maximumDouble :: Double -> Double -> Double minimumNumberDouble :: Double -> Double -> Double maximumNumberDouble :: Double -> Double -> Double -minimumFloat = minimum' -minimumDouble = minimum' minimumNumberFloat = minimumNumber minimumNumberDouble = minimumNumber -maximumFloat = maximum' -maximumDouble = maximum' maximumNumberFloat = maximumNumber maximumNumberDouble = maximumNumber diff --git a/src/Numeric/Floating/IEEE/Internal/Remainder.hs b/src/Numeric/Floating/IEEE/Internal/Remainder.hs --- a/src/Numeric/Floating/IEEE/Internal/Remainder.hs +++ b/src/Numeric/Floating/IEEE/Internal/Remainder.hs @@ -17,7 +17,7 @@ | y == 0 || isInfinite y || isNaN y || not (isFinite x) = (x - x) / y * y -- return a NaN | otherwise = let n = round (toRational x / toRational y) r = fromRational (toRational x - toRational y * fromInteger n) - in r -- if r == 0, the sign of r is the same as x + in if r == 0 then 0 * x else r -- if r == 0, return the signed zero {-# NOINLINE [1] remainder #-} #if defined(USE_FFI) diff --git a/test/ClassificationSpec.hs b/test/ClassificationSpec.hs --- a/test/ClassificationSpec.hs +++ b/test/ClassificationSpec.hs @@ -29,6 +29,7 @@ , counterexample "isSignMinus" $ isSignMinus x === (c `elem` [NegativeInfinity, NegativeNormal, NegativeSubnormal, NegativeZero]) -- isSignMinus doesn't handle negative NaNs ] where c = classify x +{-# INLINABLE prop_classify #-} {-# SPECIALIZE prop_classify :: Proxy Float -> Float -> Property #-} {-# SPECIALIZE prop_classify :: Proxy Double -> Double -> Property #-} diff --git a/test/ConversionSpec.hs b/test/ConversionSpec.hs new file mode 100644 --- /dev/null +++ b/test/ConversionSpec.hs @@ -0,0 +1,30 @@ +module ConversionSpec where +import Data.Proxy +import Numeric +import Numeric.Floating.IEEE +import Test.Hspec +import Test.Hspec.QuickCheck +import Test.QuickCheck +import Util + +default () + +prop_conversion :: (RealFloat a, RealFloat b, Show a, Show b) => Proxy a -> Proxy b -> a -> Property +prop_conversion _ proxyB x = + let y = realFloatToFrac x `asProxyTypeOf` proxyB + y' | isInfinite x = if y > 0 then 1 / 0 else -(1 / 0) + | isNaN x = 0 / 0 + | isNegativeZero x = -0 + | otherwise = fromRat (toRational x) + in y `sameFloatP` y' +{-# INLINABLE prop_conversion #-} + +{-# NOINLINE spec #-} +spec :: Spec +spec = modifyMaxSuccess (* 1000) $ do + let proxyFloat :: Proxy Float + proxyFloat = Proxy + proxyDouble :: Proxy Double + proxyDouble = Proxy + prop "Float->Double" $ forAllFloats $ prop_conversion proxyFloat proxyDouble + prop "Double->Float" $ forAllFloats $ prop_conversion proxyDouble proxyFloat diff --git a/test/FMASpec.hs b/test/FMASpec.hs --- a/test/FMASpec.hs +++ b/test/FMASpec.hs @@ -69,6 +69,7 @@ , (0x1.ffffffc000000p512, 0x1.0000002p511, -0x1p-1074, 0x1.fffffffffffffp1023) -- 0x1.ffffffc000000p512 * 0x1.0000002p511 == 0x1.fffffffffffff8p1023 (in Rational) , (-0x1.032ede48bbb28p-1022, 0x1.3cbc999ae14a8p-1, -0x1p-1074, -0x1.40accc50d63d2p-1023) , (0x1.ca903c622e5a6p-1022, 0x1.414a00c886a44p-1, 0x1.f1a8235fd56fep-1022, 0x1.88b4ec63db4f5p-1021) + , (0x1.ffffffffffffdp-511, 0x1.0000000000001p-512, 0.0, 0x1.ffffffffffffep-1023) ] casesForFloat :: [(Float, Float, Float, Float)] @@ -80,6 +81,7 @@ , (0x1.83bd78p4, -0x1.cp118, -0x1.344108p-2, -0x1.5345cap123) , (0x1p-149, 0x1.88dd0cp-1, 0x1.081ffp-127, 0x1.081ff4p-127) , (0x1.d1a9dp-126, 0x1.594da4p-1, 0x1.343de4p-126, 0x1.3725b6p-125) + , (0x1.fffffap-63, 0x1.000002p-64, 0.0, 0x1.fffffcp-127) ] testSpecialValues :: (RealFloat a, Show a) => String -> (a -> a -> a -> a) -> [(a, a, a, a)] -> Spec diff --git a/test/Float128Spec.hs b/test/Float128Spec.hs --- a/test/Float128Spec.hs +++ b/test/Float128Spec.hs @@ -1,5 +1,6 @@ {-# LANGUAGE CPP #-} {-# LANGUAGE HexFloatLiterals #-} +{-# LANGUAGE NumericUnderscores #-} {-# LANGUAGE ScopedTypeVariables #-} {-# OPTIONS_GHC -Wno-orphans #-} module Float128Spec where @@ -7,6 +8,7 @@ augmentedMultiplication_viaRational) import qualified AugmentedArithSpec import qualified ClassificationSpec +import qualified ConversionSpec import Control.Monad import Data.Function (on) import Data.Functor.Identity @@ -22,6 +24,7 @@ import Numeric.Floating.IEEE import Numeric.Floating.IEEE.Internal import Numeric.Floating.IEEE.NaN (setPayloadSignaling) +import qualified RemainderSpec import qualified RoundingSpec import qualified RoundToIntegralSpec import System.Random @@ -48,8 +51,14 @@ (x,g') = random g in (fromRational (toInteger x % 2^(16 :: Int)), g') -- TODO +fromRationalIsBuggy :: Spec -> Spec +fromRationalIsBuggy = mapSpecItem_ (allowFailure "Float128's fromRational may be incorrect") + +roundIsBuggy :: Spec -> Spec +roundIsBuggy = mapSpecItem_ (allowFailure "Float128's round may be incorrect") + spec :: Spec -spec = mapSpecItem_ (allowFailure "Float128's fromRational and round may be incorrect") $ do +spec = do let proxy :: Proxy Float128 proxy = Proxy prop "classify" $ forAllFloats $ ClassificationSpec.prop_classify proxy @@ -57,19 +66,21 @@ prop "totalOrder" $ forAllFloats2 $ ClassificationSpec.prop_totalOrder proxy prop "totalOrder (generic)" $ forAllFloats2 (ClassificationSpec.prop_totalOrder (Proxy :: Proxy (Identity Float128)) `on` Identity) prop "twoSum" $ forAllFloats2 $ TwoSumSpec.prop_twoSum proxy - prop "twoProduct" $ forAllFloats2 $ TwoSumSpec.prop_twoProduct proxy twoProduct - prop "twoProduct_generic" $ forAllFloats2 $ TwoSumSpec.prop_twoProduct proxy twoProduct_generic + fromRationalIsBuggy $ prop "twoProduct" $ forAllFloats2 $ TwoSumSpec.prop_twoProduct proxy twoProduct + fromRationalIsBuggy $ prop "twoProduct_generic" $ forAllFloats2 $ TwoSumSpec.prop_twoProduct proxy twoProduct_generic let casesForFloat128 :: [(Float128, Float128, Float128, Float128)] casesForFloat128 = [ (-0, 0, -0, -0) , (-0, -0, -0, 0) + , (0x1.ffff_ffff_ffff_ffff_ffff_ffff_fffdp-8191, 0x1.0000_0000_0000_0000_0000_0000_0001p-8192, 0.0, 0x1.ffff_ffff_ffff_ffff_ffff_ffff_ffffp-16383) -- TODO: Add more ] - FMASpec.checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForFloat128 - FMASpec.checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForFloat128 - FMASpec.checkFMA "fusedMultiplyAdd (via Rational)" fusedMultiplyAdd_viaRational casesForFloat128 + fromRationalIsBuggy $ FMASpec.checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForFloat128 + fromRationalIsBuggy $ FMASpec.checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForFloat128 + fromRationalIsBuggy $ FMASpec.checkFMA "fusedMultiplyAdd (via Rational)" fusedMultiplyAdd_viaRational casesForFloat128 prop "nextUp . nextDown == id (unless -inf)" $ forAllFloats $ NextFloatSpec.prop_nextUp_nextDown proxy prop "nextDown . nextUp == id (unless inf)" $ forAllFloats $ NextFloatSpec.prop_nextDown_nextUp proxy - prop "augmentedAddition/equality" $ forAllFloats2 $ \(x :: Float128) y -> + prop "nextTowardZero == (nextUp or nextDown, unless 0.0)" $ forAllFloats $ NextFloatSpec.prop_nextTowardZero proxy + fromRationalIsBuggy $ prop "augmentedAddition/equality" $ forAllFloats2 $ \(x :: Float128) y -> isFinite x && isFinite y ==> let (s,t) = augmentedAddition x y in isFinite s ==> isFinite t .&&. toRational s + toRational t === toRational x + toRational y @@ -77,11 +88,12 @@ augmentedAddition x y `sameFloatPairP` augmentedAddition_viaRational x y prop "augmentedMultiplication" $ forAllFloats2 $ \(x :: Float128) y -> augmentedMultiplication x y `sameFloatPairP` augmentedMultiplication_viaRational x y + fromRationalIsBuggy $ prop "remainder" $ forAllFloats2 $ RemainderSpec.prop_remainder proxy prop "fromIntegerR vs fromRationalR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_fromIntegerR_vs_fromRationalR proxy) prop "fromIntegerR vs encodeFloatR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_fromIntegerR_vs_encodeFloatR proxy) prop "fromRationalR vs encodeFloatR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_fromRationalR_vs_encodeFloatR proxy) - prop "fromRationalR vs fromRational" $ RoundingSpec.prop_fromRationalR_vs_fromRational proxy + fromRationalIsBuggy $ prop "fromRationalR vs fromRational" $ RoundingSpec.prop_fromRationalR_vs_fromRational proxy prop "scaleFloatR vs fromRationalR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_scaleFloatR_vs_fromRationalR proxy) prop "scaleFloatR vs encodeFloatR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_scaleFloatR_vs_encodeFloatR proxy) prop "result of fromIntegerR" $ \x -> RoundingSpec.prop_order proxy (fromIntegerR x) @@ -89,8 +101,14 @@ prop "result of encodeFloatR" $ \m k -> RoundingSpec.prop_order proxy (encodeFloatR m k) prop "addToOdd" $ forAllFloats2 $ RoundingSpec.prop_addToOdd proxy - prop "roundToIntegral" $ RoundToIntegralSpec.prop_roundToIntegral proxy - RoundToIntegralSpec.checkCases proxy + roundIsBuggy $ prop "roundToIntegral" $ RoundToIntegralSpec.prop_roundToIntegral proxy + roundIsBuggy $ RoundToIntegralSpec.checkCases proxy + + modifyMaxSuccess (* 1000) $ do + prop "Float->Float128" $ forAllFloats $ ConversionSpec.prop_conversion (Proxy :: Proxy Float) proxy + prop "Double->Float128" $ forAllFloats $ ConversionSpec.prop_conversion (Proxy :: Proxy Double) proxy + prop "Float128->Float" $ forAllFloats $ ConversionSpec.prop_conversion proxy (Proxy :: Proxy Float) + prop "Float128->Double" $ forAllFloats $ ConversionSpec.prop_conversion proxy (Proxy :: Proxy Double) prop "copySign" $ forAllFloats2 $ NaNSpec.prop_copySign proxy prop "isSignMinus" $ forAllFloats $ NaNSpec.prop_isSignMinus proxy diff --git a/test/HalfSpec.hs b/test/HalfSpec.hs --- a/test/HalfSpec.hs +++ b/test/HalfSpec.hs @@ -7,6 +7,7 @@ augmentedMultiplication_viaRational) import qualified AugmentedArithSpec import qualified ClassificationSpec +import qualified ConversionSpec import Control.Monad import Data.Function (on) import Data.Functor.Identity @@ -18,10 +19,12 @@ import qualified FMASpec import qualified NaNSpec import qualified NextFloatSpec +import Numeric import Numeric.Floating.IEEE import Numeric.Floating.IEEE.Internal import Numeric.Floating.IEEE.NaN (setPayloadSignaling) import Numeric.Half +import qualified RemainderSpec import qualified RoundingSpec import qualified RoundToIntegralSpec import System.Random @@ -59,8 +62,16 @@ isInfiniteIsKnownToBeBuggy = True #endif +-- https://github.com/ekmett/half/issues/41 +fromRationalIsBuggy :: Spec -> Spec +fromRationalIsBuggy = mapSpecItem_ (allowFailure "Half's fromRational may be incorrect") + +-- https://github.com/ekmett/half/issues/43 +decodeFloatIsBuggy :: Spec -> Spec +decodeFloatIsBuggy = mapSpecItem_ (allowFailure "Half's decodeFloat may be incorrect") + spec :: Spec -spec = mapSpecItem_ (allowFailure "Half's fromRational may be incorrect") $ do +spec = do let proxy :: Proxy Half proxy = Proxy prop "classify" $ forAllFloats $ isInfiniteWorkaround $ ClassificationSpec.prop_classify proxy @@ -73,21 +84,24 @@ let casesForHalf :: [(Half, Half, Half, Half)] casesForHalf = [ (-0, 0, -0, -0) , (-0, -0, -0, 0) + , (0x1.ff4p-7, 0x1.004p-8, 0, 0x1.ff8p-15) -- TODO: Add more ] - FMASpec.checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForHalf - FMASpec.checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForHalf - FMASpec.checkFMA "fusedMultiplyAdd (via Rational)" fusedMultiplyAdd_viaRational casesForHalf + decodeFloatIsBuggy $ FMASpec.checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForHalf + decodeFloatIsBuggy $ FMASpec.checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForHalf + decodeFloatIsBuggy $ FMASpec.checkFMA "fusedMultiplyAdd (via Rational)" fusedMultiplyAdd_viaRational casesForHalf prop "nextUp . nextDown == id (unless -inf)" $ forAllFloats $ NextFloatSpec.prop_nextUp_nextDown proxy prop "nextDown . nextUp == id (unless inf)" $ forAllFloats $ NextFloatSpec.prop_nextDown_nextUp proxy + prop "nextTowardZero == (nextUp or nextDown, unless 0.0)" $ forAllFloats $ NextFloatSpec.prop_nextTowardZero proxy prop "augmentedAddition/equality" $ forAllFloats2 $ \(x :: Half) y -> isFinite x && isFinite y ==> let (s,t) = augmentedAddition x y in isFinite s ==> isFinite t .&&. toRational s + toRational t === toRational x + toRational y prop "augmentedAddition" $ forAllFloats2 $ \(x :: Half) y -> augmentedAddition x y `sameFloatPairP` augmentedAddition_viaRational x y - prop "augmentedMultiplication" $ forAllFloats2 $ \(x :: Half) y -> + decodeFloatIsBuggy $ prop "augmentedMultiplication" $ forAllFloats2 $ \(x :: Half) y -> augmentedMultiplication x y `sameFloatPairP` augmentedMultiplication_viaRational x y + prop "remainder" $ forAllFloats2 $ RemainderSpec.prop_remainder proxy prop "fromIntegerR vs fromRationalR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_fromIntegerR_vs_fromRationalR proxy) prop "fromIntegerR vs encodeFloatR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_fromIntegerR_vs_encodeFloatR proxy) @@ -97,11 +111,26 @@ prop "scaleFloatR vs encodeFloatR" $ RoundingSpec.eachStrategy (RoundingSpec.prop_scaleFloatR_vs_encodeFloatR proxy) prop "result of fromIntegerR" $ \x -> RoundingSpec.prop_order proxy (fromIntegerR x) prop "result of fromRationalR" $ \x -> RoundingSpec.prop_order proxy (fromRationalR x) - prop "result of encodeFloatR" $ \m k -> RoundingSpec.prop_order proxy (encodeFloatR m k) - prop "addToOdd" $ forAllFloats2 $ RoundingSpec.prop_addToOdd proxy + decodeFloatIsBuggy $ prop "result of encodeFloatR" $ \m k -> RoundingSpec.prop_order proxy (encodeFloatR m k) + decodeFloatIsBuggy $ prop "addToOdd" $ forAllFloats2 $ RoundingSpec.prop_addToOdd proxy prop "roundToIntegral" $ RoundToIntegralSpec.prop_roundToIntegral proxy RoundToIntegralSpec.checkCases proxy + + modifyMaxSuccess (* 1000) $ do + prop "Half->Float" $ forAllFloats $ ConversionSpec.prop_conversion proxy (Proxy :: Proxy Float) + prop "Half->Double" $ forAllFloats $ ConversionSpec.prop_conversion proxy (Proxy :: Proxy Double) + prop "Float->Half" $ forAllFloats $ ConversionSpec.prop_conversion (Proxy :: Proxy Float) proxy + fromRationalIsBuggy $ prop "Double->Half" $ forAllFloats $ ConversionSpec.prop_conversion (Proxy :: Proxy Double) proxy + let casesFromDouble :: [Double] + casesFromDouble = [ 0x1.001ffffcp+0 + , 0x1.8p-25 + , 0x1.0000000000001p-25 + , -0x1.fffffffffffffp-25 + ] + forM_ casesFromDouble $ \x -> do + let label = showString "Double->Half (" . showHFloat x $ ")" + fromRationalIsBuggy $ it label $ ConversionSpec.prop_conversion (Proxy :: Proxy Double) proxy x prop "copySign" $ forAllFloats2 $ NaNSpec.prop_copySign proxy prop "isSignMinus" $ forAllFloats $ NaNSpec.prop_isSignMinus proxy diff --git a/test/NextFloatSpec.hs b/test/NextFloatSpec.hs --- a/test/NextFloatSpec.hs +++ b/test/NextFloatSpec.hs @@ -47,12 +47,21 @@ prop_nextUp_nextDown _ x = x /= (-1/0) ==> let x' = nextUp (nextDown x) in x' `sameFloatP` x .||. (isPositiveZero x .&&. isNegativeZero x') +{-# INLINABLE prop_nextUp_nextDown #-} prop_nextDown_nextUp :: (RealFloat a, Show a) => Proxy a -> a -> Property prop_nextDown_nextUp _ x = x /= (1/0) ==> let x' = nextDown (nextUp x) in x' `sameFloatP` x .||. (isNegativeZero x .&&. isPositiveZero x') +{-# INLINABLE prop_nextDown_nextUp #-} +prop_nextTowardZero :: (RealFloat a, Show a) => Proxy a -> a -> Property +prop_nextTowardZero _ x + | x == 0.0 = nextTowardZero x `sameFloatP` x + | isSignMinus x = nextTowardZero x `sameFloatP` nextUp x + | otherwise = nextTowardZero x `sameFloatP` nextDown x +{-# INLINABLE prop_nextTowardZero #-} + {-# NOINLINE spec #-} spec :: Spec spec = do @@ -66,6 +75,7 @@ #endif prop "nextUp . nextDown == id (unless -inf)" $ forAllFloats $ prop_nextUp_nextDown proxy prop "nextDown . nextUp == id (unless inf)" $ forAllFloats $ prop_nextDown_nextUp proxy + prop "nextTowardZero == (nextUp or nextDown, unless 0.0)" $ forAllFloats $ prop_nextTowardZero proxy describe "Float" $ do let proxy :: Proxy Float @@ -77,3 +87,4 @@ #endif prop "nextUp . nextDown == id (unless -inf)" $ forAllFloats $ prop_nextUp_nextDown proxy prop "nextDown . nextUp == id (unless inf)" $ forAllFloats $ prop_nextDown_nextUp proxy + prop "nextTowardZero == (nextUp or nextDown, unless 0.0)" $ forAllFloats $ prop_nextTowardZero proxy diff --git a/test/RemainderSpec.hs b/test/RemainderSpec.hs new file mode 100644 --- /dev/null +++ b/test/RemainderSpec.hs @@ -0,0 +1,30 @@ +module RemainderSpec where +import Data.Proxy +import Numeric.Floating.IEEE +import Test.Hspec +import Test.Hspec.QuickCheck +import Test.QuickCheck hiding (classify) +import Util + +prop_remainder :: (RealFloat a, Show a) => Proxy a -> a -> a -> Property +prop_remainder _ x y + | isFinite x && isFinite y && y /= 0 = + let n = round (toRational x / toRational y) + r = toRational x - toRational y * fromInteger n + r' = if r == 0 then x * 0 else fromRational r + in remainder x y `sameFloatP` r' + | isFinite x && isInfinite y = remainder x y `sameFloatP` x + | otherwise = isNaN (remainder x y) === True +{-# INLINABLE prop_remainder #-} + +{-# NOINLINE spec #-} +spec :: Spec +spec = do + describe "Double" $ do + let proxy :: Proxy Double + proxy = Proxy + prop "remainder" $ forAllFloats2 (prop_remainder proxy) + describe "Float" $ do + let proxy :: Proxy Float + proxy = Proxy + prop "remainder" $ forAllFloats2 (prop_remainder proxy) diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,11 +1,13 @@ {-# LANGUAGE CPP #-} import qualified AugmentedArithSpec import qualified ClassificationSpec +import qualified ConversionSpec import qualified FMASpec import qualified IntegerInternalsSpec import qualified MinMaxSpec import qualified NaNSpec import qualified NextFloatSpec +import qualified RemainderSpec import qualified RoundingSpec import qualified RoundToIntegralSpec import System.Environment (getArgs, withArgs) @@ -37,11 +39,13 @@ main :: IO () main = hspec $ do describe "Classification" ClassificationSpec.spec + describe "Conversion" ConversionSpec.spec describe "TwoSum" TwoSumSpec.spec describe "FMA" FMASpec.spec describe "IntegerInternals" IntegerInternalsSpec.spec describe "NextFloat" NextFloatSpec.spec describe "AugmentedArith" AugmentedArithSpec.spec + describe "Remainder" RemainderSpec.spec describe "Rounding" RoundingSpec.spec describe "RoundToIntegral" RoundToIntegralSpec.spec describe "NaN" NaNSpec.spec