fp-ieee 0.1.0.6 → 0.1.0.7
raw patch · 21 files changed
+508/−114 lines, 21 filesdep ~QuickCheckdep ~basedep ~doctestPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: QuickCheck, base, doctest, ghc-bignum
API changes (from Hackage documentation)
Files
- ChangeLog.md +14/−0
- LICENSE +1/−1
- cbits/canonicalize.c +17/−4
- cbits/half.c +106/−33
- cbits/minmax.c +76/−16
- cbits/roundeven.c +20/−2
- fp-ieee.cabal +83/−12
- src/Numeric/Floating/IEEE/Internal/Classify.hs +6/−0
- src/Numeric/Floating/IEEE/Internal/FMA.hs +4/−2
- src/Numeric/Floating/IEEE/Internal/Float128.hs +4/−4
- src/Numeric/Floating/IEEE/Internal/Half.hs +6/−6
- src/Numeric/Floating/IEEE/Internal/MinMax.hs +28/−16
- src/Numeric/Floating/IEEE/Internal/Remainder.hs +1/−1
- test/ClassificationSpec.hs +1/−0
- test/ConversionSpec.hs +30/−0
- test/FMASpec.hs +2/−0
- test/Float128Spec.hs +28/−10
- test/HalfSpec.hs +36/−7
- test/NextFloatSpec.hs +11/−0
- test/RemainderSpec.hs +30/−0
- test/Spec.hs +4/−0
ChangeLog.md view
@@ -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.
LICENSE view
@@ -1,4 +1,4 @@-Copyright ARATA Mizuki (c) 2020-2024+Copyright ARATA Mizuki (c) 2020-2026 All rights reserved.
cbits/canonicalize.c view
@@ -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; }
cbits/half.c view
@@ -1,11 +1,50 @@ #include <stdint.h> // uint16_t+#include <string.h> #include <math.h> +/*+ * 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 <x86intrin.h> -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;
cbits/minmax.c view
@@ -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; }
cbits/roundeven.c view
@@ -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; }
fp-ieee.cabal view
@@ -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 <https://github.com/minoki/haskell-floating-point/tree/master/fp-ieee#readme> 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
src/Numeric/Floating/IEEE/Internal/Classify.hs view
@@ -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.
src/Numeric/Floating/IEEE/Internal/FMA.hs view
@@ -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
src/Numeric/Floating/IEEE/Internal/Float128.hs view
@@ -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 =
src/Numeric/Floating/IEEE/Internal/Half.hs view
@@ -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
src/Numeric/Floating/IEEE/Internal/MinMax.hs view
@@ -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
src/Numeric/Floating/IEEE/Internal/Remainder.hs view
@@ -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)
test/ClassificationSpec.hs view
@@ -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 #-}
+ test/ConversionSpec.hs view
@@ -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
test/FMASpec.hs view
@@ -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
test/Float128Spec.hs view
@@ -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
test/HalfSpec.hs view
@@ -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
test/NextFloatSpec.hs view
@@ -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
+ test/RemainderSpec.hs view
@@ -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)
test/Spec.hs view
@@ -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