fp-ieee 0.1.0.2 → 0.1.0.3
raw patch · 14 files changed
+201/−31 lines, 14 filesdep ~QuickCheckdep ~basedep ~doctestPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: QuickCheck, base, doctest, ghc-bignum, hspec, hspec-core, random, tasty-bench
API changes (from Hackage documentation)
Files
- ChangeLog.md +7/−0
- README.md +1/−1
- benchmark/Benchmark.hs +1/−1
- cbits/fma.c +23/−0
- doctests.hs +1/−0
- fp-ieee.cabal +26/−16
- src/Numeric/Floating/IEEE/Internal/Augmented.hs +4/−1
- src/Numeric/Floating/IEEE/Internal/FMA.hs +85/−4
- src/Numeric/Floating/IEEE/Internal/GenericArith.hs +4/−0
- src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs +2/−2
- src/Numeric/Floating/IEEE/Internal/NextFloat.hs +7/−6
- src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs +25/−0
- test/AugmentedArithSpec.hs +2/−0
- test/FMASpec.hs +13/−0
ChangeLog.md view
@@ -1,5 +1,12 @@ # Changelog for fp-ieee +## Version 0.1.0.3 (2023-11-18)++* Allow ghc-bignum 1.3.+* Fix `IntegerInternals.roundingMode#`.+* Fix assertion in `augmentedMultiplication`.+* Use FMA primitives on GHC 9.8.+ ## Version 0.1.0.2 (2021-11-30) * Allow ghc-bignum 1.2.
README.md view
@@ -6,7 +6,7 @@ * correctly-rounding versions of `fromInteger`. * `realFloatToFrac`, which correctly handles signed zeros, infinities, and NaNs (unlike `realToFrac`). -Some operations (e.g. `fusedMultiplyAdd`) can make use of the native instruction in the architecture.+Some operations (e.g. `fusedMultiplyAdd`) can make use of the native instruction in the architecture via C FFI, or GHC 9.8's FMA primitives. For non-native targets, "Pure Haskell" mode is supported via a package flag.
benchmark/Benchmark.hs view
@@ -6,7 +6,6 @@ import Data.Coerce import Data.Functor.Identity import Data.Word-import Gauge.Main import GHC.Float (isDoubleFinite, isFloatFinite) import Numeric import Numeric.Floating.IEEE@@ -18,6 +17,7 @@ #if defined(USE_FLOAT128) import Numeric.Float128 (Float128) #endif+import Test.Tasty.Bench foreign import ccall unsafe "nextafter" c_nextafter_double :: Double -> Double -> Double
cbits/fma.c view
@@ -1,5 +1,26 @@ #include <math.h> +#if defined(__GNUC__) && defined(__FMA__)++// Make sure FMA instruction is used even if optimizations are disabled.++double hs_fusedMultiplyAddDouble(double a, double b, double c)+{+ // vfmadd132sd %xmm1, %xmm2, %xmm0+ // %xmm0 <- %xmm0 * %xmm1 + %xmm2+ __asm__("vfmadd132sd %1, %2, %0" : "+x"(a) : "x"(b), "x"(c));+ return a;+}+float hs_fusedMultiplyAddFloat(float a, float b, float c)+{+ // vfmadd132ss %xmm1, %xmm2, %xmm0+ // %xmm0 <- %xmm0 * %xmm1 + %xmm2+ __asm__("vfmadd132ss %1, %2, %0" : "+x"(a) : "x"(b), "x"(c));+ return a;+}++#else+ #if !defined(FP_FAST_FMA) #error "The compiler should define FP_FAST_FMA" #endif@@ -15,3 +36,5 @@ { return fmaf(a, b, c); }++#endif
doctests.hs view
@@ -2,6 +2,7 @@ main :: IO () main = doctest [ "-isrc"+ , "-fobject-code" -- for GHC 8.6 , "src/Numeric/Floating/IEEE/Internal/Base.hs" , "src/Numeric/Floating/IEEE/Internal/Classify.hs" , "src/Numeric/Floating/IEEE/Internal/FMA.hs"
fp-ieee.cabal view
@@ -1,7 +1,7 @@ cabal-version: 2.2 name: fp-ieee-version: 0.1.0.2+version: 0.1.0.3 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-2021 ARATA Mizuki+copyright: 2020-2023 ARATA Mizuki license: BSD-3-Clause license-file: LICENSE build-type: Simple@@ -33,7 +33,7 @@ default: False flag fma3- description: Use FMA3 instructions on x86+ description: Use FMA3 instructions on x86. On GHC 9.8 or later, this flag enables use of FMA primitives. manual: True default: False @@ -72,9 +72,12 @@ -- 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 && <5+ base >=4.12 && <4.20 if !flag(pure-hs) cpp-options: -DUSE_FFI+ if !flag(pure-hs) && os(windows)+ -- mingw-w64 is not reliable+ cpp-options: -DSOME_LIBC_FUNCTIONS_MIGHT_BE_BUGGY if flag(float128) -- Support Float128 cpp-options: -DUSE_FLOAT128@@ -141,7 +144,7 @@ integer-gmp ==1.0.* if flag(ghc-bignum) && impl(ghc >= 9.0.0) build-depends:- ghc-bignum >=1.0 && <1.3+ ghc-bignum >=1.0 && <1.4 -- Fast roundeven: needs SSE4.1 on x86 if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1) cpp-options: -DHAS_FAST_ROUNDEVEN@@ -154,20 +157,27 @@ 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)+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3) && impl(ghc < 9.8) cpp-options: -DHAS_FAST_FMA cc-options: -mfma c-sources: cbits/fma.c+ 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)+ if !flag(pure-hs) && arch(aarch64) && impl(ghc < 9.8) cpp-options: -DHAS_FAST_FMA c-sources: cbits/fma.c+ if !flag(pure-hs) && arch(aarch64) && 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) 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 -- Fast min/max: available on AArch64 if !flag(pure-hs) && arch(aarch64) cpp-options: -DHAS_FAST_MINMAX@@ -193,9 +203,11 @@ type: exitcode-stdio-1.0 main-is: doctests.hs build-depends:- doctest >=0.8- , QuickCheck+ doctest ^>=0.22.2+ , QuickCheck ^>=2.14.3 default-language: Haskell2010+ if impl(ghc >= 9.8)+ buildable: False test-suite fp-ieee-test import: deps, options@@ -217,12 +229,12 @@ test ghc-options: -threaded -rtsopts -with-rtsopts=-N -fno-ignore-asserts build-depends:- QuickCheck+ QuickCheck ^>=2.14.3 , fp-ieee- , hspec- , hspec-core+ , hspec ^>=2.11.7+ , hspec-core ^>=2.11.7 , integer-logarithms- , random+ , random ^>=1.2.1.1 if flag(half) other-modules: HalfSpec@@ -239,7 +251,5 @@ benchmark build-depends: fp-ieee- , tasty-bench- mixins:- tasty-bench (Test.Tasty.Bench as Gauge, Test.Tasty.Bench as Gauge.Main)+ , tasty-bench ^>=0.3.5 default-language: Haskell2010
src/Numeric/Floating/IEEE/Internal/Augmented.hs view
@@ -68,7 +68,10 @@ in if exy + exponent u1 >= expMin then -- The result is exact let ulpTowardZero = u1 - nextTowardZero u1- !_ = assert (2 * abs u2 <= abs ulpTowardZero) ()+ !_ = assert (case u1 of+ 0.5 -> - ulpTowardZero <= 2 * u2 && u2 <= ulpTowardZero+ -0.5 -> ulpTowardZero <= u2 && 2 * u2 <= - ulpTowardZero+ _ -> 2 * abs u2 <= abs ulpTowardZero) () (v1, v2) = if (-2) * u2 == ulpTowardZero then (u1 - ulpTowardZero, ulpTowardZero + u2) else
src/Numeric/Floating/IEEE/Internal/FMA.hs view
@@ -1,6 +1,10 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE CPP #-} {-# LANGUAGE NoImplicitPrelude #-}+#if defined(HAS_FMA_PRIM)+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE UnboxedTuples #-}+#endif module Numeric.Floating.IEEE.Internal.FMA ( isMantissaEven , twoSum@@ -25,11 +29,15 @@ isFloatBinary32, (^!)) import Numeric.Floating.IEEE.Internal.Classify (isFinite) import Numeric.Floating.IEEE.Internal.NextFloat (nextDown, nextUp)+#if defined(HAS_FMA_PRIM)+import GHC.Exts+#endif default () -- $setup -- >>> :set -XScopedTypeVariables+-- >>> import Numeric.Floating.IEEE.Internal.FMA -- Assumption: input is finite isMantissaEven :: RealFloat a => a -> Bool@@ -151,8 +159,45 @@ twoProductFloat :: Float -> Float -> (Float, Float) twoProductDouble :: Double -> Double -> (Double, Double) -#if defined(HAS_FAST_FMA)+#if defined(HAS_FMA_PRIM) && 0+-- Disabled for now: https://gitlab.haskell.org/ghc/ghc/-/issues/24160 +twoProductFloat# :: Float# -> Float# -> (# Float#, Float# #)+twoProductFloat# x y = let !r = x `timesFloat#` y+ !s = fmsubFloat# x y r+ in (# r, s #)++twoProductDouble# :: Double# -> Double# -> (# Double#, Double# #)+twoProductDouble# x y = let !r = x *## y+ !s = fmsubDouble# x y r+ in (# r, s #)++#if defined(DONT_INLINE_FMA_PRIM)+{-# NOINLINE twoProductFloat# #-}+{-# NOINLINE twoProductDouble# #-}+#else+{-# INLINE twoProductFloat# #-}+{-# INLINE twoProductDouble# #-}+#endif++twoProductFloat (F# x) (F# y) = case twoProductFloat# x y of+ (# r, s #) -> (F# r, F# s)++twoProductDouble (D# x) (D# y) = case twoProductDouble# x y of+ (# r, s #) -> (D# r, D# s)++{-# INLINE twoProductFloat #-}+{-# INLINE twoProductDouble #-}++{-# RULES+"twoProduct/Float" twoProduct = twoProductFloat+"twoProduct/Double" twoProduct = twoProductDouble+"twoProduct_nonscaling/Float" twoProduct_nonscaling = twoProductFloat+"twoProduct_nonscaling/Double" twoProduct_nonscaling = twoProductDouble+ #-}++#elif defined(HAS_FAST_FMA) || defined(HAS_FMA_PRIM)+ twoProductFloat x y = let !r = x * y !s = fusedMultiplyAddFloat x y (-r) in (r, s)@@ -255,10 +300,46 @@ | isFinite a && isFinite b = c + c -- a * b is finite, but c is Infinity or NaN | otherwise = a * b + c where- !True = isFloatBinary32 || error "fusedMultiplyAdd/Float: Float must be IEEE binary32"- !True = isDoubleBinary64 || error "fusedMultiplyAdd/Float: Double must be IEEE binary64"+ !() = if isFloatBinary32 then () else error "fusedMultiplyAdd/Float: Float must be IEEE binary32"+ !() = if isDoubleBinary64 then () else error "fusedMultiplyAdd/Float: Double must be IEEE binary64" -#if defined(HAS_FAST_FMA)+#if defined(HAS_FMA_PRIM)++#if defined(DONT_INLINE_FMA_PRIM)++fusedMultiplyAddFloat# :: Float# -> Float# -> Float# -> Float#+fusedMultiplyAddFloat# x y z = fmaddFloat# x y z+{-# NOINLINE fusedMultiplyAddFloat# #-}++fusedMultiplyAddDouble# :: Double# -> Double# -> Double# -> Double#+fusedMultiplyAddDouble# x y z = fmaddDouble# x y z+{-# NOINLINE fusedMultiplyAddDouble# #-}++fusedMultiplyAddFloat :: Float -> Float -> Float -> Float+fusedMultiplyAddFloat (F# x) (F# y) (F# z) = F# (fusedMultiplyAddFloat# x y z)++fusedMultiplyAddDouble :: Double -> Double -> Double -> Double+fusedMultiplyAddDouble (D# x) (D# y) (D# z) = D# (fusedMultiplyAddDouble# x y z)++#else++fusedMultiplyAddFloat :: Float -> Float -> Float -> Float+fusedMultiplyAddFloat (F# x) (F# y) (F# z) = F# (fmaddFloat# x y z)++fusedMultiplyAddDouble :: Double -> Double -> Double -> Double+fusedMultiplyAddDouble (D# x) (D# y) (D# z) = D# (fmaddDouble# x y z)++#endif++{-# INLINE fusedMultiplyAddFloat #-}+{-# INLINE fusedMultiplyAddDouble #-}++{-# RULES+"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat+"fusedMultiplyAdd/Double" fusedMultiplyAdd = fusedMultiplyAddDouble+ #-}++#elif defined(HAS_FAST_FMA) foreign import ccall unsafe "hs_fusedMultiplyAddFloat" fusedMultiplyAddFloat :: Float -> Float -> Float -> Float
src/Numeric/Floating/IEEE/Internal/GenericArith.hs view
@@ -8,6 +8,10 @@ default () +-- $setup+-- >>> :m + Data.Proxy+-- >>> import Numeric.Floating.IEEE.Internal.GenericArith+ infixl 6 `genericAdd`, `genericSub` infixl 7 `genericMul`, `genericDiv`
src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs view
@@ -46,7 +46,7 @@ #endif -- $setup--- >>> :m + Data.Int Test.QuickCheck+-- >>> :m + Data.Int Data.Bits Test.QuickCheck -- >>> :{ -- -- Workaround for https://github.com/sol/doctest/issues/160: -- import Numeric.Floating.IEEE.Internal.IntegerInternals@@ -208,7 +208,7 @@ <> loop s where loop 0# = EQ- loop i = case GHC.Num.BigNat.bigNatIndex# bn i of+ loop i = case GHC.Num.BigNat.bigNatIndex# bn (i -# 1#) of 0## -> loop (i -# 1#) _ -> GT
src/Numeric/Floating/IEEE/Internal/NextFloat.hs view
@@ -12,6 +12,7 @@ -- $setup -- >>> :set -XHexFloatLiterals -XNumericUnderscores+-- >>> import Numeric.Floating.IEEE.Internal.NextFloat -- | -- Returns the smallest value that is larger than the argument.@@ -185,7 +186,7 @@ w | testBit w 31 -> castWord32ToFloat (w - 1) -- negative | otherwise -> castWord32ToFloat (w + 1) -- positive where- !True = isFloatBinary32 || error "Numeric.Floating.Extra assumes Float is IEEE binary32"+ !() = if isFloatBinary32 then () else error "Numeric.Floating.IEEE assumes Float is IEEE binary32" -- | -- prop> nextUpDouble 1 == 0x1.0000_0000_0000_1p0@@ -203,7 +204,7 @@ w | testBit w 63 -> castWord64ToDouble (w - 1) -- negative | otherwise -> castWord64ToDouble (w + 1) -- positive where- !True = isDoubleBinary64 || error "Numeric.Floating.Extra assumes Double is IEEE binary64"+ !() = if isDoubleBinary64 then () else error "Numeric.Floating.IEEE assumes Double is IEEE binary64" -- | -- prop> nextDownFloat 1 == 0x1.fffffep-1@@ -221,7 +222,7 @@ w | testBit w 31 -> castWord32ToFloat (w + 1) -- negative | otherwise -> castWord32ToFloat (w - 1) -- positive where- !True = isFloatBinary32 || error "Numeric.Floating.Extra assumes Float is IEEE binary32"+ !() = if isFloatBinary32 then () else error "Numeric.Floating.IEEE assumes Float is IEEE binary32" -- | -- prop> nextDownDouble 1 == 0x1.ffff_ffff_ffff_fp-1@@ -239,7 +240,7 @@ w | testBit w 63 -> castWord64ToDouble (w + 1) -- negative | otherwise -> castWord64ToDouble (w - 1) -- positive where- !True = isDoubleBinary64 || error "Numeric.Floating.Extra assumes Double is IEEE binary64"+ !() = if isDoubleBinary64 then () else error "Numeric.Floating.IEEE assumes Double is IEEE binary64" -- | -- prop> nextTowardZeroFloat 1 == 0x1.fffffep-1@@ -258,7 +259,7 @@ 0x0000_0000 -> x -- +0 -> itself w -> castWord32ToFloat (w - 1) -- positive / negative where- !True = isFloatBinary32 || error "Numeric.Floating.Extra assumes Float is IEEE binary32"+ !() = if isFloatBinary32 then () else error "Numeric.Floating.IEEE assumes Float is IEEE binary32" -- | -- prop> nextTowardZeroDouble 1 == 0x1.ffff_ffff_ffff_fp-1@@ -277,4 +278,4 @@ 0x0000_0000_0000_0000 -> x -- +0 -> itself w -> castWord64ToDouble (w - 1) -- positive / negative where- !True = isDoubleBinary64 || error "Numeric.Floating.Extra assumes Double is IEEE binary64"+ !() = if isDoubleBinary64 then () else error "Numeric.Floating.IEEE assumes Double is IEEE binary64"
src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs view
@@ -13,12 +13,16 @@ , floor ) where import MyPrelude+#if defined(USE_FFI) && defined(SOME_LIBC_FUNCTIONS_MIGHT_BE_BUGGY)+import Numeric.Floating.IEEE.Internal.Conversion+#endif default () -- $setup -- >>> :set -XScopedTypeVariables -- >>> import Numeric.Floating.IEEE.Internal.Classify (isFinite)+-- >>> import Numeric.Floating.IEEE.Internal.RoundToIntegral -- | -- @'round'' x@ returns the nearest integral value to @x@; the even integer if @x@ is equidistant between two integers.@@ -151,8 +155,28 @@ foreign import ccall unsafe "trunc" c_truncDouble :: Double -> Double +#if defined(SOME_LIBC_FUNCTIONS_MIGHT_BE_BUGGY) {-# RULES "roundAway'/Float"+ roundAway' = c_roundFloat . canonicalizeFloat+"roundAway'/Double"+ roundAway' = c_roundDouble . canonicalizeDouble+"truncate'/Float"+ truncate' = c_truncFloat+"truncate'/Double"+ truncate' = c_truncDouble+"ceiling'/Float"+ ceiling' = c_ceilFloat+"ceiling'/Double"+ ceiling' = c_ceilDouble+"floor'/Float"+ floor' = c_floorFloat+"floor'/Double"+ floor' = c_floorDouble+ #-}+#else+{-# RULES+"roundAway'/Float" roundAway' = c_roundFloat "roundAway'/Double" roundAway' = c_roundDouble@@ -169,6 +193,7 @@ "floor'/Double" floor' = c_floorDouble #-}+#endif {- from base foreign import ccall unsafe "rintFloat"
test/AugmentedArithSpec.hs view
@@ -104,6 +104,8 @@ , (0x1.5433bcp-126, -0x1.69a04p-1, -0x1.e091e8p-127, -0x0p+0) , (0x1.c7363p-128, -0x1.c5d164p-1, -0x1.937b98p-128, -0x0p+0) , (-0x1.a31946p0, -0x1p-127, 0x1.a31944p-127, 0x0p+0)+ , (0x1.cp1, 0x1.24924ap-1, 0x1p1, 0x1.8p-24)+ , (0x1.cp1, -0x1.24924ap-1, -0x1p1, -0x1.8p-24) ] prop "augmentedMultiplication" $ testAugmented augmentedMultiplication cases prop "augmentedMultiplication_viaRational" $ testAugmented augmentedMultiplication_viaRational cases
test/FMASpec.hs view
@@ -6,6 +6,7 @@ import Data.Bits import Data.Coerce import Data.Functor.Identity+import GHC.Exts (inline) import Numeric import Numeric.Floating.IEEE import Numeric.Floating.IEEE.Internal@@ -24,6 +25,14 @@ #endif +fusedMultiplyAddInlineFloat :: Float -> Float -> Float -> Float+fusedMultiplyAddInlineFloat x y z = inline fusedMultiplyAdd x y z+{-# NOINLINE fusedMultiplyAddInlineFloat #-}++fusedMultiplyAddInlineDouble :: Double -> Double -> Double -> Double+fusedMultiplyAddInlineDouble x y z = inline fusedMultiplyAdd x y z+{-# NOINLINE fusedMultiplyAddInlineDouble #-}+ fusedMultiplyAdd_generic :: RealFloat a => a -> a -> a -> a fusedMultiplyAdd_generic x y z = runIdentity (fusedMultiplyAdd (Identity x) (Identity y) (Identity z)) @@ -88,11 +97,15 @@ spec = modifyMaxSuccess (* 100) $ do describe "Double" $ do checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForDouble+ checkFMA "fusedMultiplyAdd (monomorphic)" fusedMultiplyAddDouble casesForDouble+ checkFMA "fusedMultiplyAdd (inline)" fusedMultiplyAddInlineDouble casesForDouble checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForDouble checkFMA "fusedMultiplyAdd (via Rational)" fusedMultiplyAdd_viaRational casesForDouble checkFMA "fusedMultiplyAdd (via Integer)" fusedMultiplyAdd_viaInteger casesForDouble describe "Float" $ do checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForFloat+ checkFMA "fusedMultiplyAdd (monomorphic)" fusedMultiplyAddFloat casesForFloat+ checkFMA "fusedMultiplyAdd (inline)" fusedMultiplyAddInlineFloat casesForFloat checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForFloat checkFMA "fusedMultiplyAdd (via Rational)" fusedMultiplyAdd_viaRational casesForFloat checkFMA "fusedMultiplyAdd (via Integer)" fusedMultiplyAdd_viaInteger casesForFloat