fp-ieee (empty) → 0.1.0
raw patch · 50 files changed
+6444/−0 lines, 50 filesdep +QuickCheckdep +basedep +decimal-arithmeticsetup-changed
Dependencies added: QuickCheck, base, decimal-arithmetic, doctest, float128, fp-ieee, gauge, ghc-bignum, half, hspec, hspec-core, integer-gmp, integer-logarithms, random
Files
- ChangeLog.md +5/−0
- LICENSE +30/−0
- README.md +13/−0
- Setup.hs +2/−0
- benchmark/Benchmark.hs +494/−0
- cbits/canonicalize.c +64/−0
- cbits/fma.c +17/−0
- cbits/half.c +120/−0
- cbits/minmax.c +102/−0
- cbits/roundeven.c +48/−0
- decimal-test/NextFloatSpec.hs +42/−0
- decimal-test/Spec.hs +27/−0
- doctests.hs +13/−0
- fp-ieee.cabal +435/−0
- src/GHC/Float/Compat.hs +26/−0
- src/MyPrelude.hs +14/−0
- src/Numeric/Floating/IEEE.hs +250/−0
- src/Numeric/Floating/IEEE/Internal.hs +23/−0
- src/Numeric/Floating/IEEE/Internal/Augmented.hs +127/−0
- src/Numeric/Floating/IEEE/Internal/Base.hs +136/−0
- src/Numeric/Floating/IEEE/Internal/Classify.hs +157/−0
- src/Numeric/Floating/IEEE/Internal/Conversion.hs +68/−0
- src/Numeric/Floating/IEEE/Internal/FMA.hs +301/−0
- src/Numeric/Floating/IEEE/Internal/Float128.hs +232/−0
- src/Numeric/Floating/IEEE/Internal/GenericArith.hs +103/−0
- src/Numeric/Floating/IEEE/Internal/Half.hs +254/−0
- src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs +259/−0
- src/Numeric/Floating/IEEE/Internal/MinMax.hs +128/−0
- src/Numeric/Floating/IEEE/Internal/NaN.hs +219/−0
- src/Numeric/Floating/IEEE/Internal/NextFloat.hs +280/−0
- src/Numeric/Floating/IEEE/Internal/Remainder.hs +35/−0
- src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs +193/−0
- src/Numeric/Floating/IEEE/Internal/Rounding.hs +7/−0
- src/Numeric/Floating/IEEE/Internal/Rounding/Common.hs +165/−0
- src/Numeric/Floating/IEEE/Internal/Rounding/Encode.hs +204/−0
- src/Numeric/Floating/IEEE/Internal/Rounding/Integral.hs +316/−0
- src/Numeric/Floating/IEEE/Internal/Rounding/Rational.hs +150/−0
- src/Numeric/Floating/IEEE/NaN.hs +15/−0
- test/AugmentedArithSpec.hs +111/−0
- test/ClassificationSpec.hs +63/−0
- test/FMASpec.hs +107/−0
- test/Float128Spec.hs +108/−0
- test/HalfSpec.hs +123/−0
- test/IntegerInternalsSpec.hs +53/−0
- test/MinMaxSpec.hs +134/−0
- test/NaNSpec.hs +166/−0
- test/RoundToIntegralSpec.hs +171/−0
- test/RoundingSpec.hs +241/−0
- test/Spec.hs +54/−0
- test/TwoSumSpec.hs +39/−0
+ ChangeLog.md view
@@ -0,0 +1,5 @@+# Changelog for fp-ieee++## Version 0.1.0 (2020-12-27)++Initial release.
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright ARATA Mizuki (c) 2020++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of ARATA Mizuki nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,13 @@+# fp-ieee: IEEE 754 operations for floating-point types++This library provides IEEE 754-compliant operations, including++* `fusedMultiplyAdd`.+* 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.++For non-native targets, "Pure Haskell" mode is supported via a package flag.++Most operations require only `RealFloat` constraint, but `RealFloatNaN` is needed by some operations that access the sign and payload of NaNs.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ benchmark/Benchmark.hs view
@@ -0,0 +1,494 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE HexFloatLiterals #-}+{-# LANGUAGE NumericUnderscores #-}+import Data.Bits+import Data.Coerce+import Data.Functor.Identity+import Data.Word+import Gauge.Main+import GHC.Float (isDoubleFinite, isFloatFinite)+import Numeric+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+#if defined(USE_HALF)+import Numeric.Half hiding (isZero)+import qualified Numeric.Half+#endif+#if defined(USE_FLOAT128)+import Numeric.Float128 (Float128)+#endif++foreign import ccall unsafe "nextafter"+ c_nextafter_double :: Double -> Double -> Double+foreign import ccall unsafe "nextafterf"+ c_nextafter_float :: Float -> Float -> Float+foreign import ccall unsafe "fma"+ c_fma_double :: Double -> Double -> Double -> Double+foreign import ccall unsafe "fmaf"+ c_fma_float :: Float -> Float -> Float -> Float++class Fractional a => CFloat a where+ c_nextafter :: a -> a -> a+ c_fma :: a -> a -> a -> a++instance CFloat Double where+ c_nextafter = c_nextafter_double+ c_fma = c_fma_double++instance CFloat Float where+ c_nextafter = c_nextafter_float+ c_fma = c_fma_float++c_nextUp, c_nextDown :: (RealFloat a, CFloat a) => a -> a+c_nextUp x = c_nextafter x (1/0)+c_nextDown x = c_nextafter x (-1/0)++twoProduct_generic :: RealFloat a => a -> a -> (a, a)+twoProduct_generic x y = coerce (twoProduct (Identity x) (Identity y))++fusedMultiplyAdd_generic :: RealFloat a => a -> a -> a -> a+fusedMultiplyAdd_generic x y z = runIdentity (fusedMultiplyAdd (Identity x) (Identity y) (Identity z))++fusedMultiplyAdd_viaInteger :: RealFloat a => a -> a -> a -> a+fusedMultiplyAdd_viaInteger x y z+ | isFinite x && isFinite y && isFinite z =+ let (mx,ex) = decodeFloat x -- x == mx * b^ex, mx==0 || b^(d-1) <= abs mx < b^d+ (my,ey) = decodeFloat y -- y == my * b^ey, my==0 || b^(d-1) <= abs my < b^d+ (mz,ez) = decodeFloat z -- z == mz * b^ez, mz==0 || b^(d-1) <= abs mz < b^d+ exy = ex + ey+ ee = min ez exy+ !2 = floatRadix x+ in case mx * my `shiftL` (exy - ee) + mz `shiftL` (ez - ee) of+ 0 -> x * y + z+ m -> roundTiesToEven (encodeFloatR m ee)+ | isFinite x && isFinite y = z + z -- x * y is finite, but z is Infinity or NaN+ | otherwise = x * y + z -- either x or y is Infinity or NaN++fusedMultiplyAdd_viaRational :: RealFloat a => a -> a -> a -> a+fusedMultiplyAdd_viaRational x y z+ | isFinite x && isFinite y && isFinite z =+ case toRational x * toRational y + toRational z of+ 0 -> x * y + z+ r -> fromRational r+ | isFinite x && isFinite y = z + z -- x * is finite, but z is Infinity or NaN+ | otherwise = x * y + z -- either x or y is Infinity or NaN++main :: IO ()+main = defaultMain+ [ bgroup "FMA"+ [ let arg = (1.0, 2.0, 3.0) :: (Double, Double, Double)+ in bgroup "Double"+ [ bench "C" $ nf (\(x,y,z) -> c_fma x y z) arg+ , bench "Haskell (default)" $ nf (\(x,y,z) -> fusedMultiplyAdd x y z) arg+ , bench "Haskell (generic)" $ nf (\(x,y,z) -> fusedMultiplyAdd_generic x y z) arg+ , bench "Haskell (via Rational)" $ nf (\(x,y,z) -> fusedMultiplyAdd_viaRational x y z) arg+ , bench "Haskell (via Integer)" $ nf (\(x,y,z) -> fusedMultiplyAdd_viaInteger x y z) arg+ , bench "non-fused" $ nf (\(x,y,z) -> x * y + z) arg+ ]+ , let arg = (1.0, 2.0, 3.0) :: (Float, Float, Float)+ in bgroup "Float"+ [ bench "C" $ nf (\(x,y,z) -> c_fma x y z) arg+ , bench "Haskell (default)" $ nf (\(x,y,z) -> fusedMultiplyAdd x y z) arg+ , bench "Haskell (generic)" $ nf (\(x,y,z) -> fusedMultiplyAdd_generic x y z) arg+ , bench "Haskell (via Rational)" $ nf (\(x,y,z) -> fusedMultiplyAdd_viaRational x y z) arg+ , bench "Haskell (via Integer)" $ nf (\(x,y,z) -> fusedMultiplyAdd_viaInteger x y z) arg+ , bench "Haskell (via Double)" $ nf (\(x,y,z) -> fusedMultiplyAddFloat_viaDouble x y z) arg+ , bench "non-fused" $ nf (\(x,y,z) -> x * y + z) arg+ ]+ ]+ , bgroup "isNormal"+ [ let arg = pi :: Double+ in bgroup "Double"+ [ bench "default" $ nf isNormal arg+ , bench "generic" $ nf (isNormal . Identity) arg+ ]+ , let arg = pi :: Float+ in bgroup "Float"+ [ bench "default" $ nf isNormal arg+ , bench "generic" $ nf (isNormal . Identity) arg+ ]+ ]+ , bgroup "isFinite"+ [ let arg = pi :: Double+ in bgroup "Double"+ [ bench "default" $ nf isFinite arg+ , bench "generic" $ nf (isFinite . Identity) arg+ , bench "GHC.Float.isDoubleFinite" $ nf isDoubleFinite arg+ ]+ , let arg = pi :: Float+ in bgroup "Float"+ [ bench "default" $ nf isFinite arg+ , bench "generic" $ nf (isFinite . Identity) arg+ , bench "GHC.Float.isFloatFinite" $ nf isFloatFinite arg+ ]+ ]+ , bgroup "twoProduct"+ [ let arg :: (Double, Double)+ arg = (1.3 * 2^500, pi / 2^500)+ in bgroup "Double"+ [ bench "Haskell (default)" $ nf (uncurry twoProduct) arg+ , bench "Haskell (generic)" $ nf (uncurry twoProduct_generic) arg+ , bench "Haskell (nonscaling)" $ nf (uncurry twoProduct_nonscaling) arg+#if defined(HAS_FAST_FMA)+ , bench "FMA" $ nf (uncurry twoProductDouble) arg+#endif+ ]+ , let arg :: (Float, Float)+ arg = (1.3 * 2^50, pi / 2^50)+ in bgroup "Float"+ [ bench "Haskell (default)" $ nf (uncurry twoProduct) arg+ , bench "Haskell (generic)" $ nf (uncurry twoProduct_generic) arg+ , bench "Haskell (nonscaling)" $ nf (uncurry twoProduct_nonscaling) arg+ , bench "Haskell (via Double)" $ nf (uncurry twoProductFloat_viaDouble) arg+#if defined(HAS_FAST_FMA)+ , bench "FMA" $ nf (uncurry twoProductFloat) arg+#endif+ ]+ ]+ , bgroup "fromInteger"+ [ let x = 418237418 * 2^80 + 4811 * 2^32 + 1412+ in bgroup "large"+ [ bgroup "Double"+ [ bench "stock" $ nf (fromInteger :: Integer -> Double) x+ , bench "fromIntegerTiesToEven" $ nf (fromIntegerTiesToEven :: Integer -> Double) x+ , bench "fromIntegerTiesToAway" $ nf (fromIntegerTiesToAway :: Integer -> Double) x+ , bench "fromIntegerTowardPositive" $ nf (fromIntegerTowardPositive :: Integer -> Double) x+ , bench "fromIntegerTowardNegative" $ nf (fromIntegerTowardNegative :: Integer -> Double) x+ , bench "fromIntegerTowardZero" $ nf (fromIntegerTowardZero :: Integer -> Double) x+ ]+ , bgroup "Float"+ [ bench "stock" $ nf (fromInteger :: Integer -> Float) x+ , bench "fromIntegerTiesToEven" $ nf (fromIntegerTiesToEven :: Integer -> Float) x+ , bench "fromIntegerTiesToAway" $ nf (fromIntegerTiesToAway :: Integer -> Float) x+ , bench "fromIntegerTowardPositive" $ nf (fromIntegerTowardPositive :: Integer -> Float) x+ , bench "fromIntegerTowardNegative" $ nf (fromIntegerTowardNegative :: Integer -> Float) x+ , bench "fromIntegerTowardZero" $ nf (fromIntegerTowardZero :: Integer -> Float) x+ ]+ ]+ , let x = 3 * 2^19 + 4811 * 2^7 + 1412+ in bgroup "small"+ [ bgroup "Double"+ [ bench "stock" $ nf (fromInteger :: Integer -> Double) x+ , bench "fromIntegerTiesToEven" $ nf (fromIntegerTiesToEven :: Integer -> Double) x+ , bench "fromIntegerTiesToAway" $ nf (fromIntegerTiesToAway :: Integer -> Double) x+ , bench "fromIntegerTowardPositive" $ nf (fromIntegerTowardPositive :: Integer -> Double) x+ , bench "fromIntegerTowardNegative" $ nf (fromIntegerTowardNegative :: Integer -> Double) x+ , bench "fromIntegerTowardZero" $ nf (fromIntegerTowardZero :: Integer -> Double) x+ ]+ , bgroup "Float"+ [ bench "stock" $ nf (fromInteger :: Integer -> Float) x+ , bench "fromIntegerTiesToEven" $ nf (fromIntegerTiesToEven :: Integer -> Float) x+ , bench "fromIntegerTiesToAway" $ nf (fromIntegerTiesToAway :: Integer -> Float) x+ , bench "fromIntegerTowardPositive" $ nf (fromIntegerTowardPositive :: Integer -> Float) x+ , bench "fromIntegerTowardNegative" $ nf (fromIntegerTowardNegative :: Integer -> Float) x+ , bench "fromIntegerTowardZero" $ nf (fromIntegerTowardZero :: Integer -> Float) x+ ]+ ]+ ]+ , bgroup "fromIntegral"+ [ bgroup "Word64"+ [ let x = 0xdead_beef_1234_7777 :: Word64+ in bgroup "large"+ [ bgroup "Double"+ [ bench "stock" $ nf (fromIntegral :: Word64 -> Double) x+ , bench "fromIntegralTiesToEven" $ nf (fromIntegralTiesToEven :: Word64 -> Double) x+ , bench "fromIntegralTiesToAway" $ nf (fromIntegralTiesToAway :: Word64 -> Double) x+ , bench "fromIntegralTowardPositive" $ nf (fromIntegralTowardPositive :: Word64 -> Double) x+ , bench "fromIntegralTowardNegative" $ nf (fromIntegralTowardNegative :: Word64 -> Double) x+ , bench "fromIntegralTowardZero" $ nf (fromIntegralTowardZero :: Word64 -> Double) x+ ]+ , bgroup "Float"+ [ bench "stock" $ nf (fromIntegral :: Word64 -> Float) x+ , bench "fromIntegralTiesToEven" $ nf (fromIntegralTiesToEven :: Word64 -> Float) x+ , bench "fromIntegralTiesToAway" $ nf (fromIntegralTiesToAway :: Word64 -> Float) x+ , bench "fromIntegralTowardPositive" $ nf (fromIntegralTowardPositive :: Word64 -> Float) x+ , bench "fromIntegralTowardNegative" $ nf (fromIntegralTowardNegative :: Word64 -> Float) x+ , bench "fromIntegralTowardZero" $ nf (fromIntegralTowardZero :: Word64 -> Float) x+ ]+ ]+ , let x = 0x14_7777 :: Word64+ in bgroup "small"+ [ bgroup "Double"+ [ bench "stock" $ nf (fromIntegral :: Word64 -> Double) x+ , bench "fromIntegralTiesToEven" $ nf (fromIntegralTiesToEven :: Word64 -> Double) x+ , bench "fromIntegralTiesToAway" $ nf (fromIntegralTiesToAway :: Word64 -> Double) x+ , bench "fromIntegralTowardPositive" $ nf (fromIntegralTowardPositive :: Word64 -> Double) x+ , bench "fromIntegralTowardNegative" $ nf (fromIntegralTowardNegative :: Word64 -> Double) x+ , bench "fromIntegralTowardZero" $ nf (fromIntegralTowardZero :: Word64 -> Double) x+ ]+ , bgroup "Float"+ [ bench "stock" $ nf (fromIntegral :: Word64 -> Float) x+ , bench "fromIntegralTiesToEven" $ nf (fromIntegralTiesToEven :: Word64 -> Float) x+ , bench "fromIntegralTiesToAway" $ nf (fromIntegralTiesToAway :: Word64 -> Float) x+ , bench "fromIntegralTowardPositive" $ nf (fromIntegralTowardPositive :: Word64 -> Float) x+ , bench "fromIntegralTowardNegative" $ nf (fromIntegralTowardNegative :: Word64 -> Float) x+ , bench "fromIntegralTowardZero" $ nf (fromIntegralTowardZero :: Word64 -> Float) x+ ]+ ]+ ]+ ]+ , bgroup "fromRational"+ [ let x = (418237418 * 2^80 + 4811 * 2^32 + 1412) / (2234321954 * 2^75 + 2345234566) :: Rational+ in bgroup "large/large"+ [ bgroup "Double"+ [ bench "stock" $ nf (fromRational :: Rational -> Double) x+ , bench "fromRationalTiesToEven" $ nf (fromRationalTiesToEven :: Rational -> Double) x+ , bench "fromRationalTiesToAway" $ nf (fromRationalTiesToAway :: Rational -> Double) x+ , bench "fromRationalTowardPositive" $ nf (fromRationalTowardPositive :: Rational -> Double) x+ , bench "fromRationalTowardNegative" $ nf (fromRationalTowardNegative :: Rational -> Double) x+ , bench "fromRationalTowardZero" $ nf (fromRationalTowardZero :: Rational -> Double) x+ ]+ , bgroup "Float"+ [ bench "stock" $ nf (fromRational :: Rational -> Float) x+ , bench "fromRationalTiesToEven" $ nf (fromRationalTiesToEven :: Rational -> Float) x+ , bench "fromRationalTiesToAway" $ nf (fromRationalTiesToAway :: Rational -> Float) x+ , bench "fromRationalTowardPositive" $ nf (fromRationalTowardPositive :: Rational -> Float) x+ , bench "fromRationalTowardNegative" $ nf (fromRationalTowardNegative :: Rational -> Float) x+ , bench "fromRationalTowardZero" $ nf (fromRationalTowardZero :: Rational -> Float) x+ ]+ ]+ , let x = 355 / 113 :: Rational+ in bgroup "small/small"+ [ bgroup "Double"+ [ bench "stock" $ nf (fromRational :: Rational -> Double) x+ , bench "fromRationalTiesToEven" $ nf (fromRationalTiesToEven :: Rational -> Double) x+ , bench "fromRationalTiesToAway" $ nf (fromRationalTiesToAway :: Rational -> Double) x+ , bench "fromRationalTowardPositive" $ nf (fromRationalTowardPositive :: Rational -> Double) x+ , bench "fromRationalTowardNegative" $ nf (fromRationalTowardNegative :: Rational -> Double) x+ , bench "fromRationalTowardZero" $ nf (fromRationalTowardZero :: Rational -> Double) x+ ]+ , bgroup "Float"+ [ bench "stock" $ nf (fromRational :: Rational -> Float) x+ , bench "fromRationalTiesToEven" $ nf (fromRationalTiesToEven :: Rational -> Float) x+ , bench "fromRationalTiesToAway" $ nf (fromRationalTiesToAway :: Rational -> Float) x+ , bench "fromRationalTowardPositive" $ nf (fromRationalTowardPositive :: Rational -> Float) x+ , bench "fromRationalTowardNegative" $ nf (fromRationalTowardNegative :: Rational -> Float) x+ , bench "fromRationalTowardZero" $ nf (fromRationalTowardZero :: Rational -> Float) x+ ]+ ]+ , let x = 0x1.deafbeefcafec0ffeep100 :: Rational+ in bgroup "binary"+ [ bgroup "Double"+ [ bench "stock" $ nf (fromRational :: Rational -> Double) x+ , bench "fromRationalTiesToEven" $ nf (fromRationalTiesToEven :: Rational -> Double) x+ , bench "fromRationalTiesToAway" $ nf (fromRationalTiesToAway :: Rational -> Double) x+ , bench "fromRationalTowardPositive" $ nf (fromRationalTowardPositive :: Rational -> Double) x+ , bench "fromRationalTowardNegative" $ nf (fromRationalTowardNegative :: Rational -> Double) x+ , bench "fromRationalTowardZero" $ nf (fromRationalTowardZero :: Rational -> Double) x+ ]+ , bgroup "Float"+ [ bench "stock" $ nf (fromRational :: Rational -> Float) x+ , bench "fromRationalTiesToEven" $ nf (fromRationalTiesToEven :: Rational -> Float) x+ , bench "fromRationalTiesToAway" $ nf (fromRationalTiesToAway :: Rational -> Float) x+ , bench "fromRationalTowardPositive" $ nf (fromRationalTowardPositive :: Rational -> Float) x+ , bench "fromRationalTowardNegative" $ nf (fromRationalTowardNegative :: Rational -> Float) x+ , bench "fromRationalTowardZero" $ nf (fromRationalTowardZero :: Rational -> Float) x+ ]+ ]+ ]+ , bgroup "encodeFloat"+ [ let arg = (0xcafe_0000_abcd_7777, -25) :: (Integer, Int)+ in bgroup "Double"+ [ bench "stock" $ nf (uncurry encodeFloat :: (Integer, Int) -> Double) arg+ , bench "encodeFloatTiesToEven" $ nf (uncurry encodeFloatTiesToEven :: (Integer, Int) -> Double) arg+ , bench "encodeFloatTiesToAway" $ nf (uncurry encodeFloatTiesToAway :: (Integer, Int) -> Double) arg+ , bench "encodeFloatTowardPositive" $ nf (uncurry encodeFloatTowardPositive :: (Integer, Int) -> Double) arg+ , bench "encodeFloatTowardNegative" $ nf (uncurry encodeFloatTowardNegative :: (Integer, Int) -> Double) arg+ , bench "encodeFloatTowardZero" $ nf (uncurry encodeFloatTowardZero :: (Integer, Int) -> Double) arg+ ]+ , let arg = (0xcafe_0000_abcd_7777, -25) :: (Integer, Int)+ in bgroup "Float"+ [ bench "stock" $ nf (uncurry encodeFloat :: (Integer, Int) -> Float) arg+ , bench "encodeFloatTiesToEven" $ nf (uncurry encodeFloatTiesToEven :: (Integer, Int) -> Float) arg+ , bench "encodeFloatTiesToAway" $ nf (uncurry encodeFloatTiesToAway :: (Integer, Int) -> Float) arg+ , bench "encodeFloatTowardPositive" $ nf (uncurry encodeFloatTowardPositive :: (Integer, Int) -> Float) arg+ , bench "encodeFloatTowardNegative" $ nf (uncurry encodeFloatTowardNegative :: (Integer, Int) -> Float) arg+ , bench "encodeFloatTowardZero" $ nf (uncurry encodeFloatTowardZero :: (Integer, Int) -> Float) arg+ ]+ ]+ , bgroup "minimum"+ [ bgroup "Double"+ [ let arg = (pi, -2.3) :: (Double, Double)+ in bgroup "(pi, -2.3)"+ [ bench "stock" $ whnf (uncurry min) arg+ , bench "minimum" $ whnf (uncurry minimum') arg+ , bench "minimumNumber" $ whnf (uncurry minimumNumber) arg+ , bench "minimumMagnitude" $ whnf (uncurry minimumMagnitude) arg+ , bench "minimumMagnitudeNumber" $ whnf (uncurry minimumMagnitudeNumber) arg+ , bench "minimum (specialized)" $ whnf (uncurry minimumDouble) arg+ , bench "minimumNumber (specialized)" $ whnf (uncurry minimumNumberDouble) arg+ ]+ , let arg = (0, -0) :: (Double, Double)+ in bgroup "(0, -0)"+ [ bench "stock" $ whnf (uncurry min) arg+ , bench "minimum" $ whnf (uncurry minimum') arg+ , bench "minimumNumber" $ whnf (uncurry minimumNumber) arg+ , bench "minimumMagnitude" $ whnf (uncurry minimumMagnitude) arg+ , bench "minimumMagnitudeNumber" $ whnf (uncurry minimumMagnitudeNumber) arg+ , bench "minimum (specialized)" $ whnf (uncurry minimumDouble) arg+ , bench "minimumNumber (specialized)" $ whnf (uncurry minimumNumberDouble) arg+ ]+ ]+ ]+ , bgroup "canonicalize"+ [ let x = 0 / 0 :: Float+ in bgroup "Float"+ [ bench "Haskell" $ whnf canonicalize x+ , bench "Haskell (generic)" $ whnf canonicalize (Identity x)+ , bench "C" $ whnf canonicalizeFloat x+ , bench "identity" $ whnf id x+ ]+ , let x = 0 / 0 :: Double+ in bgroup "Double"+ [ bench "Haskell" $ whnf canonicalize x+ , bench "Haskell (generic)" $ whnf canonicalize (Identity x)+ , bench "C" $ whnf canonicalizeDouble x+ , bench "identity" $ whnf id x+ ]+ ]+ , bgroup "nextUp"+ [ let cases = [0,1,0x1.ffff_ffff_ffff_fp200] :: [Double]+ in bgroup "Double"+ [ bgroup "C"+ [ bench (showHFloat x "") $ nf c_nextUp x | x <- cases ]+ , bgroup "Haskell"+ [ bench (showHFloat x "") $ nf nextUp x | x <- cases ]+ , bgroup "Haskell (generic)"+ [ bench (showHFloat x "") $ nf nextUp (Identity x) | x <- cases ]+ ]+ , let cases = [0,1,0x1.fffffep100] :: [Float]+ in bgroup "Float"+ [ bgroup "C"+ [ bench (showHFloat x "") $ nf c_nextUp x | x <- cases ]+ , bgroup "Haskell"+ [ bench (showHFloat x "") $ nf nextUp x | x <- cases ]+ , bgroup "Haskell (generic)"+ [ bench (showHFloat x "") $ nf nextUp (Identity x) | x <- cases ]+ ]+ ]+ , bgroup "nextDown"+ [ let cases = [0,1,0x1.ffff_ffff_ffff_fp200] :: [Double]+ in bgroup "Double"+ [ bgroup "C"+ [ bench (showHFloat x "") $ nf c_nextDown x | x <- cases ]+ , bgroup "Haskell"+ [ bench (showHFloat x "") $ nf nextDown x | x <- cases ]+ , bgroup "Haskell (generic)"+ [ bench (showHFloat x "") $ nf nextDown (Identity x) | x <- cases ]+ ]+ , let cases = [0,1,0x1.fffffep100] :: [Float]+ in bgroup "Float"+ [ bgroup "C"+ [ bench (showHFloat x "") $ nf c_nextDown x | x <- cases ]+ , bgroup "Haskell"+ [ bench (showHFloat x "") $ nf nextDown x | x <- cases ]+ , bgroup "Haskell (generic)"+ [ bench (showHFloat x "") $ nf nextDown (Identity x) | x <- cases ]+ ]+ ]+#if defined(USE_HALF)+ , bgroup "Half"+ [ bgroup "from Half"+ [ let x = 1.3 :: Half+ in bgroup "to Float"+ [ bench "half" $ nf fromHalf x+#if defined(HAS_FAST_HALF_CONVERSION)+ , bench "C impl" $ nf halfToFloat x+#endif+ , bench "realToFrac" $ nf (realToFrac :: Half -> Float) x+ , bench "realFloatToFrac" $ nf (realFloatToFrac :: Half -> Float) x+ ]+ , let x = 1.3 :: Half+ in bgroup "to Double"+ [+#if defined(HAS_FAST_HALF_CONVERSION)+ bench "C impl" $ nf halfToDouble x ,+#endif+ bench "realToFrac" $ nf (realToFrac :: Half -> Double) x+ , bench "realFloatToFrac" $ nf (realFloatToFrac :: Half -> Double) x+ ]+ ]+ , bgroup "to Half"+ [ let x = 1.3 :: Float+ in bgroup "from Float"+ [ bench "half" $ nf toHalf x+#if defined(HAS_FAST_HALF_CONVERSION)+ , bench "C impl" $ nf floatToHalf x+#endif+ , bench "realToFrac" $ nf (realToFrac :: Float -> Half) x+ , bench "realFloatToFrac" $ nf (realFloatToFrac :: Float -> Half) x+ ]+ , let x = 1.3 :: Double+ in bgroup "from Double"+ [+#if defined(HAS_FAST_HALF_CONVERSION)+ bench "C impl" $ nf doubleToHalf x ,+#endif+ bench "realToFrac" $ nf (realToFrac :: Double -> Half) x+ , bench "realFloatToFrac" $ nf (realFloatToFrac :: Double -> Half) x+ ]+ ]+ , let arg = pi :: Half+ in bgroup "isNormal"+ [ bench "default" $ nf isNormal arg+ , bench "generic" $ nf (isNormal . Identity) arg+ ]+ , let arg = pi :: Half+ in bgroup "isFinite"+ [ bench "default" $ nf isFinite arg+ , bench "generic" $ nf (isFinite . Identity) arg+ ]+ , let arg = -0 :: Half+ in bgroup "isZero"+ [ bench "default" $ nf isZero arg+ , bench "generic" $ nf (isZero . Identity) arg+ , bench "Numeric.Half.isZero" $ nf Numeric.Half.isZero arg+ ]+ ]+#endif+#if defined(USE_FLOAT128)+ , bgroup "Float128"+ [ bgroup "nextUp"+ [ bench "default" $ whnf nextUp (1.23 :: Float128)+ , bench "generic" $ whnf (nextUp . Identity) (1.23 :: Float128)+ ]+ , bgroup "nextDown"+ [ bench "default" $ whnf nextDown (1.23 :: Float128)+ , bench "generic" $ whnf (nextDown . Identity) (1.23 :: Float128)+ ]+ , bgroup "nextTowardZero"+ [ bench "default" $ whnf nextTowardZero (1.23 :: Float128)+ , bench "generic" $ whnf (nextTowardZero . Identity) (1.23 :: Float128)+ ]+ , bgroup "isNormal"+ [ bench "default" $ whnf isNormal (1.23 :: Float128)+ , bench "generic" $ whnf (isNormal . Identity) (1.23 :: Float128)+ ]+ , bgroup "isFinite"+ [ bench "default" $ whnf isFinite (1.23 :: Float128)+ , bench "generic" $ whnf (isFinite . Identity) (1.23 :: Float128)+ ]+ , bgroup "classify"+ [ bench "default" $ whnf classify (1.23 :: Float128)+ , bench "generic" $ whnf (classify . Identity) (1.23 :: Float128)+ ]+ , bgroup "isMantissaEven"+ [ bench "default" $ whnf isMantissaEven (1.23 :: Float128)+ , bench "generic" $ whnf (isMantissaEven . Identity) (1.23 :: Float128)+ ]+ , bgroup "roundAway"+ [ bench "default" $ whnf roundAway' (1.23 :: Float128)+ , bench "generic" $ whnf (roundAway' . Identity) (1.23 :: Float128)+ , bench "default (as Integer)" $ whnf (roundAway :: Float128 -> Integer) (1.23 :: Float128)+ , bench "generic (as Integer)" $ whnf ((roundAway :: Identity Float128 -> Integer) . Identity) (1.23 :: Float128)+ ]+ , bgroup "floor"+ [ bench "default" $ whnf floor' (1.23 :: Float128)+ , bench "generic" $ whnf (floor' . Identity) (1.23 :: Float128)+ , bench "default (as Integer)" $ whnf (floor :: Float128 -> Integer) (1.23 :: Float128)+ , bench "generic (as Integer)" $ whnf ((floor :: Identity Float128 -> Integer) . Identity) (1.23 :: Float128)+ ]+ ]+#endif+ ]
+ cbits/canonicalize.c view
@@ -0,0 +1,64 @@+#include <math.h>++#pragma STDC FENV_ACCESS ON++#if defined(__SSE2__)++#include <x86intrin.h>++float hs_canonicalizeFloat(float x)+{+ asm volatile("mulss %1, %0" : "+x"(x) : "x"(1.0f));+ return x;+ /*+ Clang optimizes away this:+ __m128 xv = _mm_set_ss(x);+ __m128 onev = _mm_set_ss(1.0f);+ __m128 resultv = _mm_mul_ss(xv, onev);+ float result;+ _mm_store_ss(&result, resultv);+ return result;+ */+}+double hs_canonicalizeDouble(double x)+{+ asm volatile("mulsd %1, %0" : "+x"(x) : "x"(1.0));+ return x;+ /*+ Clang optimizes away this:+ __m128d xv = _mm_set_sd(x);+ __m128d onev = _mm_set_sd(1.0);+ __m128d resultv = _mm_mul_sd(xv, onev);+ double result;+ _mm_store_sd(&result, resultv);+ return result;+ */+}++#elif defined(__aarch64__)++float hs_canonicalizeFloat(float x)+{+ 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));+ return x;+}++#else++float hs_canonicalizeFloat(float x)+{+ volatile float one = 1.0f;+ return x * one;+}+double hs_canonicalizeDouble(double x)+{+ volatile double one = 1.0;+ return x * one;+}++#endif
+ cbits/fma.c view
@@ -0,0 +1,17 @@+#include <math.h>++#if !defined(FP_FAST_FMA)+#error "The compiler should define FP_FAST_FMA"+#endif+#if !defined(FP_FAST_FMAF)+#error "The compiler should define FP_FAST_FMAF"+#endif++double hs_fusedMultiplyAddDouble(double a, double b, double c)+{+ return fma(a, b, c);+}+float hs_fusedMultiplyAddFloat(float a, float b, float c)+{+ return fmaf(a, b, c);+}
+ cbits/half.c view
@@ -0,0 +1,120 @@+#include <stdint.h> // uint16_t+#include <math.h>++#if defined(__F16C__) // x86 F16C++#include <x86intrin.h>++uint16_t hs_fastFloatToHalf(float f)+{+ __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;+}++float hs_fastHalfToFloat(uint16_t c)+{+ union {+ __m128i v;+ uint16_t c;+ } u;+ u.c = c;+ __m128 w = _mm_cvtph_ps(u.v); // VCVTPH2PS+ float d;+ _mm_store_ss(&d, w);+ return d;+}++// Is this really faster than bit manipulation?+uint16_t hs_fastDoubleToHalf(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;+ }+ __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)+{+ union {+ __m128i v;+ uint16_t c;+ } u;+ u.c = c;+ __m128 w = _mm_cvtph_ps(u.v); // VCVTPH2PS+ float d;+ _mm_store_ss(&d, w);+ return (double)d;+}++#else++// Let's hope _Float16 is available++uint16_t hs_fastFloatToHalf(float x)+{+ union {+ _Float16 f;+ uint16_t u;+ } u;+ u.f = (_Float16)x;+ return u.u;+}++float hs_fastHalfToFloat(uint16_t x)+{+ union {+ _Float16 f;+ uint16_t u;+ } u;+ u.u = x;+ return (float)u.f;+}++uint16_t hs_fastDoubleToHalf(double x)+{+ union {+ _Float16 f;+ uint16_t u;+ } u;+ u.f = (_Float16)x;+ return u.u;+}++double hs_fastHalfToDouble(uint16_t x)+{+ union {+ _Float16 f;+ uint16_t u;+ } u;+ u.u = x;+ return (double)u.f;+}++#endif
+ cbits/minmax.c view
@@ -0,0 +1,102 @@++// In case of GCC, -fsignaling-nans must be set to use '*= 1.0' as canonicalization+// #if defined(__GNUC__) && !defined(__SUPPORT_SNAN__)+// #error "-fsignaling-nans must be set"+// #endif++#if defined(__aarch64__)++// Properties of minimum and maximum:+// * -0 < +0+// * If either of inputs is NaN, returns a quiet NaN.++float hs_minimumFloat(float x, float y)+{+ float result;+ 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));+ return result;+}++double hs_minimumDouble(double x, double y)+{+ double result;+ 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));+ return result;+}++// Properties of minimumNumber and maximumNumber:+// * -0 < +0+// * Treat a NaN as "lack of input".+// If both of inputs are NaNs, returns a quiet NaN.++float hs_minimumNumberFloat(float x, float y)+{+ float result;+ // FMINNM always returns a NaN if either of inputs is signaling NaN.+ // 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));+ return result;+}++float hs_maximumNumberFloat(float x, float y)+{+ float result;+ // FMAXNM always returns a NaN if either of inputs is signaling NaN.+ // 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));+ return result;+}++double hs_minimumNumberDouble(double x, double y)+{+ double result;+ // FMINNM always returns a NaN if either of inputs is signaling NaN.+ // 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));+ return result;+}++double hs_maximumNumberDouble(double x, double y)+{+ double result;+ // FMAXNM always returns a NaN if either of inputs is signaling NaN.+ // 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));+ return result;+}++#else++#error "Unsupported platform"++#endif
+ cbits/roundeven.c view
@@ -0,0 +1,48 @@+#include <math.h>+#include <fenv.h>++#if defined(__SSE4_1__) // SSE 4.1++#include <x86intrin.h>++float hs_roundevenFloat(float x)+{+ __m128 xv = _mm_set_ss(x);+ xv = _mm_round_ss(xv, xv, _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC);+ float result;+ _mm_store_ss(&result, xv);+ return result;+}++double hs_roundevenDouble(double x)+{+ __m128d xv = _mm_set_sd(x);+ xv = _mm_round_sd(xv, xv, _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC);+ double result;+ _mm_store_sd(&result, xv);+ return result;+}++#elif defined(__aarch64__) // ARMv8-A++float hs_roundevenFloat(float x)+{+ float result;+ // a floating-exception can be generated+ asm("frintn %s0, %s1" : "=w"(result) : "w"(x));+ return result;+}++double hs_roundevenDouble(double x)+{+ double result;+ // a floating-exception can be generated+ asm("frintn %d0, %d1" : "=w"(result) : "w"(x));+ return result;+}++#else++#error "Unsupported architecture"++#endif
+ decimal-test/NextFloatSpec.hs view
@@ -0,0 +1,42 @@+module NextFloatSpec where+import Data.Proxy+import Numeric.Decimal+import Numeric.Floating.IEEE+import Test.Hspec+import Test.Hspec.QuickCheck (prop)+import Test.QuickCheck+import Util (forAllFloats, sameFloatP)++isPositiveZero :: RealFloat a => a -> Bool+isPositiveZero x = x == 0 && not (isNegativeZero x)++prop_nextUp_nextDown :: (RealFloat a, Show a) => Proxy a -> a -> Property+prop_nextUp_nextDown _ x = x /= (-1/0) ==>+ let x' = nextUp (nextDown x)+ in x' `sameFloatP` x .||. (isPositiveZero x .&&. isNegativeZero x')++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')++{-# NOINLINE spec #-}+spec :: Spec+spec = do+ describe "Decimal32" $ do+ let proxy :: Proxy Decimal32+ proxy = Proxy+ prop "nextUp . nextDown == id (unless -inf)" $ forAllFloats $ prop_nextUp_nextDown proxy+ prop "nextDown . nextUp == id (unless inf)" $ forAllFloats $ prop_nextDown_nextUp proxy++ describe "Decimal64" $ do+ let proxy :: Proxy Decimal64+ proxy = Proxy+ prop "nextUp . nextDown == id (unless -inf)" $ forAllFloats $ prop_nextUp_nextDown proxy+ prop "nextDown . nextUp == id (unless inf)" $ forAllFloats $ prop_nextDown_nextUp proxy++ describe "Decimal128" $ do+ let proxy :: Proxy Decimal128+ proxy = Proxy+ prop "nextUp . nextDown == id (unless -inf)" $ forAllFloats $ prop_nextUp_nextDown proxy+ prop "nextDown . nextUp == id (unless inf)" $ forAllFloats $ prop_nextDown_nextUp proxy
+ decimal-test/Spec.hs view
@@ -0,0 +1,27 @@+{-# LANGUAGE NumericUnderscores #-}+import qualified NextFloatSpec+import Numeric.Decimal+import Numeric.Floating.IEEE+import Test.Hspec+import Test.Hspec.Core.Spec++allowFailure :: String -> Item a -> Item a+allowFailure message item@(Item { itemExample = origExample }) = item { itemExample = newExample }+ where+ newExample params around callback = do+ result <- origExample params around callback+ case result of+ Result { resultStatus = Failure loc reason } -> do+ let message' = case reason of+ NoReason -> message+ _ -> message ++ ": " ++ show reason+ return result { resultStatus = Pending loc (Just message') }+ _ -> return result++main :: IO ()+main = hspec $ do+ mapSpecItem_ (allowFailure "decimal-arithmetic's floatRange may be incorrect") $ do+ it "maxFinite :: Decimal32" $ (maxFinite :: Decimal32) == 9.999_999e96 -- 7 digits+ it "maxFinite :: Decimal64" $ (maxFinite :: Decimal64) == 9.999_999_999_999_999e384 -- 16 digits+ it "maxFinite :: Decimal128" $ (maxFinite :: Decimal128) == 9.999_999_999_999_999_999_999_999_999_999_999e6144 -- 34 digits+ describe "NextFloat" NextFloatSpec.spec
+ doctests.hs view
@@ -0,0 +1,13 @@+import Test.DocTest++main :: IO ()+main = doctest [ "-isrc"+ , "src/Numeric/Floating/IEEE/Internal/Base.hs"+ , "src/Numeric/Floating/IEEE/Internal/Classify.hs"+ , "src/Numeric/Floating/IEEE/Internal/FMA.hs"+ , "src/Numeric/Floating/IEEE/Internal/GenericArith.hs"+ , "src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs"+ , "src/Numeric/Floating/IEEE/Internal/MinMax.hs"+ , "src/Numeric/Floating/IEEE/Internal/NextFloat.hs"+ , "src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs"+ ]
+ fp-ieee.cabal view
@@ -0,0 +1,435 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.33.0.+--+-- see: https://github.com/sol/hpack+--+-- hash: f0bf89e9df957b5023d91e55d09165cecb06208750fddaf58af69fd7c2b7e35d++name: fp-ieee+version: 0.1.0+description: Please see the README on GitHub at <https://github.com/minoki/haskell-floating-point/tree/master/fp-ieee#readme>+category: Numeric, Math+homepage: https://github.com/minoki/haskell-floating-point#readme+bug-reports: https://github.com/minoki/haskell-floating-point/issues+author: ARATA Mizuki+maintainer: minorinoki@gmail.com+copyright: 2020 ARATA Mizuki+license: BSD3+license-file: LICENSE+build-type: Simple+extra-source-files:+ README.md+ ChangeLog.md++source-repository head+ type: git+ location: https://github.com/minoki/haskell-floating-point++flag f16c+ description: Use F16C instructions on x86+ manual: True+ default: False++flag float128+ description: Support Float128 via float128 package+ manual: True+ default: False++flag fma3+ description: Use FMA3 instructions on x86+ manual: True+ default: False++flag ghc-bignum+ description: Use ghc-bignum package+ manual: False+ default: True++flag half+ description: Support Half (float16) via half package+ manual: True+ default: False++flag integer-gmp+ description: Use integer-gmp package+ manual: False+ default: True++flag pure-hs+ description: Disable FFI+ manual: True+ default: False++flag sse4_1+ description: Use SSE4.1 instructions on x86+ manual: True+ default: False++library+ exposed-modules:+ Numeric.Floating.IEEE+ Numeric.Floating.IEEE.Internal+ Numeric.Floating.IEEE.NaN+ other-modules:+ GHC.Float.Compat+ MyPrelude+ Numeric.Floating.IEEE.Internal.Augmented+ Numeric.Floating.IEEE.Internal.Base+ Numeric.Floating.IEEE.Internal.Classify+ Numeric.Floating.IEEE.Internal.Conversion+ Numeric.Floating.IEEE.Internal.FMA+ Numeric.Floating.IEEE.Internal.GenericArith+ Numeric.Floating.IEEE.Internal.IntegerInternals+ Numeric.Floating.IEEE.Internal.MinMax+ Numeric.Floating.IEEE.Internal.NaN+ Numeric.Floating.IEEE.Internal.NextFloat+ Numeric.Floating.IEEE.Internal.Remainder+ Numeric.Floating.IEEE.Internal.RoundToIntegral+ Numeric.Floating.IEEE.Internal.Rounding+ Numeric.Floating.IEEE.Internal.Rounding.Common+ Numeric.Floating.IEEE.Internal.Rounding.Encode+ Numeric.Floating.IEEE.Internal.Rounding.Integral+ Numeric.Floating.IEEE.Internal.Rounding.Rational+ hs-source-dirs:+ src+ ghc-options: -Wall+ build-depends:+ base >=4.12 && <5+ , integer-logarithms >=1 && <1.1+ if arch(i386)+ ghc-options: -msse2+ cc-options: -msse2 -mfpmath=sse+ if !flag(pure-hs)+ cpp-options: -DUSE_FFI+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ cc-options: -msse4.1+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3)+ cpp-options: -DHAS_FAST_FMA+ cc-options: -mfma+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_FMA+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && !os(windows)+ cpp-options: -DUSE_C99_FMA+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_MINMAX+ c-sources:+ cbits/minmax.c+ if flag(float128)+ cpp-options: -DUSE_FLOAT128+ build-depends:+ float128 >=0.1 && <0.2+ if flag(half)+ cpp-options: -DUSE_HALF+ build-depends:+ half >=0.3 && <0.4+ if !flag(pure-hs) && flag(half) && arch(x86_64) && flag(f16c)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ cc-options: -mf16c+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && flag(half) && arch(aarch64)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && (arch(aarch64) || arch(x86_64))+ cpp-options: -DHAS_FAST_CANONICALIZE+ c-sources:+ cbits/canonicalize.c+ if flag(half)+ other-modules:+ Numeric.Floating.IEEE.Internal.Half+ if flag(float128)+ other-modules:+ Numeric.Floating.IEEE.Internal.Float128+ if flag(integer-gmp) && impl(ghc < 9.0.0)+ build-depends:+ integer-gmp >=1.0 && <1.1+ if flag(ghc-bignum) && impl(ghc >= 9.0.0)+ build-depends:+ ghc-bignum >=1.0 && <1.1+ default-language: Haskell2010++test-suite fp-ieee-decimal-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ NextFloatSpec+ hs-source-dirs:+ decimal-test+ ghc-options: -threaded -rtsopts -with-rtsopts=-N -fno-ignore-asserts+ build-depends:+ QuickCheck+ , base >=4.12 && <5+ , decimal-arithmetic+ , fp-ieee+ , hspec+ , hspec-core+ , random+ if arch(i386)+ ghc-options: -msse2+ cc-options: -msse2 -mfpmath=sse+ if !flag(pure-hs)+ cpp-options: -DUSE_FFI+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ cc-options: -msse4.1+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3)+ cpp-options: -DHAS_FAST_FMA+ cc-options: -mfma+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_FMA+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && !os(windows)+ cpp-options: -DUSE_C99_FMA+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_MINMAX+ c-sources:+ cbits/minmax.c+ if flag(float128)+ cpp-options: -DUSE_FLOAT128+ build-depends:+ float128 >=0.1 && <0.2+ if flag(half)+ cpp-options: -DUSE_HALF+ build-depends:+ half >=0.3 && <0.4+ if !flag(pure-hs) && flag(half) && arch(x86_64) && flag(f16c)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ cc-options: -mf16c+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && flag(half) && arch(aarch64)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && (arch(aarch64) || arch(x86_64))+ cpp-options: -DHAS_FAST_CANONICALIZE+ c-sources:+ cbits/canonicalize.c+ default-language: Haskell2010++test-suite fp-ieee-doctests+ type: exitcode-stdio-1.0+ main-is: doctests.hs+ other-modules:+ Paths_fp_ieee+ build-depends:+ base >=4.12 && <5+ , doctest >=0.8+ if arch(i386)+ ghc-options: -msse2+ cc-options: -msse2 -mfpmath=sse+ if !flag(pure-hs)+ cpp-options: -DUSE_FFI+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ cc-options: -msse4.1+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3)+ cpp-options: -DHAS_FAST_FMA+ cc-options: -mfma+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_FMA+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && !os(windows)+ cpp-options: -DUSE_C99_FMA+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_MINMAX+ c-sources:+ cbits/minmax.c+ if flag(float128)+ cpp-options: -DUSE_FLOAT128+ build-depends:+ float128 >=0.1 && <0.2+ if flag(half)+ cpp-options: -DUSE_HALF+ build-depends:+ half >=0.3 && <0.4+ if !flag(pure-hs) && flag(half) && arch(x86_64) && flag(f16c)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ cc-options: -mf16c+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && flag(half) && arch(aarch64)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && (arch(aarch64) || arch(x86_64))+ cpp-options: -DHAS_FAST_CANONICALIZE+ c-sources:+ cbits/canonicalize.c+ default-language: Haskell2010++test-suite fp-ieee-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ AugmentedArithSpec+ ClassificationSpec+ FMASpec+ IntegerInternalsSpec+ MinMaxSpec+ NaNSpec+ RoundingSpec+ RoundToIntegralSpec+ TwoSumSpec+ hs-source-dirs:+ test+ ghc-options: -threaded -rtsopts -with-rtsopts=-N -fno-ignore-asserts+ build-depends:+ QuickCheck+ , base >=4.12 && <5+ , fp-ieee+ , hspec+ , hspec-core+ , integer-logarithms+ , random+ if arch(i386)+ ghc-options: -msse2+ cc-options: -msse2 -mfpmath=sse+ if !flag(pure-hs)+ cpp-options: -DUSE_FFI+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ cc-options: -msse4.1+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3)+ cpp-options: -DHAS_FAST_FMA+ cc-options: -mfma+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_FMA+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && !os(windows)+ cpp-options: -DUSE_C99_FMA+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_MINMAX+ c-sources:+ cbits/minmax.c+ if flag(float128)+ cpp-options: -DUSE_FLOAT128+ build-depends:+ float128 >=0.1 && <0.2+ if flag(half)+ cpp-options: -DUSE_HALF+ build-depends:+ half >=0.3 && <0.4+ if !flag(pure-hs) && flag(half) && arch(x86_64) && flag(f16c)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ cc-options: -mf16c+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && flag(half) && arch(aarch64)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && (arch(aarch64) || arch(x86_64))+ cpp-options: -DHAS_FAST_CANONICALIZE+ c-sources:+ cbits/canonicalize.c+ if flag(half)+ other-modules:+ HalfSpec+ if flag(float128)+ other-modules:+ Float128Spec+ default-language: Haskell2010++benchmark fp-ieee-benchmark+ type: exitcode-stdio-1.0+ main-is: Benchmark.hs+ other-modules:+ Paths_fp_ieee+ hs-source-dirs:+ benchmark+ build-depends:+ base >=4.12 && <5+ , fp-ieee+ , gauge+ if arch(i386)+ ghc-options: -msse2+ cc-options: -msse2 -mfpmath=sse+ if !flag(pure-hs)+ cpp-options: -DUSE_FFI+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(sse4_1)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ cc-options: -msse4.1+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_ROUNDEVEN+ c-sources:+ cbits/roundeven.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && flag(fma3)+ cpp-options: -DHAS_FAST_FMA+ cc-options: -mfma+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_FMA+ c-sources:+ cbits/fma.c+ if !flag(pure-hs) && (arch(i386) || arch(x86_64)) && !os(windows)+ cpp-options: -DUSE_C99_FMA+ if !flag(pure-hs) && arch(aarch64)+ cpp-options: -DHAS_FAST_MINMAX+ c-sources:+ cbits/minmax.c+ if flag(float128)+ cpp-options: -DUSE_FLOAT128+ build-depends:+ float128 >=0.1 && <0.2+ if flag(half)+ cpp-options: -DUSE_HALF+ build-depends:+ half >=0.3 && <0.4+ if !flag(pure-hs) && flag(half) && arch(x86_64) && flag(f16c)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ cc-options: -mf16c+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && flag(half) && arch(aarch64)+ cpp-options: -DHAS_FAST_HALF_CONVERSION+ c-sources:+ cbits/half.c+ if !flag(pure-hs) && (arch(aarch64) || arch(x86_64))+ cpp-options: -DHAS_FAST_CANONICALIZE+ c-sources:+ cbits/canonicalize.c+ default-language: Haskell2010
+ src/GHC/Float/Compat.hs view
@@ -0,0 +1,26 @@+{-# LANGUAGE CPP #-}++-- castFloatToWord32 is buggy on GHC <= 8.8 && 64-bit systems.+-- See https://gitlab.haskell.org/ghc/ghc/issues/16617++#include "MachDeps.h"++#if MIN_VERSION_base(4,14,0) || WORD_SIZE_IN_BITS == 32++module GHC.Float.Compat (module GHC.Float) where+import GHC.Float++#else++module GHC.Float.Compat (module GHC.Float, castFloatToWord32) where+import GHC.Float hiding (castFloatToWord32)+import qualified GHC.Float as F+import Data.Bits ((.&.))+import Data.Word (Word32)++-- Let's hope the compiler is not smart enough to eliminate the bit-and...+-- Or @fromIntegral (fromIntegral x :: Int) :: Word32@ might be better?+castFloatToWord32 :: Float -> Word32+castFloatToWord32 x = F.castFloatToWord32 x .&. 0xFFFFFFFF++#endif
+ src/MyPrelude.hs view
@@ -0,0 +1,14 @@+{-+This module is the custom Prelude for this project.+You can replace Prelude's definition by a debugging-friendly one.+Examples are:++type RealFloat a = (Prelude.RealFloat a, Show a)++(^) :: (HasCallStack, Num a, Integral b) => a -> b -> a+x ^ y | y < 0 = error "Negative exponent" -- with stack trace+ | otherwise = x Prelude.^ y+-}++module MyPrelude (module Prelude) where+import Prelude
+ src/Numeric/Floating/IEEE.hs view
@@ -0,0 +1,250 @@+{-|+Module : Numeric.Floating.IEEE+Description : IEEE 754-compliant operations for floating-point numbers++This module provides IEEE 754-compliant operations for floating-point numbers.++The functions in this module assume that the given floating-point type conform to IEEE 754 format.++Since 'RealFloat' constraint is insufficient to query properties of a NaN, the functions here assumes all NaN as positive, quiet.+If you want better treatment for NaNs, use the module "Numeric.Floating.IEEE.NaN".++Since floating-point exceptions cannot be accessed from Haskell, the operations provided by this module ignore exceptional behavior.+This library assumes the default exception handling is in use.++If you are using GHC <= 8.8 on i386 target, you may need to set @-msse2@ option to get correct floating-point behavior.+-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE+ (+ -- * Standard Haskell classes+ --+ -- $stdclasses++ -- * 5.3 Homogeneous general-computational operations+ --+ -- ** 5.3.1 General operations+ round'+ , roundAway'+ , truncate'+ , ceiling'+ , floor'+ , nextUp+ , nextDown+ , nextTowardZero -- not in IEEE+ , remainder++ -- ** 5.3.2 Decimal operations (not supported)+ --+ -- | Not supported.++ -- ** 5.3.3 logBFormat operations+ --+ , scaleFloatTiesToEven+ , scaleFloatTiesToAway+ , scaleFloatTowardPositive+ , scaleFloatTowardNegative+ , scaleFloatTowardZero+ -- |+ -- The Haskell counterpart for IEEE 754 @logB@ operation is 'exponent'.+ -- Note that @logB@ and 'exponent' are different by one:+ -- @logB x = 'exponent' x - 1@+ , exponent++ -- * 5.4 formatOf general-computational operations+ --+ -- ** 5.4.1 Arithmetic operations+ --+ -- |+ -- For IEEE-compliant floating-point types, '(+)', '(-)', '(*)', '(/)', and 'sqrt' from "Prelude" should be correctly-rounding.+ -- 'fusedMultiplyAdd' is provided by this library.+ -- This library also provides \"generic\" version of the arithmetic operations, which can be useful if the target type is narrower than source.+ , (+) -- addition+ , (-) -- subtraction+ , (*) -- multiplication+ , (/) -- division+ , sqrt -- squareRoot+ , fusedMultiplyAdd+ , genericAdd+ , genericSub+ , genericMul+ , genericDiv+ -- | @genericSqrt@ is not implemented yet.+ , genericFusedMultiplyAdd+ , fromIntegerTiesToEven+ , fromIntegerTiesToAway+ , fromIntegerTowardPositive+ , fromIntegerTowardNegative+ , fromIntegerTowardZero+ , fromIntegralTiesToEven+ , fromIntegralTiesToAway+ , fromIntegralTowardPositive+ , fromIntegralTowardNegative+ , fromIntegralTowardZero+ , fromRationalTiesToEven+ , fromRationalTiesToAway+ , fromRationalTowardPositive+ , fromRationalTowardNegative+ , fromRationalTowardZero+ , round -- convertToIntegerTiesToEven+ , roundAway -- convertToIntegerTiesToAway+ , truncate -- convertToIntegerTowardZero+ , ceiling -- convertToIntegerTowardPositive+ , floor -- convertToIntegerTowardNegative++ -- ** 5.4.2 Conversion operations for floating-point formats and decimal character sequences+ --+ -- |+ -- Unfortunately, 'realToFrac' does not have a good semantics, and behaves differently with rewrite rules (consider @realToFrac (0/0 :: Float) :: Double@).+ -- As an alternative, this library provides 'realFloatToFrac', with well-defined semantics on signed zeroes, infinities and NaNs.+ -- Like 'realToFrac', 'realFloatToFrac' comes with some rewrite rules for particular types, but they should not change behavior.+ , realFloatToFrac -- convertFormat+ , canonicalize+ -- |+ -- @convertFromDecimalCharacter@: not implemented.+ --+ -- @convertToDecimalCharacter@: not implemented.++ -- * 5.4.3 Conversion operations for binary formats+ --+ -- |+ -- @convertFromHexCharacter@: not implemented.+ --+ -- @convertToHexCharacter@: 'Numeric.showHFloat' from "Numeric" can be used.++ -- * 5.5 Quiet-computational operations+ --+ -- ** 5.5.1 Sign bit operations+ --+ -- |+ -- For IEEE-compliant floating-point types, 'negate' and 'abs' from "Prelude" should comply with IEEE semantics.+ , negate+ , abs+ -- |+ -- See "Numeric.Floating.IEEE.NaN" for @copySign@.++ -- ** 5.5.2 Decimal re-encoding operations (not supported)+ --+ -- |+ -- Not supported.++ -- * 5.6 Signaling-computational operations+ --+ -- ** 5.6.1 Comparisons (not supported)+ --+ -- |+ -- This library does not support floating-point exceptions.++ -- * 5.7 Non-computational operations+ --+ -- ** 5.7.1 Conformance predicates (not supported)+ --+ -- |+ -- Not supported.++ -- ** 5.7.2 General operations+ --+ -- |+ -- Functions in this module disregards the content of NaNs: sign bit, signaling-or-quiet, and payload.+ -- All NaNs are treated as quiet, positive.+ -- To properly handle NaNs, use the typeclass and functions from "Numeric.Floating.IEEE.NaN".+ , Class(..)+ , classify -- class+ , isSignMinus+ , isNormal+ , isFinite+ , isZero+ , isDenormalized -- isSubnormal+ , isInfinite -- re-export+ , isNaN -- re-export+ -- |+ -- See "Numeric.Floating.IEEE.NaN" for @isSignaling@.+ --+ -- @isCanonical@: not supported.+ , floatRadix -- radix+ , compareByTotalOrder -- totalOrder+ , compareByTotalOrderMag -- totalOrderMag++ -- ** 5.7.3 Decimal operation (not supported)+ --+ -- |+ -- Not supported.++ -- ** 5.7.4 Operations on subsets of flags (not supported)+ --+ -- |+ -- Not supported.++ -- * 9. Recommended operations++ -- * 9.5 Augmented arithmetic operations+ , augmentedAddition+ , augmentedSubtraction+ , augmentedMultiplication++ -- * 9.6 Minimum and maximum operations+ , minimum'+ , minimumNumber+ , maximum'+ , maximumNumber+ , minimumMagnitude+ , minimumMagnitudeNumber+ , maximumMagnitude+ , maximumMagnitudeNumber++ -- * Floating-point constants+ , minPositive+ , minPositiveNormal+ , maxFinite+ ) where+import MyPrelude+import Numeric.Floating.IEEE.Internal++-- $stdclasses+--+-- This library assumes that some of the standard numeric functions correspond to the operations specified by IEEE.+-- The rounding attribute should be roundTiesToEven and the exceptional behavior should be the default one.+--+-- == 'Num'+--+-- * '(+)', '(-)', and '(*)' should be correctly-rounding.+-- * 'negate', 'abs' should comply with IEEE semantics.+-- * 'fromInteger' should be correctly-rounding, but unfortunately not for 'Float' and 'Double' (see GHC's [#17231](https://gitlab.haskell.org/ghc/ghc/-/issues/17231)).+-- This module provides a correctly-rounding alternative: 'fromIntegerTiesToEven'.+--+-- == 'Fractional'+--+-- * '(/)' should be correctly-rounding.+-- * 'fromRational' should be correctly-rounding, but some third-partiy floating-point types fail to do so.+--+-- == 'Floating'+--+-- * 'sqrt' should be correctly-rounding.+--+-- == 'RealFrac'+--+-- * 'truncate': IEEE 754 @convertToIntegerTowardZero@ operation.+-- * 'round': IEEE 754 @convertToIntegerTiesToEven@ operation; the Language Report says that this should choose the even integer if the argument is the midpoint of two successive integers.+-- * 'ceiling': IEEE 754 @convertToIntegerTowardPositive@ operation.+-- * 'floor': IEEE 754 @convertToIntegerTowardNegative@ operation.+--+-- To complete these, 'roundAway' is provided by this library.+-- Note that Haskell's 'round' is specified to be ties-to-even, whereas C's @round@ is ties-to-away.+--+-- == 'RealFloat'+--+-- This class provides information on the IEEE-compliant format.+--+-- * 'floatRadix': The base \(b\). IEEE 754 @radix@ operation.+-- * 'floatDigits': The precision \(p\).+-- * 'floatRange': The exponent range offset by 1: \((\mathit{emin}+1,\mathit{emax}+1)\)+-- * @'decodeFloat' x@: The exponent part returned is in the range \([\mathit{emin}+1-p,\mathit{emax}+1-p]\) if @x@ is normal, or in \([\mathit{emin}-2p+2,\mathit{emin}-p]\) if @x@ is subnormal.+-- * 'encodeFloat' should accept the significand in the range @[0, floatRadix x ^ floatDigits x]@. This library does not assume a particular rounding behavior when the result cannot be expressed in the target type.+-- * @'exponent' x@: The exponent offset by 1: \(\mathrm{logB}(x)+1\). Returns an integer in \([\mathit{emin}+1,\mathit{emax}+1]\) if @x@ is normal, or in \([\mathit{emin}-p+2,\mathit{emin}]\) if @x@ is subnormal.+-- * @'significand' x@: Returns the significand of @x@ as a value between \([1/b,1)\).+-- * 'scaleFloat': This library does not assume a particular rounding behavior when the result is subnormal.+-- * 'isNaN'+-- * 'isInfinite'+-- * 'isDenormalized'+-- * 'isNegativeZero'+-- * 'isIEEE' should return @True@ if you are using the type with this library.
+ src/Numeric/Floating/IEEE/Internal.hs view
@@ -0,0 +1,23 @@+{-# LANGUAGE CPP #-}+{-# OPTIONS_HADDOCK hide #-}+module Numeric.Floating.IEEE.Internal+ ( module Internal+ ) where+import Numeric.Floating.IEEE.Internal.Augmented as Internal+import Numeric.Floating.IEEE.Internal.Base as Internal hiding ((^!))+import Numeric.Floating.IEEE.Internal.Classify as Internal+import Numeric.Floating.IEEE.Internal.Conversion as Internal+import Numeric.Floating.IEEE.Internal.FMA as Internal+import Numeric.Floating.IEEE.Internal.GenericArith as Internal+import Numeric.Floating.IEEE.Internal.IntegerInternals as Internal+import Numeric.Floating.IEEE.Internal.MinMax as Internal+import Numeric.Floating.IEEE.Internal.NextFloat as Internal+import Numeric.Floating.IEEE.Internal.Remainder as Internal+import Numeric.Floating.IEEE.Internal.Rounding as Internal+import Numeric.Floating.IEEE.Internal.RoundToIntegral as Internal+#if defined(USE_HALF)+import Numeric.Floating.IEEE.Internal.Half as Internal+#endif+#if defined(USE_FLOAT128)+import Numeric.Floating.IEEE.Internal.Float128 as Internal+#endif
+ src/Numeric/Floating/IEEE/Internal/Augmented.hs view
@@ -0,0 +1,127 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.Augmented where+import Control.Exception (assert)+import MyPrelude+import Numeric.Floating.IEEE.Internal.FMA (isMantissaEven,+ twoProduct_nonscaling,+ twoSum)+import Numeric.Floating.IEEE.Internal.NextFloat (nextDown,+ nextTowardZero,+ nextUp)++default ()++-- |+-- IEEE 754 @augmentedAddition@ operation.+augmentedAddition :: RealFloat a => a -> a -> (a, a)+augmentedAddition !x !y+ | isNaN x || isInfinite x || isNaN y || isInfinite y = let !result = x + y in (result, result)+ | otherwise = let (u1, u2) = twoSum x y+ ulpTowardZero = u1 - nextTowardZero u1+ in if isNaN u2 then+ -- Handle undue overflow: e.g. 0x1.ffff_ffff_ffff_f8p1023+ handleUndueOverflow+ else+ if u2 == 0 then+ (u1, 0 * u1) -- signed zero+ else+ if (-2) * u2 == ulpTowardZero then+ (u1 - ulpTowardZero, ulpTowardZero + u2)+ else+ (u1, u2)+ where+ handleUndueOverflow =+ -- The exponents of inputs should be close enough so that neither x' nor y' underflow.+ let e = max (exponent x) (exponent y)+ x' = scaleFloat (- e) x+ y' = scaleFloat (- e) y+ (u1, u2) = twoSum x' y'+ ulpTowardZero = u1 - nextTowardZero u1+ (v1, v2) | (-2) * u2 == ulpTowardZero = (u1 - ulpTowardZero, ulpTowardZero + u2)+ | otherwise = (u1, u2)+ r1 = scaleFloat e v1+ r2 = scaleFloat e v2+ in if isInfinite r1 then+ (r1, r1) -- unavoidable overflow+ else+ assert (r2 /= 0) (r1, r2)+{-# SPECIALIZE augmentedAddition :: Float -> Float -> (Float, Float), Double -> Double -> (Double, Double) #-}++-- |+-- IEEE 754 @augmentedSubtraction@ operation.+augmentedSubtraction :: RealFloat a => a -> a -> (a, a)+augmentedSubtraction x y = augmentedAddition x (negate y)++-- |+-- IEEE 754 @augmentedMultiplication@ operation.+augmentedMultiplication :: RealFloat a => a -> a -> (a, a)+augmentedMultiplication !x !y+ | isNaN x || isInfinite x || isNaN y || isInfinite y || x * y == 0 = let !result = x * y in (result, result)+ | otherwise = let exy = exponent x + exponent y+ x' = significand x+ y' = significand y+ (u1, u2) = twoProduct_nonscaling x' y'+ !_ = assert (toRational x' * toRational y' == toRational u1 + toRational u2) ()+ -- The product is subnormal <=> exy + exponent u1 < expMin+ -- The product is inexact => exy + exponent u1 < expMin + d+ in if exy + exponent u1 >= expMin then+ -- The result is exact+ let ulpTowardZero = u1 - nextTowardZero u1+ !_ = assert (2 * abs u2 <= abs ulpTowardZero) ()+ (v1, v2) = if (-2) * u2 == ulpTowardZero then+ (u1 - ulpTowardZero, ulpTowardZero + u2)+ else+ (u1, u2)+ !_ = assert (v1 + v2 == u1 + u2) ()+ r1 = scaleFloat exy v1+ -- !_ = assert (r1 == roundTiesTowardZero (fromRationalR (toRational x * toRational y))) ()+ in if isInfinite r1 then+ (r1, r1)+ else+ if v2 == 0 then+ (r1, 0 * r1) -- signed zero+ else+ if exy >= expMin + d then+ -- The result is exact+ let r2 = scaleFloat exy v2+ in (r1, r2)+ else+ -- The upper part is normal, the lower is subnormal (and inexact)+ -- Compute 'scaleFloat exy v2' with roundTiesTowardZero+ let !r2 = scaleFloatIntoSubnormalTiesTowardZero exy v2+ -- !_ = assert (r2 == roundTiesTowardZero (fromRationalR (toRational x * toRational y - toRational r1))) ()+ in (r1, r2)+ else+ -- The upper part is subnormal (possibly inexact), and the lower is signed zero (possibly inexact)+ if u2 == 0 then+ -- u1 is exact+ let !_ = assert (toRational x' * toRational y' == toRational u1) ()+ r1 = scaleFloatIntoSubnormalTiesTowardZero exy u1+ r1' = scaleFloat (-exy) r1+ in if u1 == r1' then+ (r1, 0 * r1)+ else+ (r1, 0 * (u1 - r1'))+ else+ let u1' = scaleFloat exy u1+ v1' = scaleFloat exy (if u2 > 0 then nextUp u1 else nextDown u1)+ r1 = if u1' == v1' || not (isMantissaEven u1') then+ u1'+ else+ v1'+ r1' = scaleFloat (-exy) r1+ in (r1, 0 * (u1 - r1' + u2))+ where+ d = floatDigits x+ (expMin,_expMax) = floatRange x++ -- Compute 'scaleFloat e z' with roundTiesTowardZero+ scaleFloatIntoSubnormalTiesTowardZero e z =+ let z' = scaleFloat e z+ w' = scaleFloat e (nextTowardZero z)+ in if z' == w' || not (isMantissaEven z') then+ z'+ else+ w'+{-# SPECIALIZE augmentedMultiplication :: Float -> Float -> (Float, Float), Double -> Double -> (Double, Double) #-}
+ src/Numeric/Floating/IEEE/Internal/Base.hs view
@@ -0,0 +1,136 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.Base+ ( isFloatBinary32+ , isDoubleBinary64+ , minPositive+ , minPositiveNormal+ , maxFinite+ , (^!)+ , negateIntAsWord+ , absIntAsWord+ ) where+import Data.Bits+import MyPrelude++default ()++-- $setup+-- >>> :set -XHexFloatLiterals -XNumericUnderscores+-- >>> import Numeric.Floating.IEEE.Internal.NextFloat (nextDown)++isFloatBinary32 :: Bool+isFloatBinary32 = isIEEE x+ && floatRadix x == 2+ && floatDigits x == 24+ && floatRange x == (-125, 128)+ where x :: Float+ x = undefined++isDoubleBinary64 :: Bool+isDoubleBinary64 = isIEEE x+ && floatRadix x == 2+ && floatDigits x == 53+ && floatRange x == (-1021, 1024)+ where x :: Double+ x = undefined++-- |+-- The smallest positive value expressible in an IEEE floating-point format.+-- This value is subnormal.+--+-- >>> (minPositive :: Float) == 0x1p-149+-- True+-- >>> (minPositive :: Double) == 0x1p-1074+-- True+-- >>> nextDown (minPositive :: Float)+-- 0.0+-- >>> nextDown (minPositive :: Double)+-- 0.0+minPositive :: RealFloat a => a+minPositive = let d = floatDigits x+ (expMin,_expMax) = floatRange x+ x = encodeFloat 1 (expMin - d)+ in x+{-# INLINABLE minPositive #-}+{-# SPECIALIZE minPositive :: Float, Double #-}++-- |+-- The smallest positive normal value expressible in an IEEE floating-point format.+--+-- >>> (minPositiveNormal :: Float) == 0x1p-126+-- True+-- >>> (minPositiveNormal :: Double) == 0x1p-1022+-- True+-- >>> isDenormalized (minPositiveNormal :: Float)+-- False+-- >>> isDenormalized (minPositiveNormal :: Double)+-- False+-- >>> isDenormalized (nextDown (minPositiveNormal :: Float))+-- True+-- >>> isDenormalized (nextDown (minPositiveNormal :: Double))+-- True+minPositiveNormal :: RealFloat a => a+minPositiveNormal = let (expMin,_expMax) = floatRange x+ x = encodeFloat 1 (expMin - 1)+ in x+{-# INLINABLE minPositiveNormal #-}+{-# SPECIALIZE minPositiveNormal :: Float, Double #-}++-- |+-- The largest finite value expressible in an IEEE floating-point format.+--+-- >>> (maxFinite :: Float) == 0x1.fffffep+127+-- True+-- >>> (maxFinite :: Double) == 0x1.ffff_ffff_ffff_fp+1023+-- True+maxFinite :: RealFloat a => a+maxFinite = let d = floatDigits x+ (_expMin,expMax) = floatRange x+ r = floatRadix x+ x = encodeFloat (r ^! d - 1) (expMax - d)+ in x+{-# INLINABLE maxFinite #-}+{-# SPECIALIZE maxFinite :: Float, Double #-}++-- A variant of (^) that allows constant folding+infixr 8 ^!+(^!) :: Integer -> Int -> Integer+(^!) = (^)+{-# INLINE [0] (^!) #-}++pow_helper :: Bool -> Integer -> Int -> Integer+pow_helper _ x y = x ^ y+{-# INLINE [0] pow_helper #-}+{-# RULES+"x^!" forall x y. x ^! y = pow_helper (y > 0) x y+"pow_helper/2" forall y.+ pow_helper True 2 y = bit y+"pow_helper" forall x y.+ pow_helper True x y = if y `rem` 2 == 0 then+ (x * x) ^! (y `quot` 2)+ else+ x * (x * x) ^! (y `quot` 2)+ #-}++-- |+-- >>> negateIntAsWord minBound == fromInteger (negate (fromIntegral (minBound :: Int)))+-- True+negateIntAsWord :: Int -> Word+negateIntAsWord x = fromIntegral (negate x)++-- |+-- >>> absIntAsWord minBound == fromInteger (abs (fromIntegral (minBound :: Int)))+-- True+absIntAsWord :: Int -> Word+absIntAsWord x = fromIntegral (abs x)++{- More careful definitions:++negateIntAsWord :: Int -> Word+negateIntAsWord x | x == minBound = fromInteger (negate (fromIntegral (minBound :: Int)))+ | otherwise = fromIntegral (negate x)++absIntAsWord :: Int -> Word+absIntAsWord x | x == minBound = fromInteger (abs (fromIntegral (minBound :: Int)))+ | otherwise = fromIntegral (abs x)+-}
+ src/Numeric/Floating/IEEE/Internal/Classify.hs view
@@ -0,0 +1,157 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE NumericUnderscores #-}+module Numeric.Floating.IEEE.Internal.Classify where+import Data.Bits+import GHC.Float.Compat (castDoubleToWord64, castFloatToWord32,+ isDoubleFinite, isFloatFinite)+import MyPrelude++default ()++-- |+-- IEEE 754 @isNormal@ operation.+isNormal :: RealFloat a => a -> Bool+isNormal x = x /= 0 && not (isNaN x) && not (isInfinite x) && not (isDenormalized x)+{-# NOINLINE [1] isNormal #-}+{-# RULES+"isNormal/Float" isNormal = isFloatNormal+"isNormal/Double" isNormal = isDoubleNormal+ #-}++isFloatNormal :: Float -> Bool+isFloatNormal x = let w = castFloatToWord32 x .&. 0x7f80_0000+ in w /= 0 && w /= 0x7f80_0000++isDoubleNormal :: Double -> Bool+isDoubleNormal x = let w = castDoubleToWord64 x .&. 0x7ff0_0000_0000_0000+ in w /= 0 && w /= 0x7ff0_0000_0000_0000++-- |+-- Returns @True@ if the argument is normal, subnormal, or zero.+--+-- IEEE 754 @isFinite@ operation.+isFinite :: RealFloat a => a -> Bool+isFinite x = not (isNaN x) && not (isInfinite x)+{-# NOINLINE [1] isFinite #-}+{-# RULES+"isFinite/Float"+ isFinite = \x -> isFloatFinite x /= 0+"isFinite/Double"+ isFinite = \x -> isDoubleFinite x /= 0+ #-}++-- |+-- Returns @True@ if the argument is zero.+--+-- IEEE 754 @isZero@ operation.+isZero :: RealFloat a => a -> Bool+isZero x = x == 0++-- |+-- Returns @True@ if the argument is negative (including negative zero).+--+-- Since 'RealFloat' constraint is insufficient to query the sign of NaNs,+-- this function treats all NaNs as positive.+-- See also "Numeric.Floating.IEEE.NaN".+--+-- IEEE 754 @isSignMinus@ operation.+isSignMinus :: RealFloat a => a -> Bool+isSignMinus x = x < 0 || isNegativeZero x++-- |+-- Comparison with IEEE 754 @totalOrder@ predicate.+--+-- Since 'RealFloat' constraint is insufficient to query the sign and payload of NaNs,+-- this function treats all NaNs as positive and does not make distinction between them.+-- See also "Numeric.Floating.IEEE.NaN".+--+-- Floating-point numbers are ordered as,+-- \(-\infty < \text{negative reals} < -0 < +0 < \text{positive reals} < +\infty < \mathrm{NaN}\).+compareByTotalOrder :: RealFloat a => a -> a -> Ordering+compareByTotalOrder x y+ | x < y = LT+ | y < x = GT+ | x == y = if x == 0 then+ compare (isNegativeZero y) (isNegativeZero x)+ else+ EQ+ | otherwise = compare (isNaN x) (isNaN y) -- The sign bit and payload of NaNs are ignored+-- TODO: Specialize for Float, Double++-- |+-- Comparison with IEEE 754 @totalOrderMag@ predicate.+--+-- Equivalent to @'compareByTotalOrder' (abs x) (abs y)@.+compareByTotalOrderMag :: RealFloat a => a -> a -> Ordering+compareByTotalOrderMag x y = compareByTotalOrder (abs x) (abs y)++-- isCanonical :: a -> Bool++-- data PartialOrdering = LT | EQ | GT | UNORD++-- |+-- The classification of floating-point values.+data Class = SignalingNaN+ | QuietNaN+ | NegativeInfinity+ | NegativeNormal+ | NegativeSubnormal+ | NegativeZero+ | PositiveZero+ | PositiveSubnormal+ | PositiveNormal+ | PositiveInfinity+ deriving (Eq, Ord, Show, Read, Enum)++-- |+-- Classifies a floating-point value.+--+-- Since 'RealFloat' constraint is insufficient to query signaling status of a NaN, this function treats all NaNs as quiet.+-- See also "Numeric.Floating.IEEE.NaN".+classify :: RealFloat a => a -> Class+classify x | isNaN x = QuietNaN+ | x < 0, isInfinite x = NegativeInfinity+ | x < 0, isDenormalized x = NegativeSubnormal+ | x < 0 = NegativeNormal+ | isNegativeZero x = NegativeZero+ | x == 0 = PositiveZero+ | isDenormalized x = PositiveSubnormal+ | isInfinite x = PositiveInfinity+ | otherwise = PositiveNormal+{-# NOINLINE [1] classify #-}+{-# RULES+"classify/Float" classify = classifyFloat+"classify/Double" classify = classifyDouble+ #-}++classifyFloat :: Float -> Class+classifyFloat x = let w = castFloatToWord32 x+ s = testBit w 31 -- sign bit+ e = (w `unsafeShiftR` 23) .&. 0xff -- exponent (8 bits)+ m = w .&. 0x007f_ffff -- mantissa (23 bits without leading 1)+ in case (s, e, m) of+ (True, 0, 0) -> NegativeZero+ (False, 0, 0) -> PositiveZero+ (True, 0, _) -> NegativeSubnormal+ (False, 0, _) -> PositiveSubnormal+ (True, 0xff, 0) -> NegativeInfinity+ (False, 0xff, 0) -> PositiveInfinity+ (_, 0xff, _) -> QuietNaN -- treat all NaNs as quiet+ (True, _, _) -> NegativeNormal+ (False, _, _) -> PositiveNormal++classifyDouble :: Double -> Class+classifyDouble x = let w = castDoubleToWord64 x+ s = testBit w 63 -- sign bit+ e = (w `unsafeShiftR` 52) .&. 0x7ff -- exponent (11 bits)+ m = w .&. 0x000f_ffff_ffff_ffff -- mantissa (52 bits without leading 1)+ in case (s, e, m) of+ (True, 0, 0) -> NegativeZero+ (False, 0, 0) -> PositiveZero+ (True, 0, _) -> NegativeSubnormal+ (False, 0, _) -> PositiveSubnormal+ (True, 0x7ff, 0) -> NegativeInfinity+ (False, 0x7ff, 0) -> PositiveInfinity+ (_, 0x7ff, _) -> QuietNaN -- treat all NaNs as quiet+ (True, _, _) -> NegativeNormal+ (False, _, _) -> PositiveNormal
+ src/Numeric/Floating/IEEE/Internal/Conversion.hs view
@@ -0,0 +1,68 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.Conversion+ ( realFloatToFrac+ , canonicalize+ , canonicalizeFloat+ , canonicalizeDouble+ ) where+import GHC.Float.Compat (double2Float, float2Double)+import MyPrelude++default ()++-- |+-- Converts a floating-point value into another type.+--+-- Similar to 'realToFrac', but treats NaN, infinities, negative zero even if the rewrite rule is off.+--+-- IEEE 754 @convertFormat@ operation.+realFloatToFrac :: (RealFloat a, Fractional b) => a -> b+realFloatToFrac x | isNaN x = 0/0+ | isInfinite x = if x > 0 then 1/0 else -1/0+ | isNegativeZero x = -0+ | otherwise = realToFrac x+{-# NOINLINE [1] realFloatToFrac #-}+{-# RULES+"realFloatToFrac/a->a" realFloatToFrac = canonicalize+"realFloatToFrac/Float->Double" realFloatToFrac = float2Double+"realFloatToFrac/Double->Float" realFloatToFrac = double2Float+ #-}++-- Since GHC optimizes away '* 1.0' when the type is 'Float' or 'Double',+-- we can't canonicalize x by just 'x * 1.0'.+one :: Num a => a+one = 1+{-# NOINLINE one #-}++-- |+-- A specialized version of 'realFloatToFrac'.+--+-- The resulting value will be canonical and non-signaling.+canonicalize :: RealFloat a => a -> a+canonicalize x = x * one+{-# INLINE [1] canonicalize #-}++#if defined(HAS_FAST_CANONICALIZE)++foreign import ccall unsafe "hs_canonicalizeFloat"+ canonicalizeFloat :: Float -> Float+foreign import ccall unsafe "hs_canonicalizeDouble"+ canonicalizeDouble :: Double -> Double++{-# RULES+"canonicalize/Float" canonicalize = canonicalizeFloat+"canonicalize/Double" canonicalize = canonicalizeDouble+ #-}++#else++{-# SPECIALIZE canonicalize :: Float -> Float, Double -> Double #-}++canonicalizeFloat :: Float -> Float+canonicalizeFloat = canonicalize++canonicalizeDouble :: Double -> Double+canonicalizeDouble = canonicalize++#endif
+ src/Numeric/Floating/IEEE/Internal/FMA.hs view
@@ -0,0 +1,301 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.FMA+ ( isMantissaEven+ , twoSum+ , addToOdd+ , split+ , twoProductFloat_viaDouble+ , twoProduct+ , twoProduct_nonscaling+ , twoProductFloat+ , twoProductDouble+ , fusedMultiplyAddFloat_viaDouble+ , fusedMultiplyAdd+ , fusedMultiplyAddFloat+ , fusedMultiplyAddDouble+ ) where+import Control.Exception (assert)+import Data.Bits+import GHC.Float.Compat (castDoubleToWord64, castFloatToWord32,+ double2Float, float2Double)+import MyPrelude+import Numeric.Floating.IEEE.Internal.Base (isDoubleBinary64,+ isFloatBinary32, (^!))+import Numeric.Floating.IEEE.Internal.Classify (isFinite)+import Numeric.Floating.IEEE.Internal.NextFloat (nextDown, nextUp)++default ()++-- $setup+-- >>> :set -XScopedTypeVariables++-- Assumption: input is finite+isMantissaEven :: RealFloat a => a -> Bool+isMantissaEven 0 = True+isMantissaEven x = let !_ = assert (isFinite x) ()+ (m,n) = decodeFloat x+ d = floatDigits x+ !_ = assert (floatRadix x ^ (d - 1) <= abs m && abs m < floatRadix x ^ d) ()+ (expMin, _expMax) = floatRange x+ s = expMin - (n + d)+ !_ = assert (isDenormalized x == (s > 0)) ()+ in if s > 0 then+ even (m `shiftR` s)+ else+ even m+{-# NOINLINE [1] isMantissaEven #-}+{-# RULES+"isMantissaEven/Double"+ isMantissaEven = \x -> even (castDoubleToWord64 x)+"isMantissaEven/Float"+ isMantissaEven = \x -> even (castFloatToWord32 x)+ #-}++-- |+-- Returns @x := a + b@ and @x - \<the exact value of (a + b)\>@.+--+-- This function does not avoid undue overflow;+-- For example, the second component of+-- @twoSum (0x1.017bd555b0b1fp1022) (-0x1.fffffffffffffp1023)@+-- is a NaN.+--+-- prop> \(a :: Double) (b :: Double) -> let (_,expMax) = floatRange a in max (exponent a) (exponent b) < expMax ==> let (x, y) = twoSum a b in a + b == x && toRational a + toRational b == toRational x + toRational y+twoSum :: RealFloat a => a -> a -> (a, a)+twoSum a b =+ let x = a + b+ t = x - a+ y = (a - (x - t)) + (b - t)+ {-+ Alternative:+ y = if abs b <= abs a then+ b - (x - a)+ else+ a - (x - b)+ -}+ in (x, y)+{-# SPECIALIZE twoSum :: Float -> Float -> (Float, Float), Double -> Double -> (Double, Double) #-}++-- |+-- Addition, with round to nearest odd floating-point number.+-- Like 'twoSum', this function does not handle undue overflow.+addToOdd :: RealFloat a => a -> a -> a+addToOdd x y = let (u, v) = twoSum x y+ result | isMantissaEven u && v < 0 = nextDown u+ | isMantissaEven u && v > 0 = nextUp u+ | isMantissaEven u && isNaN v && not (isInfinite u) =+ let v' = if abs y <= abs x then+ y - (u - x)+ else+ x - (u - y)+ in if v' < 0 then+ nextDown u+ else if v' > 0 then+ nextUp u+ else+ u+ | otherwise = u+ !_ = assert (isInfinite u || toRational u == toRational x + toRational y || not (isMantissaEven result)) ()+ in result+{-# SPECIALIZE addToOdd :: Float -> Float -> Float, Double -> Double -> Double #-}++-- This function doesn't handle overflow or underflow+split :: RealFloat a => a -> (a, a)+split a =+ let c = factor * a+ x = c - (c - a)+ y = a - x+ in (x, y)+ where factor = fromInteger $ 1 + floatRadix a ^! ((floatDigits a + 1) `quot` 2)+ -- factor == 134217729 for Double, 4097 for Float+{-# SPECIALIZE split :: Float -> (Float, Float), Double -> (Double, Double) #-}++-- This function will be rewritten into fastTwoProduct{Float,Double} if fast FMA is available; the rewriting may change behavior regarding overflow.+-- TODO: subnormal behavior?+-- |+-- prop> \(a :: Double) (b :: Double) -> let (x, y) = twoProduct a b in a * b == x && fromRational (toRational a * toRational b - toRational x) == y+twoProduct :: RealFloat a => a -> a -> (a, a)+twoProduct a b =+ let eab = exponent a + exponent b+ a' = significand a+ b' = significand b+ (ah, al) = split a'+ (bh, bl) = split b'+ x = a * b -- Since 'significand' doesn't honor the sign of zero, we can't use @a' * b'@+ y' = al * bl - (scaleFloat (-eab) x - ah * bh - al * bh - ah * bl)+ in (x, scaleFloat eab y')+{-# INLINABLE [1] twoProduct #-}++twoProductFloat_viaDouble :: Float -> Float -> (Float, Float)+twoProductFloat_viaDouble a b =+ let x, y :: Float+ a', b', x' :: Double+ a' = float2Double a+ b' = float2Double b+ x' = a' * b'+ x = double2Float x'+ y = double2Float (x' - float2Double x)+ in (x, y)++-- This function will be rewritten into fastTwoProduct{Float,Double} if fast FMA is available; the rewriting may change behavior regarding overflow.+twoProduct_nonscaling :: RealFloat a => a -> a -> (a, a)+twoProduct_nonscaling a b =+ let (ah, al) = split a+ (bh, bl) = split b+ x = a * b+ y = al * bl - (x - ah * bh - al * bh - ah * bl)+ in (x, y)+{-# NOINLINE [1] twoProduct_nonscaling #-}++twoProductFloat :: Float -> Float -> (Float, Float)+twoProductDouble :: Double -> Double -> (Double, Double)++#if defined(HAS_FAST_FMA)++twoProductFloat x y = let !r = x * y+ !s = fusedMultiplyAddFloat x y (-r)+ in (r, s)++twoProductDouble x y = let !r = x * y+ !s = fusedMultiplyAddDouble x y (-r)+ in (r, s)++{-# RULES+"twoProduct/Float" twoProduct = twoProductFloat+"twoProduct/Double" twoProduct = twoProductDouble+"twoProduct_nonscaling/Float" twoProduct_nonscaling = twoProductFloat+"twoProduct_nonscaling/Double" twoProduct_nonscaling = twoProductDouble+ #-}++#else++twoProductFloat = twoProductFloat_viaDouble+{-# INLINE twoProductFloat #-}++twoProductDouble = twoProduct+{-# INLINE twoProductDouble #-}++{-# RULES+"twoProduct/Float" twoProduct = twoProductFloat_viaDouble+"twoProduct_nonscaling/Float" twoProduct_nonscaling = twoProductFloat_viaDouble+ #-}+{-# SPECIALIZE twoProduct :: Double -> Double -> (Double, Double) #-}+{-# SPECIALIZE twoProduct_nonscaling :: Double -> Double -> (Double, Double) #-}++#endif++-- |+-- @'fusedMultiplyAdd' a b c@ computes @a * b + c@ as a single, ternary operation.+-- Rounding is done only once.+--+-- May make use of hardware FMA instructions if the target architecture has it; set @fma3@ package flag on x86 systems.+--+-- IEEE 754 @fusedMultiplyAdd@ operation.+--+-- prop> \(a :: Double) (b :: Double) (c :: Double) -> fusedMultiplyAdd a b c == fromRational (toRational a * toRational b + toRational c)+fusedMultiplyAdd :: RealFloat a => a -> a -> a -> a+fusedMultiplyAdd a b c+ | isFinite a && isFinite b && isFinite c =+ let eab | a == 0 || b == 0 = fst (floatRange a) - floatDigits a -- reasonably small+ | otherwise = exponent a + exponent b+ ec | c == 0 = fst (floatRange c) - floatDigits c+ | otherwise = exponent c++ -- Avoid overflow in twoProduct+ a' = significand a+ b' = significand b+ (x', y') = twoProduct_nonscaling a' b'+ !_ = assert (toRational a' * toRational b' == toRational x' + toRational y') ()++ -- Avoid overflow in twoSum+ e = max eab ec+ x = scaleFloat (eab - e) x'+ y = scaleFloat (eab - e) y'+ c'' = scaleFloat (max (fst (floatRange c) - floatDigits c + 1) (ec - e) - ec) c -- may be inexact++ (u1,u2) = twoSum y c''+ (v1,v2) = twoSum u1 x+ w = addToOdd u2 v2+ 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) ()+ in if result0 == 0 then+ -- We need to handle the sign of zero+ if c == 0 && a /= 0 && b /= 0 then+ a * b -- let a * b underflow+ else+ a * b + c -- -0 if both a * b and c are -0+ else+ if isDenormalized result then+ -- The rounding in 'scaleFloat e result0' may yield an incorrect result.+ -- Take the slow path.+ case toRational a * toRational b + toRational c of+ 0 -> a * b + c -- This should be exact+ r -> fromRational r+ else+ result+ | isFinite a && isFinite b = c + c -- c is +-Infinity or NaN+ | otherwise = a * b + c -- Infinity or NaN+{-# INLINABLE [1] fusedMultiplyAdd #-} -- May be rewritten into a more efficient one++fusedMultiplyAddFloat_viaDouble :: Float -> Float -> Float -> Float+fusedMultiplyAddFloat_viaDouble a b c+ | isFinite a && isFinite b && isFinite c =+ let a', b', c' :: Double+ a' = float2Double a+ b' = float2Double b+ c' = float2Double c+ ab = a' * b' -- exact+ !_ = assert (toRational ab == toRational a' * toRational b') ()+ result = double2Float (addToOdd ab c')+ !_ = assert (result == fromRational (toRational a * toRational b + toRational c)) ()+ in result+ | 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 defined(HAS_FAST_FMA)++foreign import ccall unsafe "hs_fusedMultiplyAddFloat"+ fusedMultiplyAddFloat :: Float -> Float -> Float -> Float+foreign import ccall unsafe "hs_fusedMultiplyAddDouble"+ fusedMultiplyAddDouble :: Double -> Double -> Double -> Double++{-# RULES+"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat+"fusedMultiplyAdd/Double" fusedMultiplyAdd = fusedMultiplyAddDouble+ #-}++#elif defined(USE_C99_FMA)++-- libm's fma might be implemented with hardware+foreign import ccall unsafe "fmaf"+ fusedMultiplyAddFloat :: Float -> Float -> Float -> Float+foreign import ccall unsafe "fma"+ fusedMultiplyAddDouble :: Double -> Double -> Double -> Double++{-# RULES+"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat+"fusedMultiplyAdd/Double" fusedMultiplyAdd = fusedMultiplyAddDouble+ #-}++#else++fusedMultiplyAddFloat :: Float -> Float -> Float -> Float+fusedMultiplyAddFloat = fusedMultiplyAddFloat_viaDouble+{-# INLINE fusedMultiplyAddFloat #-}++fusedMultiplyAddDouble :: Double -> Double -> Double -> Double+fusedMultiplyAddDouble = fusedMultiplyAdd -- generic implementation+{-# INLINE fusedMultiplyAddDouble #-}++{-# RULES+"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat_viaDouble+ #-}+{-# SPECIALIZE fusedMultiplyAdd :: Double -> Double -> Double -> Double #-}++#endif
+ src/Numeric/Floating/IEEE/Internal/Float128.hs view
@@ -0,0 +1,232 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE NumericUnderscores #-}+{-# OPTIONS_GHC -Wno-orphans -Wno-unused-imports #-}+module Numeric.Floating.IEEE.Internal.Float128 where+import Data.Bits+import Data.Word+import GHC.Exts (Int#)+import MyPrelude+import Numeric.Float128 (Float128 (F128))+import qualified Numeric.Float128+import Numeric.Floating.IEEE.Internal.Base+import Numeric.Floating.IEEE.Internal.Classify+import Numeric.Floating.IEEE.Internal.Conversion+import Numeric.Floating.IEEE.Internal.FMA+import Numeric.Floating.IEEE.Internal.NaN (RealFloatNaN)+import qualified Numeric.Floating.IEEE.Internal.NaN as NaN+import Numeric.Floating.IEEE.Internal.NextFloat+import Numeric.Floating.IEEE.Internal.Rounding+import Numeric.Floating.IEEE.Internal.RoundToIntegral++default ()++{-+Float128:+- exponent = 15 bits+- precision = 113 bits+-}++float128ToWord64Hi, float128ToWord64Lo :: Float128 -> Word64+float128ToWord64Hi (F128 hi _lo) = hi+float128ToWord64Lo (F128 _hi lo) = lo+{-# INLINE float128ToWord64Hi #-}+{-# INLINE float128ToWord64Lo #-}++float128ToWord64Pair :: Float128 -> (Word64, Word64)+float128ToWord64Pair (F128 hi lo) = (hi, lo)+{-# INLINE float128ToWord64Pair #-}++float128FromWord64Pair :: Word64 -- ^ higher 64 bits+ -> Word64 -- ^ lower 64 bits+ -> Float128+float128FromWord64Pair hi lo = F128 hi lo+{-# INLINE float128FromWord64Pair #-}++succWord64Pair :: Word64 -> Word64 -> (Word64, Word64)+succWord64Pair hi lo | lo + 1 == 0 = (hi + 1, 0)+ | otherwise = (hi, lo + 1)++predWord64Pair :: Word64 -> Word64 -> (Word64, Word64)+predWord64Pair hi lo | lo == 0 = (hi - 1, fromInteger (-1))+ | otherwise = (hi, lo - 1)++nextUpF128 :: Float128 -> Float128+nextUpF128 x =+ case float128ToWord64Pair x of+ (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000+ , (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+ case predWord64Pair hi lo of+ (hi', lo') -> float128FromWord64Pair hi' lo'+ | otherwise -> -- positive+ case succWord64Pair hi lo of+ (hi', lo') -> float128FromWord64Pair hi' lo'++nextDownF128 :: Float128 -> Float128+nextDownF128 x =+ case float128ToWord64Pair x of+ (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000+ , (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+ case succWord64Pair hi lo of+ (hi', lo') -> float128FromWord64Pair hi' lo'+ | otherwise -> -- positive+ case predWord64Pair hi lo of+ (hi', lo') -> float128FromWord64Pair hi' lo'++nextTowardZeroF128 :: Float128 -> Float128+nextTowardZeroF128 x =+ case float128ToWord64Pair x of+ (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000+ , (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+ (hi, lo) -> -- positive / negative+ case predWord64Pair hi lo of+ (hi', lo') -> float128FromWord64Pair hi' lo'++isNormalF128 :: Float128 -> Bool+isNormalF128 x = case float128ToWord64Pair x of+ (hi, _) -> let hi' = hi .&. 0x7fff_0000_0000_0000+ in hi' /= 0 && hi' /= 0x7fff_0000_0000_0000++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++classifyF128DiscardingSignalingNaNs :: Float128 -> Class+classifyF128DiscardingSignalingNaNs x =+ let hi = float128ToWord64Hi x+ s = testBit hi 63+ e = (hi `unsafeShiftR` 48) .&. 0x7fff -- exponent (15 bits)+ m_hi = hi .&. 0x0000_ffff_ffff_ffff+ m_lo = float128ToWord64Lo x+ in case (s, e, m_hi, m_lo) of+ (True, 0, 0, 0) -> NegativeZero+ (False, 0, 0, 0) -> PositiveZero+ (True, 0, _, _) -> NegativeSubnormal+ (False, 0, _, _) -> PositiveSubnormal+ (True, 0x7fff, 0, 0) -> NegativeInfinity+ (False, 0x7fff, 0, 0) -> PositiveInfinity+ (_, 0x7fff, _, _) -> QuietNaN -- treat all NaNs as quiet+ (True, _, _, _) -> NegativeNormal+ (False, _, _, _) -> PositiveNormal++instance RealFloatNaN Float128 where+ copySign x y = let (x_hi, x_lo) = float128ToWord64Pair x+ y_hi = float128ToWord64Hi y+ in float128FromWord64Pair ((x_hi .&. 0x7fff_ffff_ffff_ffff) .|. (y_hi .&. 0x8000_0000_0000_0000)) x_lo+ isSignMinus x = let hi = float128ToWord64Hi x+ in testBit hi 63+ isSignaling x = let hi = float128ToWord64Hi x+ in isNaN x && not (testBit hi 47)++ getPayload x+ | not (isNaN x) = -1+ | otherwise = let hi = fromIntegral (float128ToWord64Hi x .&. 0x0000_7fff_ffff_ffff)+ lo = fromIntegral (float128ToWord64Lo x)+ in hi * 0x1_0000_0000_0000_0000 + lo++ setPayload x+ | 0 <= x && x <= 0x0000_7fff_ffff_ffff_ffff_ffff_ffff_ffff+ = let payloadI = round x+ hi = fromInteger (payloadI `shiftR` 64) .|. 0x7fff_8000_0000_0000+ lo = fromInteger (payloadI .&. 0xffff_ffff_ffff_ffff)+ in float128FromWord64Pair hi lo+ | otherwise = 0++ setPayloadSignaling x+ | 0 < x && x <= 0x0000_7fff_ffff_ffff_ffff_ffff_ffff_ffff+ = let payloadI = round x+ hi = fromInteger (payloadI `shiftR` 64) .|. 0x7fff_0000_0000_0000+ lo = fromInteger (payloadI .&. 0xffff_ffff_ffff_ffff)+ in float128FromWord64Pair hi lo+ | otherwise = 0++ classify x =+ let hi = float128ToWord64Hi x+ s = testBit hi 63+ e = (hi `unsafeShiftR` 48) .&. 0x7fff -- exponent (15 bits)+ m_hi = hi .&. 0x0000_ffff_ffff_ffff+ m_lo = float128ToWord64Lo x+ in case (s, e, m_hi, m_lo) of+ (True, 0, 0, 0) -> NegativeZero+ (False, 0, 0, 0) -> PositiveZero+ (True, 0, _, _) -> NegativeSubnormal+ (False, 0, _, _) -> PositiveSubnormal+ (True, 0x7fff, 0, 0) -> NegativeInfinity+ (False, 0x7fff, 0, 0) -> PositiveInfinity+ (_, 0x7fff, _, _) -> if testBit m_hi 47 then+ QuietNaN+ else+ SignalingNaN+ (True, _, _, _) -> NegativeNormal+ (False, _, _, _) -> PositiveNormal++ compareByTotalOrder x y =+ let (x_hi, x_lo) = float128ToWord64Pair x+ (y_hi, y_lo) = float128ToWord64Pair y+ in compare (testBit y_hi 63) (testBit x_hi 63) -- sign bit+ <> if testBit x_hi 63 then+ compare y_hi x_hi <> compare y_lo x_lo -- negative+ else+ compare x_hi y_hi <> compare x_lo y_lo -- positive++{-# RULES+"nextUp/Float128" nextUp = nextUpF128+"nextDown/Float128" nextDown = nextDownF128+"nextTowardZero/Float128" nextTowardZero = nextTowardZeroF128+"isNormal/F128" isNormal = isNormalF128+"isFinite/F128" isFinite = isFiniteF128+"classify/F128" classify = classifyF128DiscardingSignalingNaNs+"isMantissaEven/F128"+ isMantissaEven = \x -> case x :: Float128 of F128 _hi lo -> even lo+"roundAway'/Float128" roundAway' = Numeric.Float128.round'+"ceiling'/Float128" ceiling' = Numeric.Float128.ceiling'+"floor'/Float128" floor' = Numeric.Float128.floor'+"truncate'/Float128" truncate' = Numeric.Float128.truncate'+ #-}++-- TODO: Write directly?+{-# SPECIALIZE minPositive :: Float128 #-}+{-# SPECIALIZE minPositiveNormal :: Float128 #-}+{-# SPECIALIZE maxFinite :: Float128 #-}++-- We shouldn't need specializations of positiveWordToBinaryFloatR# as long as WORD_SIZE_IN_BITS <= 113+{-# SPECIALIZE+ fromPositiveIntegerR :: RoundingStrategy f => Bool -> Integer -> f Float128+ , Bool -> Integer -> RoundTiesToEven Float128+ , Bool -> Integer -> RoundTiesToAway Float128+ , Bool -> Integer -> RoundTowardPositive Float128+ , Bool -> Integer -> RoundTowardNegative Float128+ , Bool -> Integer -> RoundTowardZero Float128+ #-}+{-# SPECIALIZE+ fromPositiveRatioR :: RoundingStrategy f => Bool -> Integer -> Integer -> f Float128+ , Bool -> Integer -> Integer -> RoundTiesToEven Float128+ , Bool -> Integer -> Integer -> RoundTiesToAway Float128+ , Bool -> Integer -> Integer -> RoundTowardPositive Float128+ , Bool -> Integer -> Integer -> RoundTowardNegative Float128+ , Bool -> Integer -> Integer -> RoundTowardZero Float128+ #-}+{-# SPECIALIZE+ encodePositiveFloatR# :: RoundingStrategy f => Bool -> Integer -> Int# -> f Float128+ , Bool -> Integer -> Int# -> RoundTiesToEven Float128+ , Bool -> Integer -> Int# -> RoundTiesToAway Float128+ , Bool -> Integer -> Int# -> RoundTowardPositive Float128+ , Bool -> Integer -> Int# -> RoundTowardNegative Float128+ , Bool -> Integer -> Int# -> RoundTowardZero Float128+ #-}+{-# SPECIALIZE+ scaleFloatR# :: RoundingStrategy f => Int# -> Float128 -> f Float128+ , Int# -> Float128 -> RoundTiesToEven Float128+ , Int# -> Float128 -> RoundTiesToAway Float128+ , Int# -> Float128 -> RoundTowardPositive Float128+ , Int# -> Float128 -> RoundTowardNegative Float128+ , Int# -> Float128 -> RoundTowardZero Float128+ #-}
+ src/Numeric/Floating/IEEE/Internal/GenericArith.hs view
@@ -0,0 +1,103 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.GenericArith where+import Data.Proxy+import MyPrelude+import Numeric.Floating.IEEE.Internal.Classify+import Numeric.Floating.IEEE.Internal.Conversion+import Numeric.Floating.IEEE.Internal.FMA++default ()++infixl 6 `genericAdd`, `genericSub`+infixl 7 `genericMul`, `genericDiv`++-- |+-- IEEE 754 @addition@ operation.+genericAdd :: (RealFloat a, RealFloat b) => a -> a -> b+genericAdd x y | x == 0 && y == 0 = realFloatToFrac (x + y)+ | isFinite x && isFinite y = fromRational (toRational x + toRational y)+ | otherwise = realFloatToFrac (x + y)+{-# NOINLINE [1] genericAdd #-}++-- |+-- IEEE 754 @subtraction@ operation.+genericSub :: (RealFloat a, RealFloat b) => a -> a -> b+genericSub x y | x == 0 && y == 0 = realFloatToFrac (x - y)+ | isFinite x && isFinite y = fromRational (toRational x - toRational y)+ | otherwise = realFloatToFrac (x - y)+{-# NOINLINE [1] genericSub #-}++-- |+-- IEEE 754 @multiplication@ operation.+genericMul :: (RealFloat a, RealFloat b) => a -> a -> b+genericMul x y | x == 0 || y == 0 = realFloatToFrac (x * y)+ | isFinite x && isFinite y = fromRational (toRational x * toRational y)+ | otherwise = realFloatToFrac (x * y)+{-# NOINLINE [1] genericMul #-}++-- |+-- IEEE 754 @division@ operation.+genericDiv :: (RealFloat a, RealFloat b) => a -> a -> b+genericDiv x y | x == 0 || y == 0 = realFloatToFrac (x / y)+ | isFinite x && isFinite y = fromRational (toRational x / toRational y)+ | otherwise = realFloatToFrac (x / y)+{-# NOINLINE [1] genericDiv #-}++{-+-- |+-- IEEE 754 @squareRoot@ operation.+genericSqrt :: (RealFloat a, RealFloat b) => a -> b+genericSqrt x | x == 0 = realFloatToFrac x+ | x > 0, isFinite x = error "not implemented yet"+ | otherwise = realFloatToFrac (sqrt x)+-}++-- |+-- IEEE 754 @fusedMultiplyAdd@ operation.+genericFusedMultiplyAdd :: (RealFloat a, RealFloat b) => a -> a -> a -> b+genericFusedMultiplyAdd a b c+ | isFinite a && isFinite b && isFinite c = case toRational a * toRational b + toRational c of+ 0 | isNegativeZero (a * b + c) -> -0+ r -> fromRational r+ | isFinite a && isFinite b = realFloatToFrac c -- c is Infinity or NaN+ | otherwise = realFloatToFrac (a * b + c)+{-# NOINLINE [1] genericFusedMultiplyAdd #-}++{-# RULES+"genericAdd/a->a" genericAdd = (+)+"genericSub/a->a" genericSub = (-)+"genericMul/a->a" genericMul = (*)+"genericDiv/a->a" genericDiv = (/)+"genericFusedMultiplyAdd/a->a" genericFusedMultiplyAdd = fusedMultiplyAdd+ #-}++-- | Returns True if @a@ is a subtype of @b@+--+-- >>> isSubFloatingType (undefined :: Float) (undefined :: Double)+-- True+-- >>> isSubFloatingType (undefined :: Double) (undefined :: Float)+-- False+-- >>> isSubFloatingType (undefined :: Double) (undefined :: Double)+-- True+isSubFloatingType :: (RealFloat a, RealFloat b) => a -> b -> Bool+isSubFloatingType a b = ieeeA && ieeeB && baseA == baseB && eminB <= eminA && emaxA <= emaxB && digitsA <= digitsB+ where+ ieeeA = isIEEE a+ ieeeB = isIEEE b+ baseA = floatRadix a+ baseB = floatRadix b+ (eminA,emaxA) = floatRange a+ (eminB,emaxB) = floatRange b+ digitsA = floatDigits a+ digitsB = floatDigits b++-- | Returns True if @a@ is a subtype of @b@+--+-- >>> isSubFloatingTypeProxy (Proxy :: Proxy Float) (Proxy :: Proxy Double)+-- True+-- >>> isSubFloatingTypeProxy (Proxy :: Proxy Double) (Proxy :: Proxy Float)+-- False+-- >>> isSubFloatingTypeProxy (Proxy :: Proxy Double) (Proxy :: Proxy Double)+-- True+isSubFloatingTypeProxy :: (RealFloat a, RealFloat b) => Proxy a -> Proxy b -> Bool+isSubFloatingTypeProxy proxyA proxyB = isSubFloatingType (undefined `asProxyTypeOf` proxyA) (undefined `asProxyTypeOf` proxyB)
+ src/Numeric/Floating/IEEE/Internal/Half.hs view
@@ -0,0 +1,254 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# OPTIONS_GHC -Wno-orphans -Wno-unused-imports #-}+module Numeric.Floating.IEEE.Internal.Half where+import Data.Bits+import Data.Coerce+import Data.Word+import Foreign.C.Types+import GHC.Exts+import GHC.Float.Compat (float2Double)+import MyPrelude+import Numeric.Floating.IEEE.Internal.Base+import Numeric.Floating.IEEE.Internal.Classify+import Numeric.Floating.IEEE.Internal.Conversion+import Numeric.Floating.IEEE.Internal.FMA+import Numeric.Floating.IEEE.Internal.NaN (RealFloatNaN)+import qualified Numeric.Floating.IEEE.Internal.NaN as NaN+import Numeric.Floating.IEEE.Internal.NextFloat+import Numeric.Floating.IEEE.Internal.Rounding+import Numeric.Half hiding (isZero)+import qualified Numeric.Half++default ()++castHalfToWord16 :: Half -> Word16+castHalfToWord16 (Half x) = coerce x+{-# INLINE castHalfToWord16 #-}++castWord16ToHalf :: Word16 -> Half+castWord16ToHalf x = Half (coerce x)+{-# INLINE castWord16ToHalf #-}++nextUpHalf :: Half -> Half+nextUpHalf x =+ case castHalfToWord16 x of+ w | w .&. 0x7c00 == 0x7c00+ , w /= 0xfc00 -> x + x -- NaN or negative infinity -> itself+ 0x8000 -> minPositive -- -0 -> min positive+ w | testBit w 15 -> castWord16ToHalf (w - 1) -- negative+ | otherwise -> castWord16ToHalf (w + 1) -- positive++nextDownHalf :: Half -> Half+nextDownHalf x =+ case castHalfToWord16 x of+ w | w .&. 0x7c00 == 0x7c00+ , w /= 0x7c00 -> x + x -- NaN or positive infinity -> itself+ 0x0000 -> - minPositive -- +0 -> max negative+ w | testBit w 15 -> castWord16ToHalf (w + 1) -- negative+ | otherwise -> castWord16ToHalf (w - 1) -- positive++nextTowardZeroHalf :: Half -> Half+nextTowardZeroHalf x =+ case castHalfToWord16 x of+ w | w .&. 0x7c00 == 0x7c00+ , w /= 0x7fff -> x + x -- NaN -> itself+ 0x8000 -> x -- -0 -> itself+ 0x0000 -> x -- +0 -> itself+ w -> castWord16ToHalf (w - 1) -- positive / negative++isNormalHalf :: Half -> Bool+isNormalHalf x = let w = castHalfToWord16 x .&. 0x7c00+ in w /= 0 && w /= 0x7c00++isFiniteHalf :: Half -> Bool+isFiniteHalf x = let w = castHalfToWord16 x .&. 0x7c00+ in w /= 0x7c00++isSignMinusHalf :: Half -> Bool+isSignMinusHalf x = let w = castHalfToWord16 x+ in testBit w 15 && (w .&. 0x7c00 /= 0x7c00 || w .&. 0x3ff == 0) -- all NaNs are treated as positive++classifyHalf :: Half -> Class+classifyHalf x = let w = castHalfToWord16 x+ s = testBit w 15+ e = (w `unsafeShiftR` 10) .&. 0x1f -- exponent (5 bits)+ m = w .&. 0x3ff -- mantissa (10 bits without leading 1)+ in case (s, e, m) of+ (True, 0, 0) -> NegativeZero+ (False, 0, 0) -> PositiveZero+ (True, 0, _) -> NegativeSubnormal+ (False, 0, _) -> PositiveSubnormal+ (True, 0x1f, 0) -> NegativeInfinity+ (False, 0x1f, 0) -> PositiveInfinity+ (_, 0x1f, _) -> QuietNaN -- treat all NaNs as quiet+ (True, _, _) -> NegativeNormal+ (False, _, _) -> PositiveNormal++instance RealFloatNaN Half where+ copySign x y = castWord16ToHalf ((x' .&. 0x7fff) .|. (y' .&. 0x8000))+ where x' = castHalfToWord16 x+ y' = castHalfToWord16 y++ isSignMinus x = testBit (castHalfToWord16 x) 15++ isSignaling x = x' .&. 0x7c00 == 0x7c00 && x' .&. 0x7fff /= 0x7c00 && not (testBit x' 9)+ where x' = castHalfToWord16 x++ getPayload x+ | not (isNaN x) = -1+ | otherwise = fromIntegral (castHalfToWord16 x .&. 0x01ff)++ setPayload x+ | 0 <= x && x <= 0x01ff = castWord16ToHalf $ round x .|. 0x7e00+ | otherwise = 0++ setPayloadSignaling x+ | 0 < x && x <= 0x01ff = castWord16ToHalf $ round x .|. 0x7c00+ | otherwise = 0++ classify x =+ let w = castHalfToWord16 x+ s = testBit w 15+ e = (w `unsafeShiftR` 10) .&. 0x1f -- exponent (5 bits)+ m = w .&. 0x3ff -- mantissa (10 bits without leading 1)+ in case (s, e, m) of+ (True, 0, 0) -> NegativeZero+ (False, 0, 0) -> PositiveZero+ (True, 0, _) -> NegativeSubnormal+ (False, 0, _) -> PositiveSubnormal+ (True, 0x1f, 0) -> NegativeInfinity+ (False, 0x1f, 0) -> PositiveInfinity+ (_, 0x1f, _) -> if testBit w 9 then+ QuietNaN+ else+ SignalingNaN+ (True, _, _) -> NegativeNormal+ (False, _, _) -> PositiveNormal++ equalByTotalOrder x y = castHalfToWord16 x == castHalfToWord16 y++ compareByTotalOrder x y =+ let x' = castHalfToWord16 x+ y' = castHalfToWord16 y+ in compare (testBit y' 15) (testBit x' 15) -- sign bit+ <> if testBit x' 15 then+ compare y' x' -- negative+ else+ compare x' y' -- positive++{-# RULES+"nextUp/Half" nextUp = nextUpHalf+"nextDown/Half" nextDown = nextDownHalf+"nextTowardZero/Half" nextTowardZero = nextTowardZeroHalf+"isNormal/Half" isNormal = isNormalHalf+"isFinite/Half" isFinite = isFiniteHalf+"isZero/Half" isZero = Numeric.Half.isZero+"isSignMinus/Half" isSignMinus = isSignMinusHalf+"classify/Half" classify = classifyHalf+"isMantissaEven/Half" forall (x :: Half).+ isMantissaEven x = even (castHalfToWord16 x)+ #-}++{-# SPECIALIZE minPositive :: Half #-}+{-# SPECIALIZE minPositiveNormal :: Half #-}+{-# SPECIALIZE maxFinite :: Half #-}+{-# SPECIALIZE+ positiveWordToBinaryFloatR# :: RoundingStrategy f => Bool -> Word# -> f Half+ , Bool -> Word# -> RoundTiesToEven Half+ , Bool -> Word# -> RoundTiesToAway Half+ , Bool -> Word# -> RoundTowardPositive Half+ , Bool -> Word# -> RoundTowardNegative Half+ , Bool -> Word# -> RoundTowardZero Half+ #-}+{-# SPECIALIZE+ fromPositiveIntegerR :: RoundingStrategy f => Bool -> Integer -> f Half+ , Bool -> Integer -> RoundTiesToEven Half+ , Bool -> Integer -> RoundTiesToAway Half+ , Bool -> Integer -> RoundTowardPositive Half+ , Bool -> Integer -> RoundTowardNegative Half+ , Bool -> Integer -> RoundTowardZero Half+ #-}+{-# SPECIALIZE+ fromPositiveRatioR :: RoundingStrategy f => Bool -> Integer -> Integer -> f Half+ , Bool -> Integer -> Integer -> RoundTiesToEven Half+ , Bool -> Integer -> Integer -> RoundTiesToAway Half+ , Bool -> Integer -> Integer -> RoundTowardPositive Half+ , Bool -> Integer -> Integer -> RoundTowardNegative Half+ , Bool -> Integer -> Integer -> RoundTowardZero Half+ #-}+{-# SPECIALIZE+ encodePositiveFloatR# :: RoundingStrategy f => Bool -> Integer -> Int# -> f Half+ , Bool -> Integer -> Int# -> RoundTiesToEven Half+ , Bool -> Integer -> Int# -> RoundTiesToAway Half+ , Bool -> Integer -> Int# -> RoundTowardPositive Half+ , Bool -> Integer -> Int# -> RoundTowardNegative Half+ , Bool -> Integer -> Int# -> RoundTowardZero Half+ #-}+{-# SPECIALIZE+ scaleFloatR# :: RoundingStrategy f => Int# -> Half -> f Half+ , Int# -> Half -> RoundTiesToEven Half+ , Int# -> Half -> RoundTiesToAway Half+ , Int# -> Half -> RoundTowardPositive Half+ , Int# -> Half -> RoundTowardNegative Half+ , Int# -> Half -> RoundTowardZero Half+ #-}++-- Monomorphic conversion functions+halfToFloat :: Half -> Float+halfToDouble :: Half -> Double+floatToHalf :: Float -> Half+doubleToHalf :: Double -> Half++#if defined(HAS_FAST_HALF_CONVERSION)++foreign import ccall unsafe "hs_fastHalfToFloat"+ c_fastHalfToFloat :: Word16 -> Float+foreign import ccall unsafe "hs_fastHalfToDouble"+ c_fastHalfToDouble :: Word16 -> Double+foreign import ccall unsafe "hs_fastFloatToHalf"+ c_fastFloatToHalf :: Float -> Word16+foreign import ccall unsafe "hs_fastDoubleToHalf"+ c_fastDoubleToHalf :: Double -> Word16++halfToFloat = coerce c_fastHalfToFloat+{-# INLINE halfToFloat #-}++halfToDouble = coerce c_fastHalfToDouble+{-# INLINE halfToDouble #-}++floatToHalf = coerce c_fastFloatToHalf+{-# INLINE floatToHalf #-}++doubleToHalf = coerce c_fastDoubleToHalf+{-# INLINE doubleToHalf #-}++{-# RULES+"realFloatToFrac/Half->Float" realFloatToFrac = halfToFloat+"realFloatToFrac/Half->Double" realFloatToFrac = halfToDouble+"realFloatToFrac/Float->Half" realFloatToFrac = floatToHalf+"realFloatToFrac/Double->Half" realFloatToFrac = doubleToHalf+ #-}++#else++halfToFloat = fromHalf+{-# INLINE halfToFloat #-}++halfToDouble = float2Double . fromHalf+{-# INLINE halfToDouble #-}++floatToHalf = toHalf+{-# INLINE floatToHalf #-}++doubleToHalf = realFloatToFrac -- generic implementation+{-# INLINE doubleToHalf #-}++{-# RULES+"realFloatToFrac/Half->Float" realFloatToFrac = fromHalf+"realFloatToFrac/Half->Double" realFloatToFrac = (realFloatToFrac . fromHalf) :: Half -> Double+"realFloatToFrac/Float->Half" realFloatToFrac = toHalf+ #-}++#endif
+ src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs view
@@ -0,0 +1,259 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE UnboxedSums #-}+{-# LANGUAGE UnboxedTuples #-}+{-# OPTIONS_GHC -Wno-unused-imports -fobject-code #-}++#include "MachDeps.h"++module Numeric.Floating.IEEE.Internal.IntegerInternals+ ( integerToIntMaybe+ , naturalToWordMaybe+ , unsafeShiftLInteger+ , unsafeShiftRInteger+ , roundingMode+ , countTrailingZerosInteger+ , integerIsPowerOf2+ , integerLog2IsPowerOf2+ ) where+import Data.Bits+import GHC.Exts (Int#, Word#, ctz#, int2Word#, plusWord#, quotRemInt#,+ uncheckedShiftL#, word2Int#, (+#), (-#))+import GHC.Int (Int (I#))+import GHC.Word (Word (W#))+import MyPrelude+import Numeric.Floating.IEEE.Internal.Base+import Numeric.Natural+#if defined(MIN_VERSION_ghc_bignum)+import qualified GHC.Num.BigNat+import GHC.Num.Integer (Integer (IN, IP, IS))+import qualified GHC.Num.Integer+import GHC.Num.Natural (Natural (NS))+#elif defined(MIN_VERSION_integer_gmp)+import qualified GHC.Integer+import GHC.Integer.GMP.Internals (Integer (Jn#, Jp#, S#),+ indexBigNat#)+import qualified GHC.Integer.Logarithms.Internals+import GHC.Natural (Natural (NatS#))+#define IN Jn#+#define IP Jp#+#define IS S#+#define NS NatS#+#else+import Math.NumberTheory.Logarithms (integerLog2')+#endif++-- $setup+-- >>> :m + Data.Int Test.QuickCheck+-- >>> :{+-- -- Workaround for https://github.com/sol/doctest/issues/160:+-- import Numeric.Floating.IEEE.Internal.IntegerInternals+-- :}++integerToIntMaybe :: Integer -> Maybe Int+naturalToWordMaybe :: Natural -> Maybe Word++-- The instance 'Bits Integer' is not very optimized...+unsafeShiftLInteger :: Integer -> Int -> Integer+unsafeShiftRInteger :: Integer -> Int -> Integer++-- |+-- Assumption: @n > 0@, @e >= 0@, and @integerLog2 n >= e@+--+-- Returns @compare (n \`'rem'\` 2^(e+1)) (2^e)@.+roundingMode :: Integer -- ^ @n@+ -> Int -- ^ @e@+ -> Ordering++-- |+-- 'Integer' version of 'countTrailingZeros'.+-- The argument must not be zero.+--+-- prop> \(NonZero x) -> countTrailingZerosInteger (toInteger x) === countTrailingZeros (x :: Int64)+-- >>> countTrailingZerosInteger 7+-- 0+-- >>> countTrailingZerosInteger 8+-- 3+countTrailingZerosInteger :: Integer -> Int++-- |+-- Returns @Just (integerLog2 x)@ if the argument @x@ is a power of 2, and @Nothing@ otherwise.+-- The argument @x@ must be strictly positive.+integerIsPowerOf2 :: Integer -> Maybe Int++-- |+-- Returns @(integerLog2 x, isJust (integerIsPowerOf2 x))@.+-- The argument @x@ must be strictly positive.+integerLog2IsPowerOf2 :: Integer -> (Int, Bool)++#if defined(MIN_VERSION_ghc_bignum) || defined(MIN_VERSION_integer_gmp)++integerToIntMaybe (IS x) = Just (I# x)+integerToIntMaybe _ = Nothing -- relies on Integer's invariant+{-# INLINE [0] integerToIntMaybe #-}++naturalToWordMaybe (NS x) = Just (W# x)+naturalToWordMaybe _ = Nothing -- relies on Natural's invariant+{-# INLINE [0] naturalToWordMaybe #-}++integerToIntMaybe2 :: Bool -> Integer -> Maybe Int+integerToIntMaybe2 _ (IS x) = Just (I# x)+integerToIntMaybe2 _ _ = Nothing+{-# INLINE [0] integerToIntMaybe2 #-}++naturalToWordMaybe2 :: Bool -> Natural -> Maybe Word+naturalToWordMaybe2 _ (NS x) = Just (W# x)+naturalToWordMaybe2 _ _ = Nothing+{-# INLINE [0] naturalToWordMaybe2 #-}++minBoundIntAsInteger :: Integer+minBoundIntAsInteger = fromIntegral (minBound :: Int)+{-# INLINE minBoundIntAsInteger #-}++maxBoundIntAsInteger :: Integer+maxBoundIntAsInteger = fromIntegral (maxBound :: Int)+{-# INLINE maxBoundIntAsInteger #-}++maxBoundWordAsNatural :: Natural+maxBoundWordAsNatural = fromIntegral (maxBound :: Word)+{-# INLINE maxBoundWordAsNatural #-}++{-# RULES+"integerToIntMaybe" [~0] forall x.+ integerToIntMaybe x = integerToIntMaybe2 (minBoundIntAsInteger <= x && x <= maxBoundIntAsInteger) x+"integerToIntMaybe2/small" forall x.+ integerToIntMaybe2 True x = Just (fromIntegral x)+"integerToIntMaybe2/large" forall x.+ integerToIntMaybe2 False x = Nothing+"naturalToWordMaybe" [~0] forall x.+ naturalToWordMaybe x = naturalToWordMaybe2 (x <= maxBoundWordAsNatural) x+"naturalToWordIntMaybe2/small" forall x.+ naturalToWordMaybe2 True x = Just (fromIntegral x)+"naturalToWordIntMaybe2/large" forall x.+ naturalToWordMaybe2 False x = Nothing+ #-}++#else++integerToIntMaybe = toIntegralSized+naturalToWordMaybe = toIntegralSized+{-# INLINE integerToIntMaybe #-}+{-# INLINE naturalToWordMaybe #-}++#endif++#if defined(MIN_VERSION_ghc_bignum)++unsafeShiftLInteger x (I# i) = GHC.Num.Integer.integerShiftL# x (int2Word# i)+unsafeShiftRInteger x (I# i) = GHC.Num.Integer.integerShiftR# x (int2Word# i)++#elif defined(MIN_VERSION_integer_gmp)++unsafeShiftLInteger x (I# i) = GHC.Integer.shiftLInteger x i+unsafeShiftRInteger x (I# i) = GHC.Integer.shiftRInteger x i++#else++unsafeShiftLInteger = unsafeShiftL+unsafeShiftRInteger = unsafeShiftR++#endif++{-# INLINE unsafeShiftLInteger #-}+{-# INLINE unsafeShiftRInteger #-}++#if defined(MIN_VERSION_ghc_bignum) || defined(MIN_VERSION_integer_gmp)++countTrailingZerosInteger# :: Integer -> Word#+countTrailingZerosInteger# (IS x) = ctz# (int2Word# x)+countTrailingZerosInteger# (IN bn) = countTrailingZerosInteger# (IP bn)+countTrailingZerosInteger# (IP bn) = loop 0# 0##+ where+ loop i acc =+ let+#if defined(MIN_VERSION_ghc_bignum)+ !bn_i = GHC.Num.BigNat.bigNatIndex# bn i -- `i < bigNatSize# bn` must hold+#else+ !bn_i = indexBigNat# bn i -- `i < sizeOfBigNat# bn` must hold+#endif+ in case bn_i of+ 0## -> loop (i +# 1#) (acc `plusWord#` WORD_SIZE_IN_BITS##)+ w -> acc `plusWord#` ctz# w++countTrailingZerosInteger 0 = error "countTrailingZerosInteger: zero"+countTrailingZerosInteger x = I# (word2Int# (countTrailingZerosInteger# x))+{-# INLINE countTrailingZerosInteger #-}++#else++countTrailingZerosInteger 0 = error "countTrailingZerosInteger: zero"+countTrailingZerosInteger x = integerLog2' (x `xor` (x - 1))+{-# INLINE countTrailingZerosInteger #-}++#endif++#if defined(MIN_VERSION_ghc_bignum)++roundingMode# :: Integer -> Int# -> Ordering+roundingMode# (IS x) t = let !w = int2Word# x+ in compare (W# (w `uncheckedShiftL#` (WORD_SIZE_IN_BITS# -# 1# -# t))) (W# (1## `uncheckedShiftL#` (WORD_SIZE_IN_BITS# -# 1#)))+roundingMode# (IN bn) t = roundingMode# (IP bn) t -- unexpected+roundingMode# (IP bn) t = case t `quotRemInt#` WORD_SIZE_IN_BITS# of+ -- 0 <= r < WORD_SIZE_IN_BITS+ (# s, r #) -> let !w = GHC.Num.BigNat.bigNatIndex# bn s+ -- w `shiftL` (WORD_SIZE_IN_BITS - r - 1) vs. 1 `shiftL` (WORD_SIZE_IN_BITS - 1)+ in compare (W# (w `uncheckedShiftL#` (WORD_SIZE_IN_BITS# -# 1# -# r))) (W# (1## `uncheckedShiftL#` (WORD_SIZE_IN_BITS# -# 1#)))+ <> loop s+ where+ loop 0# = EQ+ loop i = case GHC.Num.BigNat.bigNatIndex# bn i of+ 0## -> loop (i -# 1#)+ _ -> GT++roundingMode x (I# t) = roundingMode# x t+{-# INLINE roundingMode #-}++integerIsPowerOf2 x = case GHC.Num.Integer.integerIsPowerOf2# x of+ (# _ | #) -> Nothing+ (# | w #) -> Just (I# (word2Int# w))+{-# INLINE integerIsPowerOf2 #-}++integerLog2IsPowerOf2 x = case GHC.Num.Integer.integerIsPowerOf2# x of+ (# _ | #) -> (I# (word2Int# (GHC.Num.Integer.integerLog2# x)), False)+ (# | w #) -> (I# (word2Int# w), True)+{-# INLINE integerLog2IsPowerOf2 #-}++#elif defined(MIN_VERSION_integer_gmp)++roundingMode x (I# t#) = case GHC.Integer.Logarithms.Internals.roundingMode# x t# of+ 0# -> LT -- round toward zero+ 1# -> EQ -- half+ _ -> GT -- 2#: round away from zero+{-# INLINE roundingMode #-}++integerIsPowerOf2 x = case GHC.Integer.Logarithms.Internals.integerLog2IsPowerOf2# x of+ (# l, 0# #) -> Just (I# l)+ (# _, _ #) -> Nothing+{-# INLINE integerIsPowerOf2 #-}++integerLog2IsPowerOf2 x = case GHC.Integer.Logarithms.Internals.integerLog2IsPowerOf2# x of+ (# l, 0# #) -> (I# l, True)+ (# l, _ #) -> (I# l, False)+{-# INLINE integerLog2IsPowerOf2 #-}++#else++roundingMode x t = compare (x .&. (bit (t + 1) - 1)) (bit t)+{-# INLINE roundingMode #-}++integerIsPowerOf2 x = if x .&. (x - 1) == 0 then+ Just (integerLog2' x)+ else+ Nothing++integerLog2IsPowerOf2 x = (integerLog2' x, x .&. (x - 1) == 0)+{-# INLINE integerLog2IsPowerOf2 #-}++#endif
+ src/Numeric/Floating/IEEE/Internal/MinMax.hs view
@@ -0,0 +1,128 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.MinMax where+import MyPrelude++default ()++-- |+-- IEEE 754 @minimum@ operation.+-- @-0@ is smaller than @+0@.+-- Propagates NaNs.+minimum' :: RealFloat a => a -> a -> a+minimum' x y | isNaN x = x + x+ | isNaN y = y + y+ | x < y || (x == y && isNegativeZero x) = x+ | otherwise = y+{-# NOINLINE [1] minimum' #-}++-- |+-- IEEE 754 @minimumNumber@ operation.+-- @-0@ is smaller than @+0@.+-- Treats NaNs as missing data.+minimumNumber :: RealFloat a => a -> a -> a+minimumNumber x y | isNaN x && isNaN y = x + x+ | x < y || isNaN y || (x == y && isNegativeZero x) = x+ | otherwise = y+{-# NOINLINE [1] minimumNumber #-}++-- |+-- IEEE 754 @maximum@ operation.+-- @-0@ is smaller than @+0@.+-- Propagates NaNs.+maximum' :: RealFloat a => a -> a -> a+maximum' x y | isNaN x = x + x+ | isNaN y = y + y+ | x < y || (x == y && isNegativeZero x) = y+ | otherwise = x+{-# NOINLINE [1] maximum' #-}++-- |+-- IEEE 754 @maximumNumber@ operation.+-- @-0@ is smaller than @+0@.+-- Treats NaNs as missing data.+maximumNumber :: RealFloat a => a -> a -> a+maximumNumber x y | isNaN x && isNaN y = x + x+ | x < y || isNaN x || (x == y && isNegativeZero x) = y+ | otherwise = x+{-# NOINLINE [1] maximumNumber #-}++-- |+-- IEEE 754 @minimumMagnitude@ operation.+minimumMagnitude :: RealFloat a => a -> a -> a+minimumMagnitude x y | abs x < abs y = x+ | abs y < abs x = y+ | otherwise = minimum' x y++-- |+-- IEEE 754 @minimumMagnitudeNumber@ operation.+minimumMagnitudeNumber :: RealFloat a => a -> a -> a+minimumMagnitudeNumber x y | abs x < abs y = x+ | abs y < abs x = y+ | otherwise = minimumNumber x y++-- |+-- IEEE 754 @maximumMagnitude@ operation.+maximumMagnitude :: RealFloat a => a -> a -> a+maximumMagnitude x y | abs x > abs y = x+ | abs y > abs x = y+ | otherwise = maximum' x y++-- |+-- IEEE 754 @maximumMagnitudeNumber@ operation.+maximumMagnitudeNumber :: RealFloat a => a -> a -> a+maximumMagnitudeNumber x y | abs x > abs y = x+ | abs y > abs x = y+ | otherwise = maximumNumber x y++#if defined(HAS_FAST_MINMAX)++foreign import ccall unsafe "hs_minimumFloat"+ 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+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++#endif
+ src/Numeric/Floating/IEEE/Internal/NaN.hs view
@@ -0,0 +1,219 @@+{-# LANGUAGE NumericUnderscores #-}+module Numeric.Floating.IEEE.Internal.NaN+ ( module Numeric.Floating.IEEE.Internal.NaN+ , Class (..)+ ) where+import Data.Bits+import GHC.Float.Compat (castDoubleToWord64, castFloatToWord32,+ castWord32ToFloat, castWord64ToDouble)+import Numeric.Floating.IEEE.Internal.Classify (Class (..))++-- | An instance of this class supports manipulation of NaN.+class RealFloat a => RealFloatNaN a where+ {-# MINIMAL (copySign | isSignMinus), (isSignaling | classify), getPayload, setPayload, setPayloadSignaling #-}++ -- 5.5.1 Sign bit operations+ -- |+ -- Returns the first operand, with the sign of the second.+ --+ -- IEEE 754 @copySign@ operation.+ copySign :: a -> a -> a+ copySign x y = if isSignMinus x == isSignMinus y then+ x+ else+ -x++ -- 5.7.2 General operations+ -- |+ -- Returns @True@ if the operand is a negative number, negative infinity, negative zero, or a NaN with negative sign bit.+ --+ -- IEEE 754 @isSignMinus@ operation.+ isSignMinus :: a -> Bool+ isSignMinus x = copySign 1.0 x < 0++ -- |+ -- Returns @True@ if the operand is a signaling NaN.+ --+ -- IEEE 754 @isSignaling@ operation.+ isSignaling :: a -> Bool+ isSignaling x = classify x == SignalingNaN++ -- 9.7 NaN payload operations++ -- |+ -- Returns the payload of a NaN.+ -- Returns @-1@ if the operand is not a NaN.+ --+ -- IEEE 754 @getPayload@ operation.+ getPayload :: a -> a++ -- |+ -- Returns a quiet NaN with a given payload.+ -- Returns a positive zero if the payload is invalid.+ --+ -- IEEE 754 @setPayload@ operation.+ setPayload :: a -> a++ -- |+ -- Returns a signaling NaN with a given payload.+ -- Returns a positive zero if the payload is invalid.+ --+ -- IEEE 754 @setPayloadSignaling@ operation.+ setPayloadSignaling :: a -> a++ -- |+ -- IEEE 754 @class@ operation.+ classify :: a -> Class+ classify = classifyDefault++ -- |+ -- Equality with IEEE 754 @totalOrder@ operation.+ equalByTotalOrder :: a -> a -> Bool+ equalByTotalOrder x y = compareByTotalOrder x y == EQ++ -- |+ -- Comparison with IEEE 754 @totalOrder@ operation.+ compareByTotalOrder :: a -> a -> Ordering+ compareByTotalOrder = compareByTotalOrderDefault++classifyDefault :: RealFloatNaN a => a -> Class+classifyDefault x+ | isNaN x = if isSignaling x then+ SignalingNaN+ else+ QuietNaN+ | x < 0, isInfinite x = NegativeInfinity+ | x < 0, isDenormalized x = NegativeSubnormal+ | x < 0 = NegativeNormal+ | isNegativeZero x = NegativeZero+ | x == 0 = PositiveZero+ | isDenormalized x = PositiveSubnormal+ | isInfinite x = PositiveInfinity+ | otherwise = PositiveNormal++compareByTotalOrderDefault :: RealFloatNaN a => a -> a -> Ordering+compareByTotalOrderDefault x y+ | x < y = LT+ | y < x = GT+ | x == y = if x == 0 then+ compare (isNegativeZero y) (isNegativeZero x)+ else+ EQ -- TODO: non-canonical?+ | otherwise = compare (isSignMinus y) (isSignMinus x)+ <> let r = compare (isNaN x) (isNaN y) -- number < +NaN+ <> compare (isSignaling y) (isSignaling x) -- +(signaling NaN) < +(quiet NaN)+ <> compare (getPayload x) (getPayload y) -- implementation-defined+ in if isSignMinus x then+ compare EQ r+ else+ r++instance RealFloatNaN Float where+ copySign x y = castWord32ToFloat ((x' .&. 0x7fff_ffff) .|. (y' .&. 0x8000_0000))+ where x' = castFloatToWord32 x+ y' = castFloatToWord32 y++ isSignMinus x = testBit (castFloatToWord32 x) 31++ isSignaling x = x' .&. 0x7f80_0000 == 0x7f80_0000 && x' .&. 0x7fff_ffff /= 0x7f80_0000 && not (testBit x' 22)+ where x' = castFloatToWord32 x++ getPayload x+ | not (isNaN x) = -1+ | otherwise = fromIntegral (castFloatToWord32 x .&. 0x007f_ffff)++ setPayload x+ | 0 <= x && x <= 0x007f_ffff = castWord32ToFloat $ round x .|. 0x7fc0_0000+ | otherwise = 0++ setPayloadSignaling x+ | 0 < x && x <= 0x007f_ffff = castWord32ToFloat $ round x .|. 0x7f80_0000+ | otherwise = 0++ classify x = let w = castFloatToWord32 x+ s = testBit w 31 -- sign bit+ e = (w `unsafeShiftR` 23) .&. 0xff -- exponent (8 bits)+ m = w .&. 0x007f_ffff -- mantissa (23 bits without leading 1)+ in case (s, e, m) of+ (True, 0, 0) -> NegativeZero+ (False, 0, 0) -> PositiveZero+ (True, 0, _) -> NegativeSubnormal+ (False, 0, _) -> PositiveSubnormal+ (True, 0xff, 0) -> NegativeInfinity+ (False, 0xff, 0) -> PositiveInfinity+ (_, 0xff, _) -> if testBit w 22 then+ QuietNaN+ else+ SignalingNaN+ (True, _, _) -> NegativeNormal+ (False, _, _) -> PositiveNormal++ equalByTotalOrder x y = castFloatToWord32 x == castFloatToWord32 y++ compareByTotalOrder x y = let x' = castFloatToWord32 x+ y' = castFloatToWord32 y+ in compare (testBit y' 31) (testBit x' 31) -- sign bit+ <> if testBit x' 31 then+ compare y' x' -- negative+ else+ compare x' y' -- positive++instance RealFloatNaN Double where+ copySign x y = castWord64ToDouble ((x' .&. 0x7fff_ffff_ffff_ffff) .|. (y' .&. 0x8000_0000_0000_0000))+ where x' = castDoubleToWord64 x+ y' = castDoubleToWord64 y++ isSignMinus x = testBit (castDoubleToWord64 x) 63++ isSignaling x = x' .&. 0x7ff0_0000_0000_0000 == 0x7ff0_0000_0000_0000 && x' .&. 0x7fff_ffff_ffff_ffff /= 0x7ff0_0000_0000_0000 && not (testBit x' 51)+ where x' = castDoubleToWord64 x++ getPayload x+ | not (isNaN x) = -1+ | otherwise = fromIntegral (castDoubleToWord64 x .&. 0x0007_ffff_ffff_ffff)++ setPayload x+ | 0 <= x && x <= 0x0007_ffff_ffff_ffff = castWord64ToDouble $ round x .|. 0x7ff8_0000_0000_0000+ | otherwise = 0++ setPayloadSignaling x+ | 0 < x && x <= 0x0007_ffff_ffff_ffff = castWord64ToDouble $ round x .|. 0x7ff0_0000_0000_0000+ | otherwise = 0++ classify x = let w = castDoubleToWord64 x+ s = testBit w 63 -- sign bit+ e = (w `unsafeShiftR` 52) .&. 0x7ff -- exponent (11 bits)+ m = w .&. 0x000f_ffff_ffff_ffff -- mantissa (52 bits without leading 1)+ in case (s, e, m) of+ (True, 0, 0) -> NegativeZero+ (False, 0, 0) -> PositiveZero+ (True, 0, _) -> NegativeSubnormal+ (False, 0, _) -> PositiveSubnormal+ (True, 0x7ff, 0) -> NegativeInfinity+ (False, 0x7ff, 0) -> PositiveInfinity+ (_, 0x7ff, _) -> if testBit w 51 then+ QuietNaN+ else+ SignalingNaN+ (True, _, _) -> NegativeNormal+ (False, _, _) -> PositiveNormal++ equalByTotalOrder x y = castDoubleToWord64 x == castDoubleToWord64 y++ compareByTotalOrder x y = let x' = castDoubleToWord64 x+ y' = castDoubleToWord64 y+ in compare (testBit y' 63) (testBit x' 63) -- sign bit+ <> if testBit x' 63 then+ compare y' x' -- negative+ else+ compare x' y' -- positive++-- | A newtype wrapper to compare floating-point numbers by @totalOrder@ predicate.+newtype TotallyOrdered a = TotallyOrdered a+ deriving (Show)++instance RealFloatNaN a => Eq (TotallyOrdered a) where+ TotallyOrdered x == TotallyOrdered y = equalByTotalOrder x y++instance RealFloatNaN a => Ord (TotallyOrdered a) where+ compare (TotallyOrdered x) (TotallyOrdered y) = compareByTotalOrder x y
+ src/Numeric/Floating/IEEE/Internal/NextFloat.hs view
@@ -0,0 +1,280 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE NumericUnderscores #-}+{-# LANGUAGE BangPatterns #-}+module Numeric.Floating.IEEE.Internal.NextFloat where+import Data.Bits+import GHC.Float.Compat (castDoubleToWord64, castFloatToWord32,+ castWord32ToFloat, castWord64ToDouble)+import MyPrelude+import Numeric.Floating.IEEE.Internal.Base++default ()++-- $setup+-- >>> :set -XHexFloatLiterals -XNumericUnderscores++-- |+-- Returns the smallest value that is larger than the argument.+--+-- IEEE 754 @nextUp@ operation.+--+-- >>> nextUp 1 == (0x1.000002p0 :: Float)+-- True+-- >>> nextUp 1 == (0x1.0000_0000_0000_1p0 :: Double)+-- True+-- >>> nextUp (1/0) == (1/0 :: Double)+-- True+-- >>> nextUp (-1/0) == (- maxFinite :: Double)+-- True+-- >>> nextUp 0 == (0x1p-1074 :: Double)+-- True+-- >>> nextUp (-0) == (0x1p-1074 :: Double)+-- True+-- >>> nextUp (-0x1p-1074) :: Double -- returns negative zero+-- -0.0+nextUp :: RealFloat a => a -> a+nextUp x | not (isIEEE x) = error "non-IEEE numbers are not supported"+ | isNaN x || (isInfinite x && x > 0) = x + x -- NaN or positive infinity+ | x >= 0 = nextUp_positive x+ | otherwise = - nextDown_positive (- x)+{-# INLINE [1] nextUp #-}++-- |+-- Returns the largest value that is smaller than the argument.+--+-- IEEE 754 @nextDown@ operation.+--+-- >>> nextDown 1 == (0x1.ffff_ffff_ffff_fp-1 :: Double)+-- True+-- >>> nextDown 1 == (0x1.fffffep-1 :: Float)+-- True+-- >>> nextDown (1/0) == (maxFinite :: Double)+-- True+-- >>> nextDown (-1/0) == (-1/0 :: Double)+-- True+-- >>> nextDown 0 == (-0x1p-1074 :: Double)+-- True+-- >>> nextDown (-0) == (-0x1p-1074 :: Double)+-- True+-- >>> nextDown 0x1p-1074 -- returns positive zero+-- 0.0+-- >>> nextDown 0x1p-1022 == (0x0.ffff_ffff_ffff_fp-1022 :: Double)+-- True+nextDown :: RealFloat a => a -> a+nextDown x | not (isIEEE x) = error "non-IEEE numbers are not supported"+ | isNaN x || (isInfinite x && x < 0) = x + x -- NaN or negative infinity+ | x >= 0 = nextDown_positive x+ | otherwise = - nextUp_positive (- x)+{-# INLINE [1] nextDown #-}++-- |+-- Returns the value whose magnitude is smaller than that of the argument, and is closest to the argument.+--+-- This operation is not in IEEE, but may be useful to some.+--+-- >>> nextTowardZero 1 == (0x1.ffff_ffff_ffff_fp-1 :: Double)+-- True+-- >>> nextTowardZero 1 == (0x1.fffffep-1 :: Float)+-- True+-- >>> nextTowardZero (1/0) == (maxFinite :: Double)+-- True+-- >>> nextTowardZero (-1/0) == (-maxFinite :: Double)+-- True+-- >>> nextTowardZero 0 :: Double -- returns positive zero+-- 0.0+-- >>> nextTowardZero (-0 :: Double) -- returns negative zero+-- -0.0+-- >>> nextTowardZero 0x1p-1074 :: Double+-- 0.0+nextTowardZero :: RealFloat a => a -> a+nextTowardZero x | not (isIEEE x) = error "non-IEEE numbers are not supported"+ | isNaN x || x == 0 = x + x -- NaN or zero+ | x >= 0 = nextDown_positive x+ | otherwise = - nextDown_positive (- x)+{-# INLINE [1] nextTowardZero #-}++nextUp_positive :: RealFloat a => a -> a+nextUp_positive x+ | isNaN x || x < 0 = error "nextUp_positive"+ | isInfinite x = x+ | x == 0 = encodeFloat 1 (expMin - d) -- min positive+ | otherwise = let m :: Integer+ e :: Int+ (m,e) = decodeFloat x+ -- x = m * 2^e, 2^(d-1) <= m < 2^d+ -- 2^expMin < x < 2^expMax+ -- 2^(expMin-d): min positive+ -- 2^(expMin - 1): min normal 0x1p-1022+ -- expMin - d <= e <= expMax - d (-1074 .. 971)+ in if expMin - d <= e then+ -- normal+ if m + 1 == base ^! d && e == expMax - d then+ 1 / 0 -- max finite -> infinity+ else+ encodeFloat (m + 1) e+ else+ -- subnormal+ let m' = if base == 2 then+ m `unsafeShiftR` (expMin - d - e)+ else+ m `quot` (base ^ (expMin - d - e))+ in encodeFloat (m' + 1) (expMin - d)+ where+ d, expMin :: Int+ base = floatRadix x+ d = floatDigits x -- 53 for Double+ (expMin,expMax) = floatRange x -- (-1021,1024) for Double+{-# INLINE nextUp_positive #-}++nextDown_positive :: RealFloat a => a -> a+nextDown_positive x+ | isNaN x || x < 0 = error "nextDown_positive"+ | isInfinite x = maxFinite+ | x == 0 = encodeFloat (-1) (expMin - d) -- max negative+ | otherwise = let m :: Integer+ e :: Int+ (m,e) = decodeFloat x+ -- x = m * 2^e, 2^(d-1) <= m < 2^d+ -- 2^expMin < x < 2^expMax+ -- 2^(expMin-d): min positive+ -- 2^(expMin - 1): min normal 0x1p-1022+ -- expMin - d <= e <= expMax - d (-1074 .. 971)+ in if expMin - d <= e then+ -- normal+ let m1 = m - 1+ in if m == base ^! (d - 1) && expMin - d /= e then+ encodeFloat (base * m - 1) (e - 1)+ else+ encodeFloat m1 e+ else+ -- subnormal+ let m' = if base == 2 then+ m `unsafeShiftR` (expMin - d - e)+ else+ m `quot` (base ^ (expMin - d - e))+ in encodeFloat (m' - 1) (expMin - d)+ where+ d, expMin :: Int+ base = floatRadix x+ d = floatDigits x -- 53 for Double+ (expMin,_expMax) = floatRange x -- (-1021,1024) for Double+{-# INLINE nextDown_positive #-}++{-# RULES+"nextUp/Float" nextUp = nextUpFloat+"nextUp/Double" nextUp = nextUpDouble+"nextDown/Float" nextDown = nextDownFloat+"nextDown/Double" nextDown = nextDownDouble+"nextTowardZero/Float" nextTowardZero = nextTowardZeroFloat+"nextTowardZero/Double" nextTowardZero = nextTowardZeroDouble+ #-}++-- |+-- prop> nextUpFloat 1 == 0x1.000002p0+-- prop> nextUpFloat (1/0) == 1/0+-- prop> nextUpFloat (-1/0) == - maxFinite+-- prop> nextUpFloat 0 == 0x1p-149+-- prop> nextUpFloat (-0) == 0x1p-149+-- prop> isNegativeZero (nextUpFloat (-0x1p-149))+nextUpFloat :: Float -> Float+nextUpFloat x =+ case castFloatToWord32 x of+ w | w .&. 0x7f80_0000 == 0x7f80_0000+ , w /= 0xff80_0000 -> x + x -- NaN or positive infinity -> itself+ 0x8000_0000 -> minPositive -- -0 -> min positive+ 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"++-- |+-- prop> nextUpDouble 1 == 0x1.0000_0000_0000_1p0+-- prop> nextUpDouble (1/0) == 1/0+-- prop> nextUpDouble (-1/0) == - maxFinite+-- prop> nextUpDouble 0 == 0x1p-1074+-- prop> nextUpDouble (-0) == 0x1p-1074+-- prop> isNegativeZero (nextUpDouble (-0x1p-1074))+nextUpDouble :: Double -> Double+nextUpDouble x =+ case castDoubleToWord64 x of+ w | w .&. 0x7ff0_0000_0000_0000 == 0x7ff0_0000_0000_0000+ , w /= 0xfff0_0000_0000_0000 -> x + x -- NaN or positive infinity -> itself+ 0x8000_0000_0000_0000 -> minPositive -- -0 -> min positive+ 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"++-- |+-- prop> nextDownFloat 1 == 0x1.fffffep-1+-- prop> nextDownFloat (1/0) == maxFinite+-- prop> nextDownFloat (-1/0) == -1/0+-- prop> nextDownFloat 0 == -0x1p-149+-- prop> nextDownFloat (-0) == -0x1p-149+-- prop> nextDownFloat 0x1p-149 == 0+nextDownFloat :: Float -> Float+nextDownFloat x =+ case castFloatToWord32 x of+ w | w .&. 0x7f80_0000 == 0x7f80_0000+ , w /= 0x7f80_0000 -> x + x -- NaN or negative infinity -> itself+ 0x0000_0000 -> - minPositive -- +0 -> max negative+ 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"++-- |+-- prop> nextDownDouble 1 == 0x1.ffff_ffff_ffff_fp-1+-- prop> nextDownDouble (1/0) == maxFinite+-- prop> nextDownDouble (-1/0) == -1/0+-- prop> nextDownDouble 0 == -0x1p-1074+-- prop> nextDownDouble (-0) == -0x1p-1074+-- prop> nextDownDouble 0x1p-1074 == 0+nextDownDouble :: Double -> Double+nextDownDouble x =+ case castDoubleToWord64 x of+ w | w .&. 0x7ff0_0000_0000_0000 == 0x7ff0_0000_0000_0000+ , w /= 0x7ff0_0000_0000_0000 -> x + x -- NaN or negative infinity -> itself+ 0x0000_0000_0000_0000 -> - minPositive -- +0 -> max negative+ 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"++-- |+-- prop> nextTowardZeroFloat 1 == 0x1.fffffep-1+-- prop> nextTowardZeroFloat (-1) == -0x1.fffffep-1+-- prop> nextTowardZeroFloat (1/0) == maxFinite+-- prop> nextTowardZeroFloat (-1/0) == -maxFinite+-- prop> nextTowardZeroFloat 0 == 0+-- prop> isNegativeZero (nextTowardZeroFloat (-0))+-- prop> nextTowardZeroFloat 0x1p-149 == 0+nextTowardZeroFloat :: Float -> Float+nextTowardZeroFloat x =+ case castFloatToWord32 x of+ w | w .&. 0x7f80_0000 == 0x7f80_0000+ , w .&. 0x007f_ffff /= 0 -> x + x -- NaN -> itself+ 0x8000_0000 -> x -- -0 -> itself+ 0x0000_0000 -> x -- +0 -> itself+ w -> castWord32ToFloat (w - 1) -- positive / negative+ where+ !True = isFloatBinary32 || error "Numeric.Floating.Extra assumes Float is IEEE binary32"++-- |+-- prop> nextTowardZeroDouble 1 == 0x1.ffff_ffff_ffff_fp-1+-- prop> nextTowardZeroDouble (-1) == -0x1.ffff_ffff_ffff_fp-1+-- prop> nextTowardZeroDouble (1/0) == maxFinite+-- prop> nextTowardZeroDouble (-1/0) == -maxFinite+-- prop> nextTowardZeroDouble 0 == 0+-- prop> isNegativeZero (nextTowardZeroDouble (-0))+-- prop> nextTowardZeroDouble 0x1p-1074 == 0+nextTowardZeroDouble :: Double -> Double+nextTowardZeroDouble x =+ case castDoubleToWord64 x of+ w | w .&. 0x7ff0_0000_0000_0000 == 0x7ff0_0000_0000_0000+ , w .&. 0x000f_ffff_ffff_ffff /= 0 -> x + x -- NaN -> itself+ 0x8000_0000_0000_0000 -> x -- -0 -> itself+ 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"
+ src/Numeric/Floating/IEEE/Internal/Remainder.hs view
@@ -0,0 +1,35 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.Remainder+ ( remainder+ ) where+import MyPrelude+import Numeric.Floating.IEEE.Internal.Classify++default ()++-- |+-- @'remainder' x y@ returns \(r=x-yn\), where \(n\) is the integer nearest the exact number \(x/y\); i.e. \(n=\mathrm{round}(x/y)\).+--+-- IEEE 754 @remainder@ operation.+remainder :: RealFloat a => a -> a -> a+remainder x y | isFinite x && isInfinite y = x+ | 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+{-# NOINLINE [1] remainder #-}++#if defined(USE_FFI)++foreign import ccall unsafe "remainderf"+ c_remainderFloat :: Float -> Float -> Float+foreign import ccall unsafe "remainder"+ c_remainderDouble :: Double -> Double -> Double++{-# RULES+"remainder/Float" remainder = c_remainderFloat+"remainder/Double" remainder = c_remainderDouble+ #-}++#endif
+ src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs view
@@ -0,0 +1,193 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.RoundToIntegral+ ( round'+ , roundAway'+ , truncate'+ , ceiling'+ , floor'+ , round+ , roundAway+ , truncate+ , ceiling+ , floor+ ) where+import MyPrelude++default ()++-- $setup+-- >>> :set -XScopedTypeVariables+-- >>> import Numeric.Floating.IEEE.Internal.Classify (isFinite)++-- |+-- @'round'' x@ returns the nearest integral value to @x@; the even integer if @x@ is equidistant between two integers.+--+-- IEEE 754 @roundToIntegralTiesToEven@ operation.+--+-- prop> \(x :: Double) -> isFinite x ==> (round' x == fromInteger (round x))+-- >>> round' (-0.5)+-- -0.0+round' :: RealFloat a => a -> a+round' x | isInfinite x || isNaN x || isNegativeZero x = x + x+round' x = case round x of+ 0 | x < 0 -> -0+ | otherwise -> 0+ n -> fromInteger n+{-# NOINLINE [1] round' #-}++-- |+-- @'roundAway'' x@ returns the nearest integral value to @x@; the one with larger magnitude is returned if @x@ is equidistant between two integers.+--+-- IEEE 754 @roundToIntegralTiesToAway@ operation.+--+-- prop> \(x :: Double) -> isFinite x ==> roundAway' x == fromInteger (roundAway x)+-- >>> roundAway' (-0.5)+-- -1.0+-- >>> roundAway' (-0.4)+-- -0.0+roundAway' :: RealFloat a => a -> a+roundAway' x | isInfinite x || isNaN x || isNegativeZero x = x + x+roundAway' x = case properFraction x of+ -- x == n + f, signum x == signum f, 0 <= abs f < 1+ (n,r) -> if abs r < 0.5 then+ -- round toward zero+ if n == 0 then+ 0.0 * r -- signed zero+ else+ fromInteger n+ else+ -- round away from zero+ if r < 0 then+ fromInteger (n - 1)+ else+ fromInteger (n + 1)+{-# NOINLINE [1] roundAway' #-}++-- |+-- @'truncate'' x@ returns the integral value nearest to @x@, and whose magnitude is not greater than that of @x@.+--+-- IEEE 754 @roundToIntegralTowardZero@ operation.+--+-- prop> \(x :: Double) -> isFinite x ==> truncate' x == fromInteger (truncate x)+-- >>> truncate' (-0.5)+-- -0.0+truncate' :: RealFloat a => a -> a+truncate' x | isInfinite x || isNaN x || isNegativeZero x = x + x+truncate' x = case truncate x of+ 0 | x < 0 -> -0+ | otherwise -> 0+ n -> fromInteger n+{-# NOINLINE [1] truncate' #-}++-- |+-- @'ceiling'' x@ returns the least integral value that is not less than @x@.+--+-- IEEE 754 @roundToIntegralTowardPositive@ operation.+--+-- prop> \(x :: Double) -> isFinite x ==> ceiling' x == fromInteger (ceiling x)+-- >>> ceiling' (-0.8)+-- -0.0+-- >>> ceiling' (-0.5)+-- -0.0+ceiling' :: RealFloat a => a -> a+ceiling' x | isInfinite x || isNaN x || isNegativeZero x = x + x+ceiling' x = case ceiling x of+ 0 | x < 0 -> -0+ | otherwise -> 0+ n -> fromInteger n+{-# NOINLINE [1] ceiling' #-}++-- |+-- @'floor'' x@ returns the greatest integral value that is not greater than @x@.+--+-- IEEE 754 @roundToIntegralTowardNegative@ operation.+--+-- prop> \(x :: Double) -> isFinite x ==> floor' x == fromInteger (floor x)+-- >>> floor' (-0.1)+-- -1.0+-- >>> floor' (-0)+-- -0.0+floor' :: RealFloat a => a -> a+floor' x | isInfinite x || isNaN x || isNegativeZero x = x + x+ | otherwise = fromInteger (floor x)+{-# NOINLINE [1] floor' #-}++-- |+-- @'roundAway' x@ returns the nearest integer to @x@; the integer with larger magnitude is returned if @x@ is equidistant between two integers.+--+-- IEEE 754 @convertToIntegerTiesToAway@ operation.+--+-- >>> roundAway 4.5+-- 5+roundAway :: (RealFrac a, Integral b) => a -> b+roundAway x = case properFraction x of+ -- x == n + f, signum x == signum f, 0 <= abs f < 1+ (n,r) -> if abs r < 0.5 then+ n+ else+ if r < 0 then+ n - 1+ else+ n + 1+{-# INLINE roundAway #-}++#ifdef USE_FFI++foreign import ccall unsafe "ceilf"+ c_ceilFloat :: Float -> Float+foreign import ccall unsafe "ceil"+ c_ceilDouble :: Double -> Double+foreign import ccall unsafe "floorf"+ c_floorFloat :: Float -> Float+foreign import ccall unsafe "floor"+ c_floorDouble :: Double -> Double+foreign import ccall unsafe "roundf"+ c_roundFloat :: Float -> Float -- ties to away+foreign import ccall unsafe "round"+ c_roundDouble :: Double -> Double -- ties to away+foreign import ccall unsafe "truncf"+ c_truncFloat :: Float -> Float+foreign import ccall unsafe "trunc"+ c_truncDouble :: Double -> Double++{-# RULES+"roundAway'/Float"+ roundAway' = c_roundFloat+"roundAway'/Double"+ roundAway' = c_roundDouble+"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+ #-}++{- from base+foreign import ccall unsafe "rintFloat"+ c_rintFloat :: Float -> Float+foreign import ccall unsafe "rintDouble"+ c_rintDouble :: Double -> Double+-}+#if defined(HAS_FAST_ROUNDEVEN)+foreign import ccall unsafe "hs_roundevenFloat"+ c_roundevenFloat :: Float -> Float+foreign import ccall unsafe "hs_roundevenDouble"+ c_roundevenDouble :: Double -> Double++{-# RULES+"round'/Float"+ round' = c_roundevenFloat+"round'/Double"+ round' = c_roundevenDouble+ #-}+#endif++#endif
+ src/Numeric/Floating/IEEE/Internal/Rounding.hs view
@@ -0,0 +1,7 @@+module Numeric.Floating.IEEE.Internal.Rounding (module M) where+import Numeric.Floating.IEEE.Internal.Rounding.Common as M+import Numeric.Floating.IEEE.Internal.Rounding.Encode as M+import Numeric.Floating.IEEE.Internal.Rounding.Integral as M+import Numeric.Floating.IEEE.Internal.Rounding.Rational as M++default ()
+ src/Numeric/Floating/IEEE/Internal/Rounding/Common.hs view
@@ -0,0 +1,165 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE NoImplicitPrelude #-}+module Numeric.Floating.IEEE.Internal.Rounding.Common where+import Control.Exception (assert)+import Data.Bits+import Data.Functor.Product+import Data.Int+import GHC.Float (expt)+import Math.NumberTheory.Logarithms (integerLog2')+import MyPrelude+import Numeric.Floating.IEEE.Internal.IntegerInternals++default ()++class Functor f => RoundingStrategy f where+ exact :: a -> f a+ inexact :: Ordering -- ^ LT -> toward-zero is the nearest, EQ -> midpoint, GT -> away-from-zero is the nearest+ -> Bool -- ^ negative (True -> negative, False -> positive)+ -> Int -- ^ parity (even -> toward-zero is even, odd -> toward-zero is odd)+ -> a -- ^ toward zero+ -> a -- ^ away from zero+ -> f a+ doRound :: Bool -- ^ exactness; if True, the Ordering must be LT+ -> Ordering -- ^ LT -> toward-zero is the nearest, EQ -> midpoint, GT -> away-from-zero is the nearest+ -> Bool -- ^ negative (True -> negative, False -> positive)+ -> Int -- ^ parity (even -> toward-zero is even, odd -> toward-zero is odd)+ -> a -- ^ toward zero+ -> a -- ^ away from zero+ -> f a+ exact x = doRound True LT False 0 x x+ inexact o neg parity zero away = doRound False o neg parity zero away++newtype RoundTiesToEven a = RoundTiesToEven { roundTiesToEven :: a }+ deriving (Functor)++instance RoundingStrategy RoundTiesToEven where+ exact = RoundTiesToEven+ inexact o _neg parity zero away = RoundTiesToEven $ case o of+ LT -> zero+ EQ | even parity -> zero+ | otherwise -> away+ GT -> away+ doRound _ex o _neg parity zero away = RoundTiesToEven $ case o of+ LT -> zero+ EQ | even parity -> zero+ | otherwise -> away+ GT -> away+ {-# INLINE exact #-}+ {-# INLINE inexact #-}+ {-# INLINE doRound #-}++newtype RoundTiesToAway a = RoundTiesToAway { roundTiesToAway :: a }+ deriving (Functor)++instance RoundingStrategy RoundTiesToAway where+ exact = RoundTiesToAway+ inexact o _neg _parity zero away = RoundTiesToAway $ case o of+ LT -> zero+ EQ -> away+ GT -> away+ doRound _ex o _neg _parity zero away = RoundTiesToAway $ case o of+ LT -> zero+ EQ -> away+ GT -> away+ {-# INLINE exact #-}+ {-# INLINE inexact #-}+ {-# INLINE doRound #-}++newtype RoundTowardPositive a = RoundTowardPositive { roundTowardPositive :: a }+ deriving (Functor)++instance RoundingStrategy RoundTowardPositive where+ exact = RoundTowardPositive+ inexact _o neg _parity zero away | neg = RoundTowardPositive zero+ | otherwise = RoundTowardPositive away+ doRound ex _o neg _parity zero away | ex || neg = RoundTowardPositive zero+ | otherwise = RoundTowardPositive away+ {-# INLINE exact #-}+ {-# INLINE inexact #-}+ {-# INLINE doRound #-}++newtype RoundTowardNegative a = RoundTowardNegative { roundTowardNegative :: a }+ deriving (Functor)++instance RoundingStrategy RoundTowardNegative where+ exact = RoundTowardNegative+ inexact _o neg _parity zero away | neg = RoundTowardNegative away+ | otherwise = RoundTowardNegative zero+ doRound ex _o neg _parity zero away | not ex && neg = RoundTowardNegative away+ | otherwise = RoundTowardNegative zero+ {-# INLINE exact #-}+ {-# INLINE inexact #-}+ {-# INLINE doRound #-}++newtype RoundTowardZero a = RoundTowardZero { roundTowardZero :: a }+ deriving (Functor)++instance RoundingStrategy RoundTowardZero where+ exact = RoundTowardZero+ inexact _o _neg _parity zero _away = RoundTowardZero zero+ doRound _ex _o _neg _parity zero _away = RoundTowardZero zero+ {-# INLINE exact #-}+ {-# INLINE inexact #-}+ {-# INLINE doRound #-}++instance (RoundingStrategy f, RoundingStrategy g) => RoundingStrategy (Product f g) where+ exact x = Pair (exact x) (exact x)+ inexact o neg parity zero away = Pair (inexact o neg parity zero away) (inexact o neg parity zero away)+ doRound ex o neg parity zero away = Pair (doRound ex o neg parity zero away) (doRound ex o neg parity zero away)+ {-# INLINE exact #-}+ {-# INLINE inexact #-}+ {-# INLINE doRound #-}++{-+from GHC.Float:+expt :: Integer -> Int -> Integer+expt base n = base ^ n+-}++quotRemByExpt :: Integer -- ^ the dividend @x@+ -> Integer -- ^ base+ -> Int -- ^ the exponent @e@ (must be non-negative)+ -> (Integer, Integer) -- ^ @x \`'quotRem'\` (base ^ e)@+quotRemByExpt x 2 n = assert (n >= 0) (x `unsafeShiftRInteger` n, x .&. (bit n - 1))+quotRemByExpt x base n = x `quotRem` expt base n+{-# INLINE quotRemByExpt #-}++multiplyByExpt :: Integer -- ^ the multiplicand @x@+ -> Integer -- ^ base+ -> Int -- ^ the exponent @e@ (must be non-negative)+ -> Integer -- ^ @x * base ^ e@+multiplyByExpt x 2 n = assert (n >= 0) (x `unsafeShiftLInteger` n)+multiplyByExpt x base n = x * expt base n+{-# INLINE multiplyByExpt #-}++isDivisibleByExpt :: Integer -- ^ the dividend @x@+ -> Integer -- ^ the base+ -> Int -- ^ the exponent @e@ (must be non-negative)+ -> Integer -- ^ the remainder @r@ (must be @x \`'rem'\` (base ^ e)@)+ -> Bool -- ^ @r == 0@+isDivisibleByExpt x 2 e r = assert (r == x `rem` (2 ^ e)) $ x == 0 || Numeric.Floating.IEEE.Internal.IntegerInternals.countTrailingZerosInteger x >= e+isDivisibleByExpt x base e r = assert (r == x `rem` (base ^ e)) (r == 0)+{-# INLINE isDivisibleByExpt #-}++-- |+-- Assumption: @n >= 0@, @e >= 0@, and @r == n \`'rem'\` base^(e+1)@+--+-- Returns @compare r (base^e)@.+compareWithExpt :: Integer -- ^ base+ -> Integer -- ^ the number @n@ (must be non-negative)+ -> Integer -- ^ the remainder @r@ (must be @n \`'rem'\' base^(e+1)@)+ -> Int -- ^ the exponent @e@ (must be non-negative)+ -> Ordering+compareWithExpt 2 n r e = assert (r == n `rem` expt 2 (e+1)) $+ if n == 0 || integerLog2' n < e then+ -- If integerLog2 n < e (i.e. n < 2^e), it is trivial+ LT+ else+ -- In this branch, n > 0 && integerLog2' n >= e+ let result = Numeric.Floating.IEEE.Internal.IntegerInternals.roundingMode n e+ !_ = assert (result == compare r (expt 2 e)) ()+ in result+compareWithExpt base n r e = assert (r == n `rem` expt base (e+1)) $ compare r (expt base e)+{-# INLINE compareWithExpt #-}
+ src/Numeric/Floating/IEEE/Internal/Rounding/Encode.hs view
@@ -0,0 +1,204 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Numeric.Floating.IEEE.Internal.Rounding.Encode where+import Control.Exception (assert)+import Data.Functor.Product+import Data.Int+import GHC.Exts+import Math.NumberTheory.Logarithms (integerLog2', integerLogBase')+import MyPrelude+import Numeric.Floating.IEEE.Internal.Base+import Numeric.Floating.IEEE.Internal.Classify (isFinite)+import Numeric.Floating.IEEE.Internal.Rounding.Common++default ()++encodeFloatTiesToEven, encodeFloatTiesToAway, encodeFloatTowardPositive, encodeFloatTowardNegative, encodeFloatTowardZero :: RealFloat a => Integer -> Int -> a+encodeFloatTiesToEven m = roundTiesToEven . encodeFloatR m+encodeFloatTiesToAway m = roundTiesToAway . encodeFloatR m+encodeFloatTowardPositive m = roundTowardPositive . encodeFloatR m+encodeFloatTowardNegative m = roundTowardNegative . encodeFloatR m+encodeFloatTowardZero m = roundTowardZero . encodeFloatR m+{-# INLINE encodeFloatTiesToEven #-}+{-# INLINE encodeFloatTiesToAway #-}+{-# INLINE encodeFloatTowardPositive #-}+{-# INLINE encodeFloatTowardNegative #-}+{-# INLINE encodeFloatTowardZero #-}++encodeFloatR :: (RealFloat a, RoundingStrategy f) => Integer -> Int -> f a+encodeFloatR 0 !_ = exact 0+encodeFloatR m n | m < 0 = negate <$> encodePositiveFloatR True (- m) n+ | otherwise = encodePositiveFloatR False m n+{-# INLINE encodeFloatR #-}++-- Avoid cross-module specialization issue with manual worker/wrapper transformation+encodePositiveFloatR :: (RealFloat a, RoundingStrategy f) => Bool -> Integer -> Int -> f a+encodePositiveFloatR neg m (I# n#) = encodePositiveFloatR# neg m n#+{-# INLINE encodePositiveFloatR #-}++encodePositiveFloatR# :: forall f a. (RealFloat a, RoundingStrategy f) => Bool -> Integer -> Int# -> f a+encodePositiveFloatR# !neg !m n# = assert (m > 0) result+ where+ n = I# n#+ result = let k = if base == 2 then+ integerLog2' m+ else+ integerLogBase' base m+ -- base^k <= m < base^(k+1)+ -- base^^(k+n) <= m * base^^n < base^^(k+n+1)+ in if expMin <= k + n + 1 && k + n + 1 <= expMax then+ -- normal+ -- base^(fDigits-1) <= m / base^(k-fDigits+1) < base^fDigits+ if k < fDigits then+ -- m < base^(k+1) <= base^fDigits+ exact $ encodeFloat m n+ else+ -- k >= fDigits+ let (q,r) = quotRemByExpt m base (k - fDigits + 1)+ -- m = q * base^^(k-fDigits+1) + r+ -- base^(fDigits-1) <= q = m `quot` (base^^(k-fDigits+1)) < base^fDigits+ -- m * base^^n = q * base^^(n+k-fDigits+1) + r * base^^n+ towardzero_or_exact = encodeFloat q (n + k - fDigits + 1)+ awayfromzero = encodeFloat (q + 1) (n + k - fDigits + 1)+ parity = fromInteger q :: Int+ in doRound+ (isDivisibleByExpt m base (k - fDigits + 1) r) -- exactness (r == 0)+ (compareWithExpt base m r (k - fDigits))+ -- (compare r (expt base (k - fDigits)))+ neg+ parity+ towardzero_or_exact+ awayfromzero+ else+ if expMax <= k + n then+ -- overflow+ let inf = 1 / 0+ in inexact GT neg 1 maxFinite inf+ else -- k + n + 1 < expMin+ -- subnormal+ if expMin - fDigits <= n then+ -- k <= expMin - n <= fDigits+ exact $ encodeFloat m n+ else -- n < expMin - fDigits+ -- k <= expMin - n, fDigits < expMin - n+ let (q,r) = quotRemByExpt m base (expMin - fDigits - n)+ -- m = q * base^(expMin-fDigits-n) + r+ -- q <= m * base^^(n-expMin+fDigits) < q+1+ -- q * base^^(expMin-fDigits) <= m * base^^n < (q+1) * base^^(expMin-fDigits)+ !_ = assert (toRational q * toRational base^^(expMin-fDigits) <= toRational m * toRational base^^n) ()+ !_ = assert (toRational m * toRational base^^n < toRational (q+1) * toRational base^^(expMin-fDigits)) ()+ towardzero_or_exact = encodeFloat q (expMin - fDigits)+ awayfromzero = encodeFloat (q + 1) (expMin - fDigits)+ parity = fromInteger q :: Int+ in doRound+ (isDivisibleByExpt m base (expMin - fDigits - n) r) -- exactness (r == 0)+ (compareWithExpt base m r (expMin - fDigits - n - 1))+ -- (compare r (expt base (expMin - fDigits - n - 1)))+ neg+ parity+ towardzero_or_exact+ awayfromzero++ !base = floatRadix (undefined :: a)+ !fDigits = floatDigits (undefined :: a) -- 53 for Double+ (!expMin, !expMax) = floatRange (undefined :: a) -- (-1021, 1024) for Double+{-# INLINABLE [0] encodePositiveFloatR# #-}+{-# SPECIALIZE+ encodePositiveFloatR# :: RealFloat a => Bool -> Integer -> Int# -> RoundTiesToEven a+ , RealFloat a => Bool -> Integer -> Int# -> RoundTiesToAway a+ , RealFloat a => Bool -> Integer -> Int# -> RoundTowardPositive a+ , RealFloat a => Bool -> Integer -> Int# -> RoundTowardZero a+ , RealFloat a => Bool -> Integer -> Int# -> Product RoundTowardNegative RoundTowardPositive a+ , RoundingStrategy f => Bool -> Integer -> Int# -> f Double+ , RoundingStrategy f => Bool -> Integer -> Int# -> f Float+ , Bool -> Integer -> Int# -> RoundTiesToEven Double+ , Bool -> Integer -> Int# -> RoundTiesToAway Double+ , Bool -> Integer -> Int# -> RoundTowardPositive Double+ , Bool -> Integer -> Int# -> RoundTowardZero Double+ , Bool -> Integer -> Int# -> RoundTiesToEven Float+ , Bool -> Integer -> Int# -> RoundTiesToAway Float+ , Bool -> Integer -> Int# -> RoundTowardPositive Float+ , Bool -> Integer -> Int# -> RoundTowardZero Float+ , Bool -> Integer -> Int# -> Product RoundTowardNegative RoundTowardPositive Double+ , Bool -> Integer -> Int# -> Product RoundTowardNegative RoundTowardPositive Float+ #-}+{-# RULES+"encodePositiveFloatR#/RoundTowardNegative"+ encodePositiveFloatR# = \neg x y -> RoundTowardNegative (roundTowardPositive (encodePositiveFloatR# (not neg) x y))+ #-}++-- |+-- IEEE 754 @scaleB@ operation, with each rounding attributes.+scaleFloatTiesToEven, scaleFloatTiesToAway, scaleFloatTowardPositive, scaleFloatTowardNegative, scaleFloatTowardZero :: RealFloat a => Int -> a -> a+scaleFloatTiesToEven e = roundTiesToEven . scaleFloatR e+scaleFloatTiesToAway e = roundTiesToAway . scaleFloatR e+scaleFloatTowardPositive e = roundTowardPositive . scaleFloatR e+scaleFloatTowardNegative e = roundTowardNegative . scaleFloatR e+scaleFloatTowardZero e = roundTowardZero . scaleFloatR e+{-# INLINE scaleFloatTiesToEven #-}+{-# INLINE scaleFloatTiesToAway #-}+{-# INLINE scaleFloatTowardPositive #-}+{-# INLINE scaleFloatTowardNegative #-}+{-# INLINE scaleFloatTowardZero #-}++scaleFloatR :: (RealFloat a, RoundingStrategy f) => Int -> a -> f a+scaleFloatR (I# e#) x = scaleFloatR# e# x+{-# INLINE scaleFloatR #-}++scaleFloatR# :: (RealFloat a, RoundingStrategy f) => Int# -> a -> f a+scaleFloatR# e# x+ | x /= 0, isFinite x =+ let e = I# e#+ (m,n) = decodeFloat x+ -- x = m * base^^n, expMin <= n <= expMax+ -- base^(fDigits-1) <= abs m < base^fDigits+ -- base^(fDigits+n+e-1) <= abs x * base^^e < base^(fDigits+n+e)+ in if expMin - fDigits <= n + e && n + e <= expMax - fDigits then+ -- normal+ exact $ encodeFloat m (n + e)+ else+ if expMax - fDigits < n + e then+ -- infinity+ (signum x *) <$> inexact GT (x < 0) 1 maxFinite (1 / 0)+ else+ -- subnormal+ let !_ = assert (e + n < expMin - fDigits) ()+ m' = abs m+ (q,r) = quotRemByExpt m' base (expMin - fDigits - (e + n))+ towardzero_or_exact = encodeFloat q (expMin - fDigits)+ awayfromzero = encodeFloat (q + 1) (expMin - fDigits)+ parity = fromInteger q :: Int+ in (signum x *) <$> doRound+ (isDivisibleByExpt m' base (expMin - fDigits - (e + n)) r)+ (compareWithExpt base m' r (expMin - fDigits - (e + n) - 1))+ (x < 0)+ parity+ towardzero_or_exact+ awayfromzero+ | otherwise = exact (x + x) -- +-0, +-Infinity, NaN+ where+ base = floatRadix x+ (expMin,expMax) = floatRange x+ fDigits = floatDigits x+{-# INLINABLE [0] scaleFloatR# #-}+{-# SPECIALIZE+ scaleFloatR# :: RealFloat a => Int# -> a -> RoundTiesToEven a+ , RealFloat a => Int# -> a -> RoundTiesToAway a+ , RealFloat a => Int# -> a -> RoundTowardPositive a+ , RealFloat a => Int# -> a -> RoundTowardNegative a+ , RealFloat a => Int# -> a -> RoundTowardZero a+ , RoundingStrategy f => Int# -> Double -> f Double+ , RoundingStrategy f => Int# -> Float -> f Float+ , Int# -> Double -> RoundTiesToEven Double+ , Int# -> Double -> RoundTiesToAway Double+ , Int# -> Double -> RoundTowardPositive Double+ , Int# -> Double -> RoundTowardNegative Double+ , Int# -> Double -> RoundTowardZero Double+ , Int# -> Float -> RoundTiesToEven Float+ , Int# -> Float -> RoundTiesToAway Float+ , Int# -> Float -> RoundTowardPositive Float+ , Int# -> Float -> RoundTowardNegative Float+ , Int# -> Float -> RoundTowardZero Float+ #-}
+ src/Numeric/Floating/IEEE/Internal/Rounding/Integral.hs view
@@ -0,0 +1,316 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+module Numeric.Floating.IEEE.Internal.Rounding.Integral where+import Control.Exception (assert)+import Data.Bits+import Data.Functor.Product+import Data.Int+import Data.Proxy+import Data.Word+import GHC.Exts+import Math.NumberTheory.Logarithms (integerLog2', integerLogBase',+ wordLog2')+import MyPrelude+import Numeric.Floating.IEEE.Internal.Base+import Numeric.Floating.IEEE.Internal.IntegerInternals+import Numeric.Floating.IEEE.Internal.Rounding.Common++default ()++-- |+-- IEEE 754 @convertFromInt@ operation, with each rounding attributes.+fromIntegerTiesToEven, fromIntegerTiesToAway, fromIntegerTowardPositive, fromIntegerTowardNegative, fromIntegerTowardZero :: RealFloat a => Integer -> a+fromIntegerTiesToEven = roundTiesToEven . fromIntegerR+fromIntegerTiesToAway = roundTiesToAway . fromIntegerR+fromIntegerTowardPositive = roundTowardPositive . fromIntegerR+fromIntegerTowardNegative = roundTowardNegative . fromIntegerR+fromIntegerTowardZero = roundTowardZero . fromIntegerR+{-# INLINE fromIntegerTiesToEven #-}+{-# INLINE fromIntegerTiesToAway #-}+{-# INLINE fromIntegerTowardPositive #-}+{-# INLINE fromIntegerTowardNegative #-}+{-# INLINE fromIntegerTowardZero #-}++-- |+-- IEEE 754 @convertFromInt@ operation, with each rounding attributes.+fromIntegralTiesToEven, fromIntegralTiesToAway, fromIntegralTowardPositive, fromIntegralTowardNegative, fromIntegralTowardZero :: (Integral i, RealFloat a) => i -> a+fromIntegralTiesToEven = roundTiesToEven . fromIntegralR+fromIntegralTiesToAway = roundTiesToAway . fromIntegralR+fromIntegralTowardPositive = roundTowardPositive . fromIntegralR+fromIntegralTowardNegative = roundTowardNegative . fromIntegralR+fromIntegralTowardZero = roundTowardZero . fromIntegralR+{-# INLINE fromIntegralTiesToEven #-}+{-# INLINE fromIntegralTiesToAway #-}+{-# INLINE fromIntegralTowardPositive #-}+{-# INLINE fromIntegralTowardNegative #-}+{-# INLINE fromIntegralTowardZero #-}++fromIntegerR :: (RealFloat a, RoundingStrategy f) => Integer -> f a+fromIntegerR n = case integerToIntMaybe n of+ Just x -> fromIntegralRBits x+ Nothing | n < 0 -> negate <$> fromPositiveIntegerR True (- n)+ | otherwise -> fromPositiveIntegerR False n+{-# INLINE fromIntegerR #-}++fromIntegralR :: (Integral i, RealFloat a, RoundingStrategy f) => i -> f a+fromIntegralR x = fromIntegerR (toInteger x)+{-# INLINE [0] fromIntegralR #-}+{-# RULES+"fromIntegralR/Integer->a" fromIntegralR = fromIntegerR+"fromIntegralR/Int->a" fromIntegralR = fromIntegralRBits @Int+"fromIntegralR/Int8->a" fromIntegralR = fromIntegralRBits @Int8+"fromIntegralR/Int16->a" fromIntegralR = fromIntegralRBits @Int16+"fromIntegralR/Int32->a" fromIntegralR = fromIntegralRBits @Int32+"fromIntegralR/Int64->a" fromIntegralR = fromIntegralRBits @Int64+"fromIntegralR/Word->a" fromIntegralR = fromIntegralRBits @Word+"fromIntegralR/Word8->a" fromIntegralR = fromIntegralRBits @Word8+"fromIntegralR/Word16->a" fromIntegralR = fromIntegralRBits @Word16+"fromIntegralR/Word32->a" fromIntegralR = fromIntegralRBits @Word32+"fromIntegralR/Word64->a" fromIntegralR = fromIntegralRBits @Word64+ #-}++fromIntegralRBits :: forall i f a. (Integral i, Bits i, RealFloat a, RoundingStrategy f) => i -> f a+fromIntegralRBits x+ -- Small enough: fromIntegral should be sufficient+ | ieee+ , let resultI = fromIntegral x+ , let (min', max') = boundsForExactConversion (Proxy :: Proxy a)+ , maybe True (<= x) min'+ , maybe True (x <=) max'+ = exact resultI++ -- Signed, and not small enough: Test if the value fits in Int+ | ieee+ , base == 2+ , signed+ , Just y <- toIntegralSized x :: Maybe Int+ = if y < 0 then+ negate <$> positiveWordToBinaryFloatR True (negateIntAsWord y)+ else+ -- We can assume x /= 0+ positiveWordToBinaryFloatR False (fromIntegral y)++ -- Unsigned, and not small enough: Test if the value fits in Word+ | ieee+ , base == 2+ , not signed+ , Just y <- toIntegralSized x :: Maybe Word+ = -- We can assume x /= 0+ positiveWordToBinaryFloatR False y++ -- General case: Convert via Integer+ | otherwise = result+ where+ result | x == 0 = exact 0+ | x < 0 = negate <$> fromPositiveIntegerR True (- toInteger x)+ | otherwise = fromPositiveIntegerR False (toInteger x)+ signed = isSigned x+ ieee = isIEEE (undefined :: a)+ base = floatRadix (undefined :: a)+{-# INLINE fromIntegralRBits #-}++-- |+-- >>> boundsForExactConversion (Proxy :: Proxy Double) :: (Maybe Integer, Maybe Integer) -- (Just (-2^53),Just (2^53))+-- (Just (-9007199254740992),Just 9007199254740992)+-- >>> boundsForExactConversion (Proxy :: Proxy Double) :: (Maybe Int32, Maybe Int32) -- the conversion is always exact+-- (Nothing,Nothing)+-- >>> boundsForExactConversion (Proxy :: Proxy Float) :: (Maybe Word, Maybe Word) -- (Nothing,Just (2^24))+-- (Nothing,Just 16777216)+boundsForExactConversion :: forall a i. (Integral i, Bits i, RealFloat a) => Proxy a -> (Maybe i, Maybe i)+boundsForExactConversion _ = assert ieee (minI, maxI)+ where+ maxInteger = base ^! digits+ minInteger = - maxInteger+ minI = case minBoundAsInteger (undefined :: i) of+ Just minBound' | minInteger <= minBound' -> Nothing -- all negative integers can be expressed in the target floating-type: no check for lower-bound is needed+ _ -> Just (fromInteger minInteger)+ maxI = case maxBoundAsInteger (undefined :: i) of+ Just maxBound' | maxBound' <= maxInteger -> Nothing -- all positive integral values can be expressed in the target floating-type: no check for upper-bound is needed+ _ -> Just (fromInteger maxInteger)+ ieee = isIEEE (undefined :: a)+ base = floatRadix (undefined :: a)+ digits = floatDigits (undefined :: a)+{-# INLINE boundsForExactConversion #-}++minBoundAsInteger :: Bits i => i -> Maybe Integer+minBoundAsInteger dummyI = if isSigned dummyI then+ case bitSizeMaybe dummyI of+ Just bits -> Just (- bit (bits-1))+ Nothing -> Nothing+ else+ Just 0+{-# INLINE [1] minBoundAsInteger #-}+{-# RULES+"minBoundAsInteger/Int" minBoundAsInteger = (\_ -> Just (toInteger (minBound :: Int))) :: Int -> Maybe Integer+"minBoundAsInteger/Int8" minBoundAsInteger = (\_ -> Just (toInteger (minBound :: Int8))) :: Int8 -> Maybe Integer+"minBoundAsInteger/Int16" minBoundAsInteger = (\_ -> Just (toInteger (minBound :: Int16))) :: Int16 -> Maybe Integer+"minBoundAsInteger/Int32" minBoundAsInteger = (\_ -> Just (toInteger (minBound :: Int32))) :: Int32 -> Maybe Integer+"minBoundAsInteger/Int64" minBoundAsInteger = (\_ -> Just (toInteger (minBound :: Int64))) :: Int64 -> Maybe Integer+"minBoundAsInteger/Word" minBoundAsInteger = (\_ -> Just 0) :: Word -> Maybe Integer+"minBoundAsInteger/Word8" minBoundAsInteger = (\_ -> Just 0) :: Word8 -> Maybe Integer+"minBoundAsInteger/Word16" minBoundAsInteger = (\_ -> Just 0) :: Word16 -> Maybe Integer+"minBoundAsInteger/Word32" minBoundAsInteger = (\_ -> Just 0) :: Word32 -> Maybe Integer+"minBoundAsInteger/Word64" minBoundAsInteger = (\_ -> Just 0) :: Word64 -> Maybe Integer+ #-}++maxBoundAsInteger :: Bits i => i -> Maybe Integer+maxBoundAsInteger dummyI = case bitSizeMaybe dummyI of+ Just bits | isSigned dummyI -> Just (bit (bits-1) - 1)+ | otherwise -> Just (bit bits - 1)+ Nothing -> Nothing+{-# INLINE [1] maxBoundAsInteger #-}+{-# RULES+"maxBoundAsInteger/Int" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Int))) :: Int -> Maybe Integer+"maxBoundAsInteger/Int8" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Int8))) :: Int8 -> Maybe Integer+"maxBoundAsInteger/Int16" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Int16))) :: Int16 -> Maybe Integer+"maxBoundAsInteger/Int32" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Int32))) :: Int32 -> Maybe Integer+"maxBoundAsInteger/Int64" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Int64))) :: Int64 -> Maybe Integer+"maxBoundAsInteger/Word" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Word))) :: Word -> Maybe Integer+"maxBoundAsInteger/Word8" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Word8))) :: Word8 -> Maybe Integer+"maxBoundAsInteger/Word16" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Word16))) :: Word16 -> Maybe Integer+"maxBoundAsInteger/Word32" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Word32))) :: Word32 -> Maybe Integer+"maxBoundAsInteger/Word64" maxBoundAsInteger = (\_ -> Just (toInteger (maxBound :: Word64))) :: Word64 -> Maybe Integer+ #-}++-- Avoid cross-module specialization issue with manual worker/wrapper transformation+positiveWordToBinaryFloatR :: (RealFloat a, RoundingStrategy f) => Bool -> Word -> f a+positiveWordToBinaryFloatR neg (W# n#) = positiveWordToBinaryFloatR# neg n#+{-# INLINE positiveWordToBinaryFloatR #-}++positiveWordToBinaryFloatR# :: forall f a. (RealFloat a, RoundingStrategy f) => Bool -> Word# -> f a+positiveWordToBinaryFloatR# !neg n# = result+ where+ n = W# n#+ result = let k = wordLog2' n -- floor (log2 n)+ -- 2^k <= n < 2^(k+1) <= 2^(finiteBitSize n)+ -- k <= finiteBitSize n - 1+ in if k < fDigits then+ exact $ fromIntegral n+ else+ -- expMax <= k implies expMax <= finiteBitSize n - 1+ if expMax <= finiteBitSize n - 1 && k >= expMax then+ -- overflow+ let inf = 1 / 0+ in inexact GT neg 1 maxFinite inf+ else+ -- k >= fDigits+ let e = k - fDigits + 1 -- 1 <= e <= finiteBitSize n - fDigits+ q = n `unsafeShiftR` e -- q <= n / 2^e = 2^(log2 n - (floor (log2 n) - fDigits + 1)) < 2^fDigits+ r = n .&. ((1 `unsafeShiftL` e) - 1)+ -- (q, r) = n `quotRem` (base^e)+ -- base^(fDigits - 1) <= q < base^fDigits, 0 <= r < base^(k-fDigits+1)+ towardzero_or_exact = fromIntegral (q `unsafeShiftL` e)+ -- Although (q `unsafeShiftL` e) fits in Word, ((q + 1) `unsafeShiftL` e) may overflow.+ -- fDigits + e = k + 1 <= WORD_SIZE_IN_BITS+ -- Equality holds when wordLog2' n == WORD_SIZE_IN_BITS - 1, i.e. 2^(WORD_SIZE_IN_BITS - 1) <= n.+ -- In particular,+ -- * When q + 1 < 2^fDigits, (q + 1) * 2^e < 2^(fDigits + e) = 2^(k + 1) <= 2^WORD_SIZE_IN_BITS, so (q + 1) * 2^e does not overflow.+ -- * When k + 1 < WORD_SIZE_IN_BITS, (q + 1) * 2^e <= 2^(fDigits + e) = 2^(k+1) < 2^WORD_SIZE_IN_BITS, so (q + 1) * 2^e does not overflow.+ -- * q + 1 <= 2^fDigits and k + 1 <= WORD_SIZE_IN_BITS always hold.+ -- * Therefore, ((q + 1) `unsafeShiftL` e) overflows only if q + 1 == 2^fDigits && k + 1 == WORD_SIZE_IN_BITS+ awayfromzero = if q + 1 == (1 `unsafeShiftL` fDigits) && k == finiteBitSize n - 1 then+ -- (q + 1) `shiftL` e = 2^(fDigits + e) = 2^(k+1) = 2^(finiteBitSize n)+ encodeFloat 1 (finiteBitSize n)+ else+ fromIntegral ((q + 1) `unsafeShiftL` e)+ parity = fromIntegral q :: Int+ in doRound+ (r == 0) -- exactness+ (compare r (1 `unsafeShiftL` (e - 1)))+ neg+ parity+ towardzero_or_exact+ awayfromzero++ !fDigits = floatDigits (undefined :: a) -- 53 for Double+ (_expMin, !expMax) = floatRange (undefined :: a) -- (-1021, 1024) for Double+{-# INLINABLE [0] positiveWordToBinaryFloatR# #-}+{-# SPECIALIZE+ positiveWordToBinaryFloatR# :: RoundingStrategy f => Bool -> Word# -> f Float+ , RoundingStrategy f => Bool -> Word# -> f Double+ , RealFloat a => Bool -> Word# -> RoundTiesToEven a+ , RealFloat a => Bool -> Word# -> RoundTiesToAway a+ , RealFloat a => Bool -> Word# -> RoundTowardPositive a+ , RealFloat a => Bool -> Word# -> RoundTowardZero a+ , RealFloat a => Bool -> Word# -> Product RoundTowardNegative RoundTowardPositive a+ , Bool -> Word# -> RoundTiesToEven Float+ , Bool -> Word# -> RoundTiesToAway Float+ , Bool -> Word# -> RoundTowardPositive Float+ , Bool -> Word# -> RoundTowardZero Float+ , Bool -> Word# -> RoundTiesToEven Double+ , Bool -> Word# -> RoundTiesToAway Double+ , Bool -> Word# -> RoundTowardPositive Double+ , Bool -> Word# -> RoundTowardZero Double+ , Bool -> Word# -> Product RoundTowardNegative RoundTowardPositive Float+ , Bool -> Word# -> Product RoundTowardNegative RoundTowardPositive Double+ #-}+{-# RULES+"positiveWordToBinaryFloatR#/RoundTowardNegative"+ positiveWordToBinaryFloatR# = \neg x -> RoundTowardNegative (roundTowardPositive (positiveWordToBinaryFloatR# (not neg) x))+ #-}++-- n > 0+fromPositiveIntegerR :: forall f a. (RealFloat a, RoundingStrategy f) => Bool -> Integer -> f a+fromPositiveIntegerR !neg !n = assert (n > 0) result+ where+ result = let k = if base == 2 then+ integerLog2' n+ else+ integerLogBase' base n -- floor (logBase base n)+ -- base^k <= n < base^(k+1)+ in if k < fDigits then+ exact $ fromInteger n+ else+ if k >= expMax then+ -- overflow+ let inf = 1 / 0+ in inexact GT neg 1 maxFinite inf+ else+ -- k >= fDigits+ let e = k - fDigits + 1+ -- k >= e (assuming fDigits >= 1)+ -- Therefore, base^e <= n+ (q, r) = quotRemByExpt n base e -- n `quotRem` (base^e)+ -- base^(fDigits - 1) <= q < base^fDigits, 0 <= r < base^(k-fDigits+1)+ towardzero_or_exact = encodeFloat q e+ awayfromzero = encodeFloat (q + 1) e+ parity = fromInteger q :: Int+ in doRound+ (isDivisibleByExpt n base e r) -- exactness (r == 0)+ (compareWithExpt base n r (e - 1))+ -- (compare r (expt base (e - 1)))+ neg+ parity+ towardzero_or_exact+ awayfromzero++ !base = floatRadix (undefined :: a) -- 2 or 10+ !fDigits = floatDigits (undefined :: a) -- 53 for Double+ (_expMin, !expMax) = floatRange (undefined :: a) -- (-1021, 1024) for Double+{-# INLINABLE [0] fromPositiveIntegerR #-}+{-# SPECIALIZE+ fromPositiveIntegerR :: RealFloat a => Bool -> Integer -> RoundTiesToEven a+ , RealFloat a => Bool -> Integer -> RoundTiesToAway a+ , RealFloat a => Bool -> Integer -> RoundTowardPositive a+ , RealFloat a => Bool -> Integer -> RoundTowardZero a+ , RealFloat a => Bool -> Integer -> Product RoundTowardNegative RoundTowardPositive a+ , RoundingStrategy f => Bool -> Integer -> f Double+ , RoundingStrategy f => Bool -> Integer -> f Float+ , Bool -> Integer -> RoundTiesToEven Double+ , Bool -> Integer -> RoundTiesToAway Double+ , Bool -> Integer -> RoundTowardPositive Double+ , Bool -> Integer -> RoundTowardZero Double+ , Bool -> Integer -> RoundTiesToEven Float+ , Bool -> Integer -> RoundTiesToAway Float+ , Bool -> Integer -> RoundTowardPositive Float+ , Bool -> Integer -> RoundTowardZero Float+ , Bool -> Integer -> Product RoundTowardNegative RoundTowardPositive Double+ , Bool -> Integer -> Product RoundTowardNegative RoundTowardPositive Float+ #-}+{-# RULES+"fromPositiveIntegerR/RoundTowardNegative"+ fromPositiveIntegerR = \neg x -> RoundTowardNegative (roundTowardPositive (fromPositiveIntegerR (not neg) x))+ #-}
+ src/Numeric/Floating/IEEE/Internal/Rounding/Rational.hs view
@@ -0,0 +1,150 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Numeric.Floating.IEEE.Internal.Rounding.Rational where+import Control.Exception (assert)+import Data.Functor.Product+import Data.Ratio+import GHC.Float (expt)+import Math.NumberTheory.Logarithms (integerLog2', integerLogBase')+import MyPrelude+import Numeric.Floating.IEEE.Internal.Base+import Numeric.Floating.IEEE.Internal.Rounding.Common++default ()++-- |+-- Conversion from a rational number to floating-point value, with each rounding attributes.+fromRationalTiesToEven, fromRationalTiesToAway, fromRationalTowardPositive, fromRationalTowardNegative, fromRationalTowardZero :: RealFloat a => Rational -> a+fromRationalTiesToEven = roundTiesToEven . fromRationalR+fromRationalTiesToAway = roundTiesToAway . fromRationalR+fromRationalTowardPositive = roundTowardPositive . fromRationalR+fromRationalTowardNegative = roundTowardNegative . fromRationalR+fromRationalTowardZero = roundTowardZero . fromRationalR+{-# INLINE fromRationalTiesToEven #-}+{-# INLINE fromRationalTiesToAway #-}+{-# INLINE fromRationalTowardPositive #-}+{-# INLINE fromRationalTowardNegative #-}+{-# INLINE fromRationalTowardZero #-}++fromRationalR :: (RealFloat a, RoundingStrategy f) => Rational -> f a+fromRationalR x = fromRatioR (numerator x) (denominator x)+{-# INLINE fromRationalR #-}++fromRatioR :: (RealFloat a, RoundingStrategy f)+ => Integer -- ^ numerator+ -> Integer -- ^ denominator+ -> f a+fromRatioR 0 !_ = exact 0+fromRatioR n 0 | n > 0 = exact (1 / 0) -- positive infinity+ | otherwise = exact (- 1 / 0) -- negative infinity+fromRatioR n d | d < 0 = error "fromRatio: negative denominator"+ | n < 0 = negate <$> fromPositiveRatioR True (- n) d+ | otherwise = fromPositiveRatioR False n d+{-# INLINE fromRatioR #-}++fromPositiveRatioR :: forall f a. (RealFloat a, RoundingStrategy f)+ => Bool -- ^ True if the result will be negated+ -> Integer -- ^ numerator (> 0)+ -> Integer -- ^ denominator (> 0)+ -> f a+fromPositiveRatioR !neg !n !d = assert (n > 0 && d > 0) result+ where+ result = let e0 :: Int+ e0 = if base == 2 then+ integerLog2' n - integerLog2' d - fDigits+ else+ integerLogBase' base n - integerLogBase' base d - fDigits+ q0, r0, d0 :: Integer+ (!d0, (!q0, !r0)) =+ if e0 >= 0 then+ -- n = q0 * (d * base^e0) + r0, 0 <= r0 < d * base^e0+ let d_ = multiplyByExpt d base e0+ in (d_, n `quotRem` d_)+ else+ -- n * base^(-e0) = q0 * d + r0, 0 <= r0 < d+ (d, (multiplyByExpt n base (-e0)) `quotRem` d)+ -- Invariant: n / d * base^^(-e0) = q0 + r0 / d0+ !_ = assert (n % d * fromInteger base^^(-e0) == fromInteger q0 + r0 % d0) ()+ !_ = assert (base^(fDigits-1) <= q0 && q0 < base^(fDigits+1)) ()++ q, r, d' :: Integer+ e :: Int+ (!q, !r, !d', !e) =+ if q0 < expt base fDigits then+ -- base^(fDigits-1) <= q0 < base^fDigits+ (q0, r0, d0, e0)+ else+ -- base^fDigits <= q0 < base^(fDigits+1)+ let (q', r') = q0 `quotRem` base+ in (q', r' * d0 + r0, base * d0, e0 + 1)+ -- Invariant: n / d * 2^^(-e) = q + r / d', base^(fDigits-1) <= q < base^fDigits, 0 <= r < d'+ !_ = assert (n % d * fromInteger base^^(-e) == fromInteger q + r % d') ()+ -- base^(e+fDigits-1) <= q * base^^e <= n/d < (q+1) * base^^e <= base^(e+fDigits)+ -- In particular, base^(fDigits-1) <= q < base^fDigits+ in if expMin <= e + fDigits && e + fDigits <= expMax then+ -- normal: base^^(expMin-1) <= n/d < base^expMax+ let towardzero_or_exact = encodeFloat q e+ awayfromzero = encodeFloat (q + 1) e -- may be infinity+ parity = fromInteger q :: Int+ in doRound+ (r == 0)+ (compare (base * r) d')+ neg+ parity+ towardzero_or_exact+ awayfromzero+ else+ if expMax < e + fDigits then+ -- overflow+ let inf = 1 / 0+ in inexact GT neg 1 maxFinite inf+ else+ -- subnormal: 0 < n/d < base^^(expMin-1)+ -- e + fDigits < expMin+ let (q', r') = quotRemByExpt q base (expMin - fDigits - e)+ !_ = assert (q == q' * base^(expMin-fDigits-e) + r' && 0 <= r' && r' < base^(expMin-fDigits-e)) ()+ -- q = q' * base^(expMin-fDigits-e) + r', 0 <= r' < base^(expMin-fDigits-e)+ -- n / d * base^^(-e) = q' * base^(expMin-fDigits-e) + r' + r / d'+ -- n / d = q' * base^^(expMin - fDigits) + (r' + r / d') * base^^e+ !_ = assert (n % d == fromInteger q' * fromInteger base^^(expMin - fDigits) + (fromInteger r' + r % d') * fromInteger base^^e) ()+ -- rounding direction: (r' + r / d') * base^^e vs. base^^(expMin-fDigits-1)+ towardzero = encodeFloat q' (expMin - fDigits)+ awayfromzero = encodeFloat (q' + 1) (expMin - fDigits)+ parity = fromInteger q' :: Int+ in doRound+ (r == 0 && r' == 0)+ (compareWithExpt base q r' (expMin - fDigits - e - 1) <> if r == 0 then EQ else GT)+ -- (compare r' (expt base (expMin - fDigits - e - 1)) <> if r == 0 then EQ else GT)+ neg+ parity+ towardzero+ awayfromzero++ !base = floatRadix (undefined :: a)+ !fDigits = floatDigits (undefined :: a) -- 53 for Double+ (!expMin, !expMax) = floatRange (undefined :: a) -- (-1021, 1024) for Double+{-# INLINABLE [0] fromPositiveRatioR #-}+{-# SPECIALIZE+ fromPositiveRatioR :: RealFloat a => Bool -> Integer -> Integer -> RoundTiesToEven a+ , RealFloat a => Bool -> Integer -> Integer -> RoundTiesToAway a+ , RealFloat a => Bool -> Integer -> Integer -> RoundTowardPositive a+ , RealFloat a => Bool -> Integer -> Integer -> RoundTowardZero a+ , RealFloat a => Bool -> Integer -> Integer -> Product RoundTowardNegative RoundTowardPositive a+ , RoundingStrategy f => Bool -> Integer -> Integer -> f Double+ , RoundingStrategy f => Bool -> Integer -> Integer -> f Float+ , Bool -> Integer -> Integer -> RoundTiesToEven Double+ , Bool -> Integer -> Integer -> RoundTiesToAway Double+ , Bool -> Integer -> Integer -> RoundTowardPositive Double+ , Bool -> Integer -> Integer -> RoundTowardZero Double+ , Bool -> Integer -> Integer -> RoundTiesToEven Float+ , Bool -> Integer -> Integer -> RoundTiesToAway Float+ , Bool -> Integer -> Integer -> RoundTowardPositive Float+ , Bool -> Integer -> Integer -> RoundTowardZero Float+ , Bool -> Integer -> Integer -> Product RoundTowardNegative RoundTowardPositive Double+ , Bool -> Integer -> Integer -> Product RoundTowardNegative RoundTowardPositive Float+ #-}+{-# RULES+"fromPositiveRatioR/RoundTowardNegative"+ fromPositiveRatioR = \neg x y -> RoundTowardNegative (roundTowardPositive (fromPositiveRatioR (not neg) x y))+ #-}
+ src/Numeric/Floating/IEEE/NaN.hs view
@@ -0,0 +1,15 @@+{-|+Module : Numeric.Floating.IEEE.NaN+Description : Accessing the sign and payload of NaNs++This module provides the typeclass for NaN manipulation: 'RealFloatNaN'.++In addition to 'Float' and 'Double', a couple of floating-point types provided by third-party libraries can be supported via package flags: @Half@ via @half@ and @Float128@ via @float128@.+-}+module Numeric.Floating.IEEE.NaN+ ( RealFloatNaN(..)+ , Class(..)+ , TotallyOrdered(..)+ ) where+import Numeric.Floating.IEEE.Internal ()+import Numeric.Floating.IEEE.Internal.NaN
+ test/AugmentedArithSpec.hs view
@@ -0,0 +1,111 @@+{-# LANGUAGE HexFloatLiterals #-}+{-# LANGUAGE ScopedTypeVariables #-}+module AugmentedArithSpec where+import Control.Monad+import Numeric+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import RoundingSpec (RoundTiesTowardZero (..))+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import Util++augmentedAddition_viaRational :: (RealFloat a, Show a) => a -> a -> (a, a)+augmentedAddition_viaRational x y+ | isFinite x && isFinite y && (x /= 0 || y /= 0) =+ let z :: Rational+ z = toRational x + toRational y+ z' = roundTiesTowardZero (fromRationalR z) `asTypeOf` x+ in if isInfinite z' then+ (z', z')+ else+ let w :: Rational+ w = z - toRational z'+ w' = roundTiesTowardZero (fromRationalR w) `asTypeOf` x+ in if w == 0 then+ (z', 0 * z')+ else+ (z', w')+ | otherwise = let z = x + y+ in (z, z)++augmentedMultiplication_viaRational :: (RealFloat a, Show a) => a -> a -> (a, a)+augmentedMultiplication_viaRational x y+ | isFinite x && isFinite y && x * y /= 0 =+ let z :: Rational+ z = toRational x * toRational y+ z' = roundTiesTowardZero (fromRationalR z) `asTypeOf` x+ in if isInfinite z' then+ (z', z')+ else+ let w :: Rational+ w = z - toRational z'+ w' = roundTiesTowardZero (fromRationalR w) `asTypeOf` x+ in if w == 0 then+ (z', 0 * z')+ else+ (z', w')+ | otherwise = let z = x * y+ in (z, z)++testAugmented :: (RealFloat a, Show a) => (a -> a -> (a, a)) -> [(a, a, a, a)] -> Property+testAugmented f cases = conjoin+ [ let label = showHFloat a . showChar ' ' . showHFloat b $ ""+ in counterexample label $ f a b `sameFloatPairP` (r1,r2)+ | (a,b,r1,r2) <- cases+ ]++{-# NOINLINE spec #-}+spec :: Spec+spec = modifyMaxSuccess (* 100) $ do+ describe "Double" $ do+ do -- augmentedAddition+ prop "augmentedAddition/equality" $ forAllFloats2 $ \(x :: Double) y ->+ isFinite x && isFinite y ==>+ let (s,t) = augmentedAddition x y+ in isFinite s ==> isFinite t .&&. toRational s + toRational t === toRational x + toRational y+ let cases :: [(Double, Double, Double, Double)]+ cases = [ (-0, -0, -0, -0)+ ]+ prop "augmentedAddition" $ testAugmented augmentedAddition cases+ prop "augmentedAddition_viaRational" $ testAugmented augmentedAddition_viaRational cases+ prop "augmentedAddition" $ forAllFloats2 $ \(x :: Double) y ->+ augmentedAddition x y `sameFloatPairP` augmentedAddition_viaRational x y++ do -- augmentedMultiplication+ let cases :: [(Double, Double, Double, Double)]+ cases = [ (-0x1.3deed726aad4p-1023, 0x1.e179bde0a1dd2p-1, -0x1.2afa79f9d38c6p-1023, 0x0p+0)+ , (-0x1.8eb0e02044f68p-1022, -0x1.c93b83a5751c8p-2, 0x1.640b37f1b9d02p-1023,-0x0p+0)+ , (0x1.b877a1cd61478p-1023, -0x1.7a77bb9df06dap-1, -0x1.459753aa4d2bep-1023, -0x0p+0)+ , (-0x1.d25f2402fe726p-1, -0x1.0b42f4e9eb842p-1, 0x1.e6e335433c1f9p-2, -0x1.bb70c80f1834p-58)+ ]+ prop "augmentedMultiplication" $ testAugmented augmentedMultiplication cases+ prop "augmentedMultiplication_viaRational" $ testAugmented augmentedMultiplication_viaRational cases+ prop "augmentedMultiplication" $ forAllFloats2 $ \(x :: Double) y ->+ augmentedMultiplication x y `sameFloatPairP` augmentedMultiplication_viaRational x y++ describe "Float" $ do+ do -- augmentedAddition+ prop "augmentedAddition/equality" $ forAllFloats2 $ \(x :: Float) y ->+ isFinite x && isFinite y ==>+ let (s,t) = augmentedAddition x y+ in isFinite s ==> isFinite t .&&. toRational s + toRational t === toRational x + toRational y+ let cases :: [(Float, Float, Float, Float)]+ cases = [(-0, -0, -0, -0)]+ prop "augmentedAddition" $ testAugmented augmentedAddition cases+ prop "augmentedAddition_viaRational" $ testAugmented augmentedAddition_viaRational cases+ prop "augmentedAddition" $ forAllFloats2 $ \(x :: Float) y ->+ augmentedAddition x y `sameFloatPairP` augmentedAddition_viaRational x y++ do -- augmentedMultiplication+ let cases :: [(Float, Float, Float, Float)]+ cases = [ (0x1.b8508p-130, -0x1.93994p-4, -0x1.5b17p-133, -0x0p+0)+ , (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)+ ]+ prop "augmentedMultiplication" $ testAugmented augmentedMultiplication cases+ prop "augmentedMultiplication_viaRational" $ testAugmented augmentedMultiplication_viaRational cases+ prop "augmentedMultiplication" $ forAllFloats2 $ \(x :: Float) y ->+ augmentedMultiplication x y `sameFloatPairP` augmentedMultiplication_viaRational x y
+ test/ClassificationSpec.hs view
@@ -0,0 +1,63 @@+module ClassificationSpec where+import Data.Function (on)+import Data.Functor.Identity+import Data.Proxy+import Numeric.Floating.IEEE+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import Util++default ()++prop_classify :: (RealFloat a, Show a) => Proxy a -> a -> Property+prop_classify _ x = conjoin+ [ counterexample "NegativeInfinity" $ (c == NegativeInfinity) === (x < 0 && isInfinite x)+ , counterexample "NegativeNormal" $ (c == NegativeNormal) === (x < 0 && isNormal x)+ , counterexample "NegativeSubnormal" $ (c == NegativeSubnormal) === (x < 0 && isDenormalized x)+ , counterexample "NegativeZero" $ (c == NegativeZero) === (isNegativeZero x)+ , counterexample "PositiveZero" $ (c == PositiveZero) === (x == 0 && not (isNegativeZero x))+ , counterexample "PositiveSubnormal" $ (c == PositiveSubnormal) === (x > 0 && isDenormalized x)+ , counterexample "PositiveNormal" $ (c == PositiveNormal) === (x > 0 && isNormal x)+ , counterexample "PositiveInfinity" $ (c == PositiveInfinity) === (x > 0 && isInfinite x)+ , counterexample "isNaN" $ isNaN x === (c == SignalingNaN || c == QuietNaN)+ , counterexample "isInfinite" $ isInfinite x === (c == NegativeInfinity || c == PositiveInfinity)+ , counterexample "isNormal" $ isNormal x === (c == NegativeNormal || c == PositiveNormal)+ , counterexample "isDenormalized" $ isDenormalized x === (c == NegativeSubnormal || c == PositiveSubnormal)+ , counterexample "isZero" $ isZero x === (c == NegativeZero || c == PositiveZero)+ , counterexample "isFinite" $ isFinite x === (c `elem` [NegativeNormal, NegativeSubnormal, NegativeZero, PositiveZero, PositiveSubnormal, PositiveNormal])+ , counterexample "isSignMinus" $ isSignMinus x === (c `elem` [NegativeInfinity, NegativeNormal, NegativeSubnormal, NegativeZero]) -- isSignMinus doesn't handle negative NaNs+ ]+ where c = classify x+{-# SPECIALIZE prop_classify :: Proxy Float -> Float -> Property, Proxy Double -> Double -> Property #-}++prop_totalOrder :: RealFloat a => Proxy a -> a -> a -> Property+prop_totalOrder proxy x y = let cmp_x_y = compareByTotalOrder x y+ cmp_y_x = compareByTotalOrder y x+ in cmp_x_y === compare EQ cmp_y_x+ .&&. (if x < y then cmp_x_y === LT else property True)+ .&&. (if y < x then cmp_x_y === GT else property True)+{-# SPECIALIZE prop_totalOrder :: Proxy Float -> Float -> Float -> Property, Proxy Double -> Double -> Double -> Property #-}++spec :: Spec+spec = do+ describe "Double" $ do+ let proxy :: Proxy Double+ proxy = Proxy+ prop "classify" $ forAllFloats $ prop_classify proxy+ prop "totalOrder" $ forAllFloats2 $ prop_totalOrder proxy+ describe "Double (generic)" $ do+ let proxy :: Proxy (Identity Double)+ proxy = Proxy+ prop "classify" $ forAllFloats $ prop_classify proxy . Identity+ prop "totalOrder" $ forAllFloats2 (prop_totalOrder proxy `on` Identity)+ describe "Float" $ do+ let proxy :: Proxy Float+ proxy = Proxy+ prop "classify" $ forAllFloats $ prop_classify proxy+ prop "totalOrder" $ forAllFloats2 $ prop_totalOrder proxy+ describe "Float (generic)" $ do+ let proxy :: Proxy (Identity Float)+ proxy = Proxy+ prop "classify" $ forAllFloats $ prop_classify proxy . Identity+ prop "totalOrder" $ forAllFloats2 (prop_totalOrder proxy `on` Identity)
+ test/FMASpec.hs view
@@ -0,0 +1,107 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE HexFloatLiterals #-}+module FMASpec where+import Control.Monad+import Data.Bits+import Data.Coerce+import Data.Functor.Identity+import Numeric+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import System.Random+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck+import Util (forAllFloats3, sameFloatP)++#if defined(USE_FFI)++foreign import ccall unsafe "fma"+ c_fma_double :: Double -> Double -> Double -> Double+foreign import ccall unsafe "fmaf"+ c_fma_float :: Float -> Float -> Float -> Float++#endif++fusedMultiplyAdd_generic :: RealFloat a => a -> a -> a -> a+fusedMultiplyAdd_generic x y z = runIdentity (fusedMultiplyAdd (Identity x) (Identity y) (Identity z))++fusedMultiplyAdd_viaInteger :: RealFloat a => a -> a -> a -> a+fusedMultiplyAdd_viaInteger x y z+ | isFinite x && isFinite y && isFinite z =+ let (mx,ex) = decodeFloat x -- x == mx * b^ex, mx==0 || b^(d-1) <= abs mx < b^d+ (my,ey) = decodeFloat y -- y == my * b^ey, my==0 || b^(d-1) <= abs my < b^d+ (mz,ez) = decodeFloat z -- z == mz * b^ez, mz==0 || b^(d-1) <= abs mz < b^d+ exy = ex + ey+ ee = min ez exy+ !2 = floatRadix x+ in case mx * my `shiftL` (exy - ee) + mz `shiftL` (ez - ee) of+ 0 -> x * y + z+ m -> roundTiesToEven (encodeFloatR m ee)+ | isFinite x && isFinite y = z + z -- x * y is finite, but z is Infinity or NaN+ | otherwise = x * y + z -- either x or y is Infinity or NaN++fusedMultiplyAdd_viaRational :: RealFloat a => a -> a -> a -> a+fusedMultiplyAdd_viaRational x y z+ | isFinite x && isFinite y && isFinite z =+ case toRational x * toRational y + toRational z of+ 0 -> x * y + z+ r -> fromRational r+ | isFinite x && isFinite y = z + z -- x * is finite, but z is Infinity or NaN+ | otherwise = x * y + z -- either x or y is Infinity or NaN++casesForDouble :: [(Double, Double, Double, Double)]+casesForDouble =+ [ (0x1.af7da9fc47b3ep-1, 0x1p-1074, -0x1p-1074, -0)+ , (0x1p512, 0x1p512, -0x1p1023, 0x1p1023)+ , (0x1.0000000000008p500, 0x1.1p500, 0x1p-1074, 0x1.1000000000009p1000)+ , (0x1.0000000000001p500, 0x1.8p500, -0x1p-1074, 0x1.8000000000001p1000)+ , (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)+ ]++casesForFloat :: [(Float, Float, Float, Float)]+casesForFloat =+ [ (16777215, 268435520, 63.5, 0x1.000002p52)+ , (0x1.84ae30p125, 0x1.6p-141, 0x1p-149, 0x1.0b37c2p-15)+ , (0x1.000010p50, 0x1.1p50, 0x1p-149, 0x1.100012p100)+ , (0x1.000002p50, 0x1.8p50, -0x1p-149, 0x1.800002p100)+ , (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)+ ]++testSpecialValues :: (RealFloat a, Show a) => String -> (a -> a -> a -> a) -> [(a, a, a, a)] -> Spec+testSpecialValues name f cases = forM_ cases $ \(a,b,c,result) -> do+ let label = showString name . showChar ' ' . showHFloat a . showChar ' ' . showHFloat b . showChar ' ' . showHFloat c . showString " should be " . showHFloat result $ ""+ it label $ f a b c `sameFloatP` result++checkFMA :: (RealFloat a, Show a, Arbitrary a, Random a) => String -> (a -> a -> a -> a) -> [(a, a, a, a)] -> Spec+checkFMA name f cases = do+ prop name $ forAllFloats3 $ \a b c -> do+ f a b c `sameFloatP` fusedMultiplyAdd_viaRational a b c+ testSpecialValues name f cases++spec :: Spec+spec = modifyMaxSuccess (* 100) $ do+ describe "Double" $ do+ checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd 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 (generic)" fusedMultiplyAdd_generic casesForFloat+ checkFMA "fusedMultiplyAdd (via Rational)" fusedMultiplyAdd_viaRational casesForFloat+ checkFMA "fusedMultiplyAdd (via Integer)" fusedMultiplyAdd_viaInteger casesForFloat+ checkFMA "fusedMultiplyAdd (via Double)" fusedMultiplyAddFloat_viaDouble casesForFloat+#if defined(USE_FFI)+ describe "Extra" $ do+ describe "Double" $ do+ checkFMA "C fma" c_fma_double casesForDouble+ describe "Float" $ do+ checkFMA "C fmaf" c_fma_float casesForFloat+#endif+{-# NOINLINE spec #-}
+ test/Float128Spec.hs view
@@ -0,0 +1,108 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE HexFloatLiterals #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wno-orphans #-}+module Float128Spec where+import AugmentedArithSpec (augmentedAddition_viaRational,+ augmentedMultiplication_viaRational)+import qualified AugmentedArithSpec+import qualified ClassificationSpec+import Control.Monad+import Data.Function (on)+import Data.Functor.Identity+import Data.Int+import Data.Proxy+import Data.Ratio+import FMASpec (fusedMultiplyAdd_generic,+ fusedMultiplyAdd_viaRational)+import qualified FMASpec+import qualified NaNSpec+import qualified NextFloatSpec+import Numeric.Float128+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import Numeric.Floating.IEEE.NaN (setPayloadSignaling)+import qualified RoundingSpec+import qualified RoundToIntegralSpec+import System.Random+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import TwoSumSpec (twoProduct_generic)+import qualified TwoSumSpec+import Util++-- orphan instances+instance Arbitrary Float128 where+ arbitrary = arbitrarySizedFractional+ shrink = shrinkDecimal++instance Random Float128 where+ -- Float128:+ -- emin = -14, emax = 15+ -- precision = 11 bits+ -- maxFinite = 0xffe0 (65504)+ randomR (lo,hi) g = let (x,g') = random g+ in (lo + x * (hi - lo), g') -- TODO: avoid overflow+ random g = let x :: Int64+ (x,g') = random g+ in (fromRational (toInteger x % 2^(16 :: Int)), g') -- TODO++spec :: Spec+spec = mapSpecItem_ (allowFailure "Float128's fromRational and round may be incorrect") $ do+ let proxy :: Proxy Float128+ proxy = Proxy+ prop "classify" $ forAllFloats $ ClassificationSpec.prop_classify proxy+ prop "classify (generic)" $ forAllFloats $ ClassificationSpec.prop_classify (Proxy :: Proxy (Identity Float128)) . Identity+ 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+ let casesForFloat128 :: [(Float128, Float128, Float128, Float128)]+ casesForFloat128 = [ (-0, 0, -0, -0)+ , (-0, -0, -0, 0)+ -- TODO: Add more+ ]+ FMASpec.checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForFloat128+ FMASpec.checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForFloat128+ 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 ->+ 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 :: Float128) y ->+ augmentedAddition x y `sameFloatPairP` augmentedAddition_viaRational x y+ prop "augmentedMultiplication" $ forAllFloats2 $ \(x :: Float128) y ->+ augmentedMultiplication x y `sameFloatPairP` augmentedMultiplication_viaRational x y++ 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+ 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)+ 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++ prop "roundToIntegral" $ RoundToIntegralSpec.prop_roundToIntegral proxy+ RoundToIntegralSpec.checkCases proxy++ prop "copySign" $ forAllFloats2 $ NaNSpec.prop_copySign proxy+ prop "isSignMinus" $ forAllFloats $ NaNSpec.prop_isSignMinus proxy+ prop "isSignaling" $ NaNSpec.prop_isSignaling proxy+ prop "setPayload/getPayload" $ NaNSpec.prop_setPayload_getPayload proxy+ prop "setPayload/0" $ NaNSpec.prop_setPayload proxy 0+ prop "setPayload/0x1p9" $ NaNSpec.prop_setPayload proxy 0x1p9+ prop "setPayload/Int" $ NaNSpec.prop_setPayload proxy . (fromIntegral :: Int -> Float128)+ prop "setPayloadSignaling/0" $ NaNSpec.prop_setPayloadSignaling proxy 0+ prop "setPayloadSignaling/0x1p9" $ NaNSpec.prop_setPayloadSignaling proxy 0x1p9+ prop "setPayloadSignaling/Int" $ NaNSpec.prop_setPayloadSignaling proxy . (fromIntegral :: Int -> Float128)+ prop "classify" $ forAllFloats $ NaNSpec.prop_classify proxy+ prop "classify (signaling NaN)" $ NaNSpec.prop_classify proxy (setPayloadSignaling 123)+ prop "signaling NaN propagation" $ NaNSpec.prop_signalingNaN proxy+ prop "totalOrder" $ forAllFloats2 $ NaNSpec.prop_totalOrder proxy
+ test/HalfSpec.hs view
@@ -0,0 +1,123 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE HexFloatLiterals #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wno-orphans #-}+module HalfSpec where+import AugmentedArithSpec (augmentedAddition_viaRational,+ augmentedMultiplication_viaRational)+import qualified AugmentedArithSpec+import qualified ClassificationSpec+import Control.Monad+import Data.Function (on)+import Data.Functor.Identity+import Data.Int+import Data.Proxy+import Data.Ratio+import FMASpec (fusedMultiplyAdd_generic,+ fusedMultiplyAdd_viaRational)+import qualified FMASpec+import qualified NaNSpec+import qualified NextFloatSpec+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import Numeric.Floating.IEEE.NaN (setPayloadSignaling)+import Numeric.Half+import qualified RoundingSpec+import qualified RoundToIntegralSpec+import System.Random+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import TwoSumSpec (twoProduct_generic)+import qualified TwoSumSpec+import Util++-- orphan instances+instance Arbitrary Half where+ arbitrary = arbitrarySizedFractional+ shrink = shrinkDecimal++instance Random Half where+ -- Half:+ -- emin = -14, emax = 15+ -- precision = 11 bits+ -- maxFinite = 0xffe0 (65504)+ randomR (lo,hi) g = let (x,g') = random g+ in (lo + x * (hi - lo), g') -- TODO: avoid overflow+ random g = let x :: Int32+ (x,g') = random g+ in (fromRational (toInteger x % 2^(16 :: Int)), g')++isInfiniteWorkaround :: (Half -> Property) -> (Half -> Property)+isInfiniteIsKnownToBeBuggy :: Bool+#if MIN_VERSION_half(0,3,1)+-- I hope https://github.com/ekmett/half/issues/23 is fixed before the next releaso+isInfiniteWorkaround = id+isInfiniteIsKnownToBeBuggy = False+#else+isInfiniteWorkaround f x = not (isNaN x) ==> f x+isInfiniteIsKnownToBeBuggy = True+#endif++spec :: Spec+spec = mapSpecItem_ (allowFailure "Half's fromRational may be incorrect") $ do+ let proxy :: Proxy Half+ proxy = Proxy+ prop "classify" $ forAllFloats $ isInfiniteWorkaround $ ClassificationSpec.prop_classify proxy+ prop "classify (generic)" $ forAllFloats $ isInfiniteWorkaround $ ClassificationSpec.prop_classify (Proxy :: Proxy (Identity Half)) . Identity+ prop "totalOrder" $ forAllFloats2 $ ClassificationSpec.prop_totalOrder proxy+ prop "totalOrder (generic)" $ forAllFloats2 (ClassificationSpec.prop_totalOrder (Proxy :: Proxy (Identity Half)) `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+ let casesForHalf :: [(Half, Half, Half, Half)]+ casesForHalf = [ (-0, 0, -0, -0)+ , (-0, -0, -0, 0)+ -- TODO: Add more+ ]+ FMASpec.checkFMA "fusedMultiplyAdd (default)" fusedMultiplyAdd casesForHalf+ FMASpec.checkFMA "fusedMultiplyAdd (generic)" fusedMultiplyAdd_generic casesForHalf+ 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 "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 ->+ augmentedMultiplication x y `sameFloatPairP` augmentedMultiplication_viaRational x y++ 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+ 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)+ 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++ prop "roundToIntegral" $ RoundToIntegralSpec.prop_roundToIntegral proxy+ RoundToIntegralSpec.checkCases proxy++ prop "copySign" $ forAllFloats2 $ NaNSpec.prop_copySign proxy+ prop "isSignMinus" $ forAllFloats $ NaNSpec.prop_isSignMinus proxy+ prop "isSignaling" $ NaNSpec.prop_isSignaling proxy+ prop "setPayload/getPayload" $ NaNSpec.prop_setPayload_getPayload proxy+ prop "setPayload/0" $ NaNSpec.prop_setPayload proxy 0+ prop "setPayload/0x1p9" $ NaNSpec.prop_setPayload proxy 0x1p9+ prop "setPayload/Int" $ NaNSpec.prop_setPayload proxy . (fromIntegral :: Int -> Half)+ prop "setPayloadSignaling/0" $ NaNSpec.prop_setPayloadSignaling proxy 0+ prop "setPayloadSignaling/0x1p9" $ NaNSpec.prop_setPayloadSignaling proxy 0x1p9+ prop "setPayloadSignaling/Int" $ NaNSpec.prop_setPayloadSignaling proxy . (fromIntegral :: Int -> Half)+ prop "classify" $ forAllFloats $ isInfiniteWorkaround $ NaNSpec.prop_classify proxy+ when (not isInfiniteIsKnownToBeBuggy) $ do+ prop "classify (signaling NaN)" $ NaNSpec.prop_classify proxy (setPayloadSignaling 123)+ prop "signaling NaN propagation" $ NaNSpec.prop_signalingNaN proxy+ prop "totalOrder" $ forAllFloats2 $ NaNSpec.prop_totalOrder proxy++ when isInfiniteIsKnownToBeBuggy $ do+ runIO $ putStrLn "Half's isInfinite is known to be buggy on this version. Some tests were skipped."
+ test/IntegerInternalsSpec.hs view
@@ -0,0 +1,53 @@+module IntegerInternalsSpec (spec) where+import Data.Bits+import Data.Int+import Data.Maybe+import Math.NumberTheory.Logarithms+import Numeric.Floating.IEEE.Internal+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import Util++default ()++noinline :: a -> a+noinline = id+{-# NOINLINE noinline #-}++{-# NOINLINE spec #-}+spec :: Spec+spec = do+ describe "integerToIntMaybe" $ do+ it "0" $ integerToIntMaybe 0 `shouldBe` Just 0+ it "123" $ integerToIntMaybe 123 `shouldBe` Just 123+ it "minBound :: Int" $ integerToIntMaybe (toInteger (minBound :: Int)) `shouldBe` Just minBound+ it "maxBound :: Int" $ integerToIntMaybe (toInteger (maxBound :: Int)) `shouldBe` Just maxBound+ it "(minBound :: Int) - 1" $ integerToIntMaybe (toInteger (minBound :: Int) - 1) `shouldBe` Nothing+ it "(maxBound :: Int) + 1" $ integerToIntMaybe (toInteger (maxBound :: Int) + 1) `shouldBe` Nothing+ prop "small integer" $ \x -> integerToIntMaybe (toInteger x) `shouldBe` Just x++ describe "integerToIntMaybe/noinline" $ do+ it "0" $ noinline integerToIntMaybe 0 `shouldBe` Just 0+ it "123" $ noinline integerToIntMaybe 123 `shouldBe` Just 123+ it "minBound :: Int" $ noinline integerToIntMaybe (toInteger (minBound :: Int)) `shouldBe` Just minBound+ it "maxBound :: Int" $ noinline integerToIntMaybe (toInteger (maxBound :: Int)) `shouldBe` Just maxBound+ it "(minBound :: Int) - 1" $ noinline integerToIntMaybe (toInteger (minBound :: Int) - 1) `shouldBe` Nothing+ it "(maxBound :: Int) + 1" $ noinline integerToIntMaybe (toInteger (maxBound :: Int) + 1) `shouldBe` Nothing+ prop "small integer" $ \x -> noinline integerToIntMaybe (toInteger x) `shouldBe` Just x++ prop "unsafeShiftLInteger" $ \x (NonNegative y) -> unsafeShiftLInteger x y `shouldBe` shiftL x y+ prop "unsafeShiftRInteger" $ \x (NonNegative y) -> unsafeShiftRInteger x y `shouldBe` shiftR x y++ describe "roundingMode" $ do+ prop "prop" $ \(Positive n) -> forAll (choose (0, integerLog2 n)) $ \e -> integerLog2 n >= e ==> roundingMode n e `shouldBe` compare (n `rem` 2^(e+1 :: Int)) (2^e)++ describe "countTrailingZerosInteger" $ do+ prop "test with Int64" $ \(NonZero x) -> countTrailingZerosInteger (fromIntegral x) == countTrailingZeros (x :: Int64)++ describe "integerIsPowerOf2" $ do+ prop "power of 2" $ \(NonNegative x) -> integerIsPowerOf2 (2^(x :: Int)) `shouldBe` Just x+ prop "(power of 2) + 1" $ \(Positive x) -> integerIsPowerOf2 (2^(x :: Int) + 1) `shouldBe` Nothing+ prop "(power of 2) - 1" $ \(Positive x) -> integerIsPowerOf2 (2^(x+1 :: Int) - 1) `shouldBe` Nothing++ prop "integerLog2IsPowerOf2" $ \(Positive x) -> integerLog2IsPowerOf2 x `shouldBe` (integerLog2 x, isJust (integerIsPowerOf2 x))
+ test/MinMaxSpec.hs view
@@ -0,0 +1,134 @@+module MinMaxSpec where+import Data.Coerce+import Data.Functor.Identity+import Data.Proxy+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import Numeric.Floating.IEEE.NaN (RealFloatNaN(..))+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck+import Util++default ()++isQuietNaN :: RealFloatNaN a => a -> Bool+isQuietNaN x = isNaN x && not (isSignaling x)++prop_minimum :: RealFloatNaN a => Proxy a -> (a -> a -> a) -> Property+prop_minimum _ m =+ let sNaN = setPayloadSignaling 1+ qNaN = setPayload 1+ in conjoin+ [ counterexample "(1,3)" $ m 1 3 `sameFloatP` 1+ , counterexample "(1,-1)" $ m 1 (-1) `sameFloatP` (-1)+ , counterexample "(0,0)" $ m 0 0 `sameFloatP` 0+ , counterexample "(0,-0)" $ m 0 (-0) `sameFloatP` (-0)+ , counterexample "(-0,0)" $ m (-0) 0 `sameFloatP` (-0)+ , counterexample "(-0,-0)" $ m (-0) (-0) `sameFloatP` (-0)+ , counterexample "(sNaN,sNaN)" $ isQuietNaN (m sNaN sNaN)+ , counterexample "(sNaN,qNaN)" $ isQuietNaN (m sNaN qNaN)+ , counterexample "(qNaN,sNaN)" $ isQuietNaN (m qNaN sNaN)+ , counterexample "(qNaN,qNaN)" $ isQuietNaN (m qNaN qNaN)+ , counterexample "(sNaN,1.0)" $ isQuietNaN (m sNaN 1.0)+ , counterexample "(1.0,sNaN)" $ isQuietNaN (m 1.0 sNaN)+ , counterexample "(qNaN,1.0)" $ isQuietNaN (m qNaN 1.0)+ , counterexample "(1.0,qNaN)" $ isQuietNaN (m 1.0 qNaN)+ ]++prop_maximum :: RealFloatNaN a => Proxy a -> (a -> a -> a) -> Property+prop_maximum _ m =+ let sNaN = setPayloadSignaling 1+ qNaN = setPayload 1+ in conjoin+ [ counterexample "(1,3)" $ m 1 3 `sameFloatP` 3+ , counterexample "(1,-1)" $ m 1 (-1) `sameFloatP` 1+ , counterexample "(0,0)" $ m 0 0 `sameFloatP` 0+ , counterexample "(0,-0)" $ m 0 (-0) `sameFloatP` 0+ , counterexample "(-0,0)" $ m (-0) 0 `sameFloatP` 0+ , counterexample "(-0,-0)" $ m (-0) (-0) `sameFloatP` (-0)+ , counterexample "(sNaN,sNaN)" $ isQuietNaN (m sNaN sNaN)+ , counterexample "(sNaN,qNaN)" $ isQuietNaN (m sNaN qNaN)+ , counterexample "(qNaN,sNaN)" $ isQuietNaN (m qNaN sNaN)+ , counterexample "(qNaN,qNaN)" $ isQuietNaN (m qNaN qNaN)+ , counterexample "(sNaN,1.0)" $ isQuietNaN (m sNaN 1.0)+ , counterexample "(1.0,sNaN)" $ isQuietNaN (m 1.0 sNaN)+ , counterexample "(qNaN,1.0)" $ isQuietNaN (m qNaN 1.0)+ , counterexample "(1.0,qNaN)" $ isQuietNaN (m 1.0 qNaN)+ ]++prop_minimumNumber :: RealFloatNaN a => Proxy a -> (a -> a -> a) -> Property+prop_minimumNumber _ m =+ let sNaN = setPayloadSignaling 1+ qNaN = setPayload 1+ in conjoin+ [ counterexample "(1,3)" $ m 1 3 `sameFloatP` 1+ , counterexample "(1,-1)" $ m 1 (-1) `sameFloatP` (-1)+ , counterexample "(0,0)" $ m 0 0 `sameFloatP` 0+ , counterexample "(0,-0)" $ m 0 (-0) `sameFloatP` (-0)+ , counterexample "(-0,0)" $ m (-0) 0 `sameFloatP` (-0)+ , counterexample "(-0,-0)" $ m (-0) (-0) `sameFloatP` (-0)+ , counterexample "(sNaN,sNaN)" $ isQuietNaN (m sNaN sNaN)+ , counterexample "(sNaN,qNaN)" $ isQuietNaN (m sNaN qNaN)+ , counterexample "(qNaN,sNaN)" $ isQuietNaN (m qNaN sNaN)+ , counterexample "(qNaN,qNaN)" $ isQuietNaN (m qNaN qNaN)+ , counterexample "(sNaN,1.0)" $ m sNaN 1.0 `sameFloatP` 1.0+ , counterexample "(1.0,sNaN)" $ m 1.0 sNaN `sameFloatP` 1.0+ , counterexample "(qNaN,1.0)" $ m qNaN 1.0 `sameFloatP` 1.0+ , counterexample "(1.0,qNaN)" $ m 1.0 qNaN `sameFloatP` 1.0+ ]++prop_maximumNumber :: RealFloatNaN a => Proxy a -> (a -> a -> a) -> Property+prop_maximumNumber _ m =+ let sNaN = setPayloadSignaling 1+ qNaN = setPayload 1+ in conjoin+ [ counterexample "(1,3)" $ m 1 3 `sameFloatP` 3+ , counterexample "(1,-1)" $ m 1 (-1) `sameFloatP` 1+ , counterexample "(0,0)" $ m 0 0 `sameFloatP` 0+ , counterexample "(0,-0)" $ m 0 (-0) `sameFloatP` 0+ , counterexample "(-0,0)" $ m (-0) 0 `sameFloatP` 0+ , counterexample "(-0,-0)" $ m (-0) (-0) `sameFloatP` (-0)+ , counterexample "(sNaN,sNaN)" $ isQuietNaN (m sNaN sNaN)+ , counterexample "(sNaN,qNaN)" $ isQuietNaN (m sNaN qNaN)+ , counterexample "(qNaN,sNaN)" $ isQuietNaN (m qNaN sNaN)+ , counterexample "(qNaN,qNaN)" $ isQuietNaN (m qNaN qNaN)+ , counterexample "(sNaN,1.0)" $ m sNaN 1.0 `sameFloatP` 1.0+ , counterexample "(1.0,sNaN)" $ m 1.0 sNaN `sameFloatP` 1.0+ , counterexample "(qNaN,1.0)" $ m qNaN 1.0 `sameFloatP` 1.0+ , counterexample "(1.0,qNaN)" $ m 1.0 qNaN `sameFloatP` 1.0+ ]++{-# NOINLINE spec #-}+spec :: Spec+spec = do+ describe "Float" $ do+ let proxy :: Proxy Float+ proxy = Proxy+ prop "minimum'" $ prop_minimum proxy minimum'+ prop "minimum' (generic)" $ prop_minimum proxy (coerce (minimum' :: Identity Float -> Identity Float -> Identity Float))+ prop "minimumFloat" $ prop_minimum proxy minimumFloat+ prop "minimumNumber" $ prop_minimumNumber proxy minimumNumber+ prop "minimumNumber (generic)" $ prop_minimumNumber proxy (coerce (minimumNumber :: Identity Float -> Identity Float -> Identity Float))+ prop "minimumNumberFloat" $ prop_minimumNumber proxy minimumNumberFloat+ prop "maximum'" $ prop_maximum proxy maximum'+ prop "maximum' (generic)" $ prop_maximum proxy (coerce (maximum' :: Identity Float -> Identity Float -> Identity Float))+ prop "maximumFloat" $ prop_maximum proxy maximumFloat+ prop "maximumNumber" $ prop_maximumNumber proxy maximumNumber+ prop "maximumNumber (generic)" $ prop_maximumNumber proxy (coerce (maximumNumber :: Identity Float -> Identity Float -> Identity Float))+ prop "maximumNumberFloat" $ prop_maximumNumber proxy maximumNumberFloat+ describe "Double" $ do+ let proxy :: Proxy Double+ proxy = Proxy+ prop "minimum'" $ prop_minimum proxy minimum'+ prop "minimum' (generic)" $ prop_minimum proxy (coerce (minimum' :: Identity Double -> Identity Double -> Identity Double))+ prop "minimumDouble" $ prop_minimum proxy minimumDouble+ prop "minimumNumber" $ prop_minimumNumber proxy minimumNumber+ prop "minimumNumber (generic)" $ prop_minimumNumber proxy (coerce (minimumNumber :: Identity Double -> Identity Double -> Identity Double))+ prop "minimumNumberDouble" $ prop_minimumNumber proxy minimumNumberDouble+ prop "maximum'" $ prop_maximum proxy maximum'+ prop "maximum' (generic)" $ prop_maximum proxy (coerce (maximum' :: Identity Double -> Identity Double -> Identity Double))+ prop "maximumDouble" $ prop_maximum proxy maximumDouble+ prop "maximumNumber" $ prop_maximumNumber proxy maximumNumber+ prop "maximumNumber (generic)" $ prop_maximumNumber proxy (coerce (maximumNumber :: Identity Double -> Identity Double -> Identity Double))+ prop "maximumNumberDouble" $ prop_maximumNumber proxy maximumNumberDouble
+ test/NaNSpec.hs view
@@ -0,0 +1,166 @@+{-# LANGUAGE HexFloatLiterals #-}+module NaNSpec where+import Data.Proxy+import Numeric.Floating.IEEE hiding (classify, compareByTotalOrder,+ isSignMinus)+import Numeric.Floating.IEEE.NaN+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import Util++default ()++prop_copySign :: (RealFloatNaN a) => Proxy a -> a -> a -> Property+prop_copySign _ x y = let x' = copySign x y+ in isSignMinus x' === isSignMinus y++prop_isSignMinus :: (RealFloatNaN a) => Proxy a -> a -> Property+prop_isSignMinus _ x = isSignMinus (negate x) === not (isSignMinus x)++prop_isSignaling :: (RealFloatNaN a) => Proxy a -> Bool+prop_isSignaling proxy = let nan = (0 / 0) `asProxyTypeOf` proxy+ -- common floating-point operations should generate a quiet NaN+ in not (isSignaling nan)++prop_setPayload_getPayload :: (RealFloatNaN a) => Proxy a -> Property+prop_setPayload_getPayload proxy =+ let nan = (0 / 0) `asProxyTypeOf` proxy+ nan2 = setPayload (getPayload nan)+ in classify nan2 /= PositiveZero ==> compareByTotalOrder (abs nan) nan2 === EQ++prop_setPayload :: (RealFloatNaN a, Show a) => Proxy a -> a -> Property+prop_setPayload _ payload =+ let snan = setPayload payload+ in classify snan === PositiveZero .||. (not (isSignaling snan) .&&. classify snan === QuietNaN)++prop_setPayloadSignaling :: (RealFloatNaN a, Show a) => Proxy a -> a -> Property+prop_setPayloadSignaling _ payload =+ let snan = setPayloadSignaling payload+ in classify snan === PositiveZero .||. (isSignaling snan .&&. classify snan === SignalingNaN)++prop_classify :: (RealFloatNaN a, Show a) => Proxy a -> a -> Property+prop_classify _ x = conjoin+ [ counterexample "NegativeInfinity" $ (c == NegativeInfinity) === (x < 0 && isInfinite x)+ , counterexample "NegativeNormal" $ (c == NegativeNormal) === (x < 0 && isNormal x)+ , counterexample "NegativeSubnormal" $ (c == NegativeSubnormal) === (x < 0 && isDenormalized x)+ , counterexample "NegativeZero" $ (c == NegativeZero) === (isNegativeZero x)+ , counterexample "PositiveZero" $ (c == PositiveZero) === (x == 0 && not (isNegativeZero x))+ , counterexample "PositiveSubnormal" $ (c == PositiveSubnormal) === (x > 0 && isDenormalized x)+ , counterexample "PositiveNormal" $ (c == PositiveNormal) === (x > 0 && isNormal x)+ , counterexample "PositiveInfinity" $ (c == PositiveInfinity) === (x > 0 && isInfinite x)+ , counterexample "isNaN" $ isNaN x === (c == SignalingNaN || c == QuietNaN)+ , counterexample "isSignaling" $ isSignaling x === (c == SignalingNaN)+ , counterexample "isSignaling implies isNaN" $ if isSignaling x then isNaN x else True+ , counterexample "isInfinite" $ isInfinite x === (c == NegativeInfinity || c == PositiveInfinity)+ , counterexample "isNormal" $ isNormal x === (c == NegativeNormal || c == PositiveNormal)+ , counterexample "isDenormalized" $ isDenormalized x === (c == NegativeSubnormal || c == PositiveSubnormal)+ , counterexample "isZero" $ isZero x === (c == NegativeZero || c == PositiveZero)+ , counterexample "isFinite" $ isFinite x === (c `elem` [NegativeNormal, NegativeSubnormal, NegativeZero, PositiveZero, PositiveSubnormal, PositiveNormal])+ , counterexample "isSignMinus" $ if isSignMinus x then+ c `elem` [NegativeInfinity, NegativeNormal, NegativeSubnormal, NegativeZero, QuietNaN, SignalingNaN]+ else+ c `elem` [PositiveInfinity, PositiveNormal, PositiveSubnormal, PositiveZero, QuietNaN, SignalingNaN]+ -- , counterexample "class method" $ classify x === classifyDefault x+ ]+ where c = classify x+{-# SPECIALIZE prop_classify :: Proxy Float -> Float -> Property, Proxy Double -> Double -> Property #-}++isQuietNaN :: (RealFloatNaN a) => a -> Bool+isQuietNaN x = isNaN x && not (isSignaling x)++prop_signalingNaN :: (RealFloatNaN a, Show a) => Proxy a -> Property+prop_signalingNaN proxy =+ let snan = setPayloadSignaling 123 `asProxyTypeOf` proxy -- Assume 123 is a valid payload+ qnan = setPayload 123 `asProxyTypeOf` proxy -- Assume 123 is a valid payload+ in conjoin+ [ counterexample "setPayloadSignaling produces a signaling NaN" $ isSignaling snan+ , counterexample "round'" $ isQuietNaN (round' snan)+ , counterexample "roundAway'" $ isQuietNaN (roundAway' snan)+ , counterexample "truncate'" $ isQuietNaN (truncate' snan)+ , counterexample "ceiling'" $ isQuietNaN (ceiling' snan)+ , counterexample "floor'" $ isQuietNaN (floor' snan)+ , counterexample "nextUp" $ isQuietNaN (nextUp snan)+ , counterexample "nextDown" $ isQuietNaN (nextDown snan)+ , counterexample "nextTowardZero" $ isQuietNaN (nextTowardZero snan)+ -- , counterexample "remainder" $ isQuietNaN (remainder snan snan)+ -- , counterexample "scaleFloat" $ isQuietNaN (scaleFloat 1 snan)+ , counterexample "+" $ isQuietNaN (snan + snan)+ , counterexample "-" $ isQuietNaN (snan - snan)+ , counterexample "*" $ isQuietNaN (snan * snan)+ , counterexample "/" $ isQuietNaN (snan / snan)+ , counterexample "sqrt" $ isQuietNaN (sqrt snan)+ , counterexample "fusedMultiplyAdd" $ isQuietNaN (fusedMultiplyAdd snan snan snan)+ , counterexample "fusedMultiplyAdd" $ isQuietNaN (fusedMultiplyAdd 0 0 snan)+ , counterexample "negate" $ isSignaling (negate snan)+ , counterexample "abs" $ isSignaling (abs snan)+ , counterexample "augmentedAddition" $ case augmentedAddition snan snan of (x, y) -> isQuietNaN x .&&. isQuietNaN y+ , counterexample "augmentedSubtraction" $ case augmentedSubtraction snan snan of (x, y) -> isQuietNaN x .&&. isQuietNaN y+ , counterexample "augmentedMultiplication" $ case augmentedMultiplication snan snan of (x, y) -> isQuietNaN x .&&. isQuietNaN y+ , counterexample "minimum" $ isQuietNaN (minimum' snan snan)+ , counterexample "minimumNumber" $ isQuietNaN (minimumNumber snan snan)+ , counterexample "maximum" $ isQuietNaN (maximum' snan snan)+ , counterexample "maximumNumber" $ isQuietNaN (maximumNumber snan snan)+ , counterexample "minimumMagnitude" $ isQuietNaN (minimumMagnitude snan snan)+ , counterexample "minimumMagnitudeNumber" $ isQuietNaN (minimumMagnitudeNumber snan snan)+ , counterexample "maximumMagnitude" $ isQuietNaN (maximumMagnitude snan snan)+ , counterexample "maximumMagnitudeNumber" $ isQuietNaN (maximumMagnitudeNumber snan snan)+ , counterexample "canonicalize" $ isQuietNaN (canonicalize snan)+ , counterexample "realFloatToFrac" $ isQuietNaN (realFloatToFrac snan `asProxyTypeOf` proxy)+ ]+{-# INLINE prop_signalingNaN #-}++prop_totalOrder :: RealFloatNaN a => Proxy a -> a -> a -> Property+prop_totalOrder proxy x y = let cmp_x_y = compareByTotalOrder x y+ cmp_y_x = compareByTotalOrder y x+ eq = equalByTotalOrder x y+ -- cmp_reference = compareByTotalOrderDefault x y+ in cmp_x_y === compare EQ cmp_y_x+ .&&. (cmp_x_y == EQ) === eq+ -- .&&. cmp_x_y === cmp_reference+ .&&. (if x < y then cmp_x_y === LT else property True)+ .&&. (if y < x then cmp_x_y === GT else property True)+ .&&. equalByTotalOrder x x+ .&&. equalByTotalOrder y y++{-# NOINLINE spec #-}+spec :: Spec+spec = do+ describe "Float" $ do+ let proxy :: Proxy Float+ proxy = Proxy+ let snan = setPayloadSignaling 123 `asProxyTypeOf` proxy -- Assume 123 is a valid payload+ prop "copySign" $ forAllFloats2 $ prop_copySign proxy+ prop "isSignMinus" $ forAllFloats $ prop_isSignMinus proxy+ prop "isSignaling" $ prop_isSignaling proxy+ prop "setPayload/getPayload" $ prop_setPayload_getPayload proxy+ prop "setPayload/0" $ prop_setPayload proxy 0+ prop "setPayload/0x1p24" $ prop_setPayload proxy 0x1p24+ prop "setPayload/Int" $ prop_setPayload proxy . (fromIntegral :: Int -> Float)+ prop "setPayloadSignaling/0" $ prop_setPayloadSignaling proxy 0+ prop "setPayloadSignaling/0x1p24" $ prop_setPayloadSignaling proxy 0x1p24+ prop "setPayloadSignaling/Int" $ prop_setPayloadSignaling proxy . (fromIntegral :: Int -> Float)+ prop "classify" $ forAllFloats $ prop_classify proxy+ prop "classify (signaling NaN)" $ prop_classify proxy (setPayloadSignaling 123)+ prop "signaling NaN propagation" $ prop_signalingNaN proxy+ prop "totalOrder" $ forAllFloats2 $ prop_totalOrder proxy+ prop "canonicalize" $ isQuietNaN (canonicalize snan)+ describe "Double" $ do+ let proxy :: Proxy Double+ proxy = Proxy+ let snan = setPayloadSignaling 123 `asProxyTypeOf` proxy -- Assume 123 is a valid payload+ prop "copySign" $ forAllFloats2 $ prop_copySign proxy+ prop "isSignMinus" $ forAllFloats $ prop_isSignMinus proxy+ prop "isSignaling" $ prop_isSignaling proxy+ prop "setPayload/getPayload" $ prop_setPayload_getPayload proxy+ prop "setPayload/0" $ prop_setPayload proxy 0+ prop "setPayload/0x1p53" $ prop_setPayload proxy 0x1p53+ prop "setPayload/Int" $ prop_setPayload proxy . (fromIntegral :: Int -> Double)+ prop "setPayloadSignaling/0" $ prop_setPayloadSignaling proxy 0+ prop "setPayloadSignaling/0x1p53" $ prop_setPayloadSignaling proxy 0x1p53+ prop "setPayloadSignaling/Int" $ prop_setPayloadSignaling proxy . (fromIntegral :: Int -> Double)+ prop "classify" $ forAllFloats $ prop_classify proxy+ prop "classify (signaling NaN)" $ prop_classify proxy (setPayloadSignaling 123)+ prop "signaling NaN propagation" $ prop_signalingNaN proxy+ prop "totalOrder" $ forAllFloats2 $ prop_totalOrder proxy+ prop "canonicalize" $ isQuietNaN (canonicalize snan)
+ test/RoundToIntegralSpec.hs view
@@ -0,0 +1,171 @@+module RoundToIntegralSpec where+import Data.Proxy+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import Util++prop_roundToIntegral :: (RealFloat a, Show a) => Proxy a -> a -> Property+prop_roundToIntegral _ x = isFinite x ==>+ let tiesToEven = round' x+ tiesToEvenInt = round x :: Integer+ tiesToAway = roundAway' x+ tiesToAwayInt = roundAway x :: Integer+ towardPositive = ceiling' x+ towardPositiveInt = ceiling x :: Integer+ towardNegative = floor' x+ towardNegativeInt = floor x :: Integer+ towardZero = truncate' x+ towardZeroInt = truncate x :: Integer+ sameInteger f i = round f === i .&&. f === fromInteger i+ in conjoin+ [ counterexample "tiesToEven" $ isFinite tiesToEven .&&. sameInteger tiesToEven tiesToEvenInt+ , counterexample "tiesToAway" $ isFinite tiesToAway .&&. sameInteger tiesToAway tiesToAwayInt+ , counterexample "towardPositive" $ isFinite towardPositive .&&. sameInteger towardPositive towardPositiveInt+ , counterexample "towardNegative" $ isFinite towardNegative .&&. sameInteger towardNegative towardNegativeInt+ , counterexample "towardZero" $ isFinite towardZero .&&. sameInteger towardZero towardZeroInt+ , counterexample "towardNegative <= original value" $ towardNegative <= x+ , counterexample "towardNegative <= tiesToEven" $ towardNegative <= tiesToEven+ , counterexample "towardNegative <= tiesToAway" $ towardNegative <= tiesToAway+ , counterexample "towardNegative <= towardPositive" $ towardNegative <= towardPositive+ , counterexample "towardNegative <= towardZero" $ towardNegative <= towardZero+ , counterexample "original value <= towardPositive" $ x <= towardPositive+ , counterexample "tiesToEven <= towardPositive" $ tiesToEven <= towardPositive+ , counterexample "tiesToAway <= towardPositive" $ tiesToAway <= towardPositive+ , counterexample "towardZero <= towardPositive" $ towardZero <= towardPositive+ , counterexample "abs towardZero <= abs (original value)" $ abs towardZero <= abs x+ , counterexample "abs towardZero <= abs tiesToEven" $ abs towardZero <= abs tiesToEven+ , counterexample "abs towardZero <= abs tiesToAway" $ abs towardZero <= abs tiesToAway+ , counterexample "abs towardZero <= abs towardPositive" $ abs towardZero <= abs towardPositive+ , counterexample "abs towardZero <= abs towardNegative" $ abs towardZero <= abs towardNegative+ ]++data RoundResult a = RoundResult { resultTiesToEven :: a+ , resultTiesToAway :: a+ , resultTowardPositive :: a+ , resultTowardNegative :: a+ , resultTowardZero :: a+ }++checkBehavior :: RealFloat a => Proxy a -> a -> RoundResult a -> RoundResult Integer -> Spec+checkBehavior _ x result resultI = do+ it "tiesToEven" $ round' x `sameFloatP` resultTiesToEven result+ it "tiesToEven (Integer)" $ round x `shouldBe` resultTiesToEven resultI+ it "tiesToAway" $ roundAway' x `sameFloatP` resultTiesToAway result+ it "tiesToAway (Integer)" $ roundAway x `shouldBe` resultTiesToAway resultI+ it "ceiling" $ ceiling' x `sameFloatP` resultTowardPositive result+ it "ceiling (Integer)" $ ceiling x `shouldBe` resultTowardPositive resultI+ it "floor" $ floor' x `sameFloatP` resultTowardNegative result+ it "floor (Integer)" $ floor x `shouldBe` resultTowardNegative resultI+ it "truncate" $ truncate' x `sameFloatP` resultTowardZero result+ it "truncate (Integer)" $ truncate x `shouldBe` resultTowardZero resultI++checkCases :: RealFloat a => Proxy a -> Spec+checkCases proxy = do+ describe "0.5" $ checkBehavior proxy 0.5+ RoundResult { resultTiesToEven = 0.0+ , resultTiesToAway = 1.0+ , resultTowardPositive = 1.0+ , resultTowardNegative = 0.0+ , resultTowardZero = 0.0+ }+ RoundResult { resultTiesToEven = 0+ , resultTiesToAway = 1+ , resultTowardPositive = 1+ , resultTowardNegative = 0+ , resultTowardZero = 0+ }+ describe "0.25" $ checkBehavior proxy 0.25+ RoundResult { resultTiesToEven = 0.0+ , resultTiesToAway = 0.0+ , resultTowardPositive = 1.0+ , resultTowardNegative = 0.0+ , resultTowardZero = 0.0+ }+ RoundResult { resultTiesToEven = 0+ , resultTiesToAway = 0+ , resultTowardPositive = 1+ , resultTowardNegative = 0+ , resultTowardZero = 0+ }+ describe "-0.25" $ checkBehavior proxy (-0.25)+ RoundResult { resultTiesToEven = -0.0+ , resultTiesToAway = -0.0+ , resultTowardPositive = -0.0+ , resultTowardNegative = -1.0+ , resultTowardZero = -0.0+ }+ RoundResult { resultTiesToEven = 0+ , resultTiesToAway = 0+ , resultTowardPositive = 0+ , resultTowardNegative = -1+ , resultTowardZero = 0+ }+ describe "-0.5" $ checkBehavior proxy (-0.5)+ RoundResult { resultTiesToEven = -0.0+ , resultTiesToAway = -1.0+ , resultTowardPositive = -0.0+ , resultTowardNegative = -1.0+ , resultTowardZero = -0.0+ }+ RoundResult { resultTiesToEven = 0+ , resultTiesToAway = -1+ , resultTowardPositive = 0+ , resultTowardNegative = -1+ , resultTowardZero = 0+ }+ describe "4.5" $ checkBehavior proxy 4.5+ RoundResult { resultTiesToEven = 4.0+ , resultTiesToAway = 5.0+ , resultTowardPositive = 5.0+ , resultTowardNegative = 4.0+ , resultTowardZero = 4.0+ }+ RoundResult { resultTiesToEven = 4+ , resultTiesToAway = 5+ , resultTowardPositive = 5+ , resultTowardNegative = 4+ , resultTowardZero = 4+ }+ describe "-5.5" $ checkBehavior proxy (-5.5)+ RoundResult { resultTiesToEven = -6.0+ , resultTiesToAway = -6.0+ , resultTowardPositive = -5.0+ , resultTowardNegative = -6.0+ , resultTowardZero = -5.0+ }+ RoundResult { resultTiesToEven = -6+ , resultTiesToAway = -6+ , resultTowardPositive = -5+ , resultTowardNegative = -6+ , resultTowardZero = -5+ }+ describe "-6.5" $ checkBehavior proxy (-6.5)+ RoundResult { resultTiesToEven = -6.0+ , resultTiesToAway = -7.0+ , resultTowardPositive = -6.0+ , resultTowardNegative = -7.0+ , resultTowardZero = -6.0+ }+ RoundResult { resultTiesToEven = -6+ , resultTiesToAway = -7+ , resultTowardPositive = -6+ , resultTowardNegative = -7+ , resultTowardZero = -6+ }++{-# NOINLINE spec #-}+spec :: Spec+spec = do+ describe "Double" $ do+ let proxy :: Proxy Double+ proxy = Proxy+ prop "roundToIntegral" $ prop_roundToIntegral proxy+ checkCases proxy+ describe "Float" $ do+ let proxy :: Proxy Double+ proxy = Proxy+ prop "roundToIntegral" $ prop_roundToIntegral proxy+ checkCases proxy
+ test/RoundingSpec.hs view
@@ -0,0 +1,241 @@+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE HexFloatLiterals #-}+{-# LANGUAGE NumericUnderscores #-}+{-# LANGUAGE RankNTypes #-}+module RoundingSpec where+import Control.Monad+import Data.Int+import Data.Proxy+import Data.Ratio+import Data.Word+import Numeric+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (classify)+import Util++newtype RoundTiesTowardZero a = RoundTiesTowardZero { roundTiesTowardZero :: a }+ deriving (Functor)++instance RoundingStrategy RoundTiesTowardZero where+ exact = RoundTiesTowardZero+ inexact o _neg _parity zero away = RoundTiesTowardZero $ case o of+ LT -> zero+ EQ -> zero+ GT -> away+ doRound _exact o _neg _parity zero away = RoundTiesTowardZero $ case o of+ LT -> zero+ EQ -> zero+ GT -> away++newtype RoundToOdd a = RoundToOdd { roundToOdd :: a }+ deriving (Functor)++instance RoundingStrategy RoundToOdd where+ exact = RoundToOdd+ inexact _o _neg parity zero away | even parity = RoundToOdd away+ | otherwise = RoundToOdd zero+ doRound exact _o _neg parity zero away | not exact && even parity = RoundToOdd away+ | otherwise = RoundToOdd zero++newtype Exactness a = Exactness { isExact :: Bool }+ deriving (Functor)++instance RoundingStrategy Exactness where+ exact _ = Exactness True+ inexact _o _neg _parity _zero _away = Exactness False+ doRound exact _o _neg _parity _zero _away = Exactness exact++prop_fromIntegerR_vs_fromIntegralR :: (RealFloat a, RoundingStrategy f, Integral i) => Proxy a -> Proxy i -> (f a -> a) -> i -> Property+prop_fromIntegerR_vs_fromIntegralR _ _ f m =+ let x = f (fromIntegerR (toInteger m))+ y = f (fromIntegralR m)+ in x `sameFloatP` y++prop_fromIntegerR_vs_fromRationalR :: (RealFloat a, RoundingStrategy f) => Proxy a -> (f a -> a) -> Integer -> Property+prop_fromIntegerR_vs_fromRationalR _ f m =+ let x = f (fromIntegerR m)+ y = f (fromRationalR (m % 1))+ in x `sameFloatP` y++prop_fromIntegerR_vs_encodeFloatR :: (RealFloat a, RoundingStrategy f) => Proxy a -> (f a -> a) -> Integer -> NonNegative Int -> Property+prop_fromIntegerR_vs_encodeFloatR _ f m (NonNegative k) =+ let x = f (fromIntegerR m)+ y = f (encodeFloatR (m * floatRadix x ^ k) (-k))+ in x `sameFloatP` y++prop_fromRationalR_vs_encodeFloatR :: (RealFloat a, RoundingStrategy f) => Proxy a -> (f a -> a) -> Integer -> Int -> Property+prop_fromRationalR_vs_encodeFloatR _ f m k =+ let x = f (fromRationalR (fromInteger m * fromInteger (floatRadix x) ^^ k))+ y = f (encodeFloatR m k)+ in x `sameFloatP` y++prop_fromRationalR_vs_fromRational :: RealFloat a => Proxy a -> Rational -> Property+prop_fromRationalR_vs_fromRational proxy q =+ let x = roundTiesToEven (fromRationalR q) `asProxyTypeOf` proxy+ y = fromRational q `asProxyTypeOf` proxy+ in x `sameFloatP` y++prop_scaleFloatR_vs_fromRationalR :: (RealFloat a, RoundingStrategy f) => Proxy a -> (f a -> a) -> Int -> a -> Property+prop_scaleFloatR_vs_fromRationalR proxy f e x = isFinite x && not (isNegativeZero x) ==>+ let base = floatRadix x+ y = f (scaleFloatR e x)+ z = f (fromRationalR (toRational x * fromInteger base^^e))+ in y `sameFloatP` z++prop_scaleFloatR_vs_encodeFloatR :: (RealFloat a, RoundingStrategy f) => Proxy a -> (f a -> a) -> Int -> a -> Property+prop_scaleFloatR_vs_encodeFloatR proxy f e x = isFinite x && not (isNegativeZero x) ==>+ let base = floatRadix x+ (m,n) = decodeFloat x+ y = f (scaleFloatR e x)+ z = f (encodeFloatR m (n + e))+ in y `sameFloatP` z++prop_encodeFloatR_roundtrip :: (RealFloat a, RoundingStrategy f) => Proxy a -> a -> (f a -> a) -> Property+prop_encodeFloatR_roundtrip proxy x rounding = isFinite x && not (isNegativeZero x) ==>+ let (m,n) = decodeFloat x+ in rounding (encodeFloatR m n) `sameFloatP` x++prop_order :: RealFloat a => Proxy a -> (forall f. RoundingStrategy f => f a) -> Property+prop_order _ result =+ let tiesToEven = roundTiesToEven result+ tiesToAway = roundTiesToAway result+ tiesTowardZero = roundTiesTowardZero result+ up = roundTowardPositive result+ down = roundTowardNegative result+ zero = roundTowardZero result+ toOdd = roundToOdd result+ in if isExact result then+ counterexample "exact case" $ conjoin+ [ counterexample "tiesToAway == tiesToEven" $ tiesToAway `sameFloatP` tiesToEven+ , counterexample "tiesTowardZero == tiesToEven" $ tiesTowardZero `sameFloatP` tiesToEven+ , counterexample "upward == tiesToEven" $ up `sameFloatP` tiesToEven+ , counterexample "downward == tiesToEven" $ down `sameFloatP` tiesToEven+ , counterexample "towardZero == tiesToEven" $ zero `sameFloatP` tiesToEven+ , counterexample "toOdd == tiesToEven" $ toOdd `sameFloatP` tiesToEven+ ]+ else+ counterexample "inexact case" $ conjoin+ [ counterexample "down < up" $ down < up+ , counterexample "down <= tiesToEven" $ down <= tiesToEven+ , counterexample "down <= tiesToAway" $ down <= tiesToAway+ , counterexample "down <= tiesTowardZero" $ down <= tiesTowardZero+ , counterexample "down <= towardZero" $ down <= zero+ , counterexample "down <= odd" $ down <= toOdd+ , counterexample "tiesToEven <= up" $ tiesToEven <= up+ , counterexample "tiesToAway <= up" $ tiesToAway <= up+ , counterexample "tiesTowardZero <= up" $ tiesTowardZero <= up+ , counterexample "towardZero <= up" $ zero <= up+ , counterexample "odd <= up" $ toOdd <= up+ , counterexample "nextUp down == up" $ nextUp down `sameFloatP` up+ , counterexample "down == nextDown up" $ down `sameFloatP` nextDown up+ , counterexample "abs towardZero < max (abs down) (abs up)" $ abs zero < max (abs down) (abs up)+ , counterexample "not (isMantissaEven toOdd)" $ not (isMantissaEven toOdd)+ ]++prop_addToOdd :: RealFloat a => Proxy a -> a -> a -> Property+prop_addToOdd _ x y = isFinite x && isFinite y && isFinite (x + y) ==>+ let z = addToOdd x y+ w = if x == 0 && y == 0 then+ x + y+ else+ roundToOdd (fromRationalR (toRational x + toRational y))+ in z `sameFloatP` w++eachStrategy :: Testable prop => (forall f. RoundingStrategy f => (f a -> a) -> prop) -> Property+eachStrategy p = conjoin+ [ counterexample "roundTiesToEven" (p roundTiesToEven)+ , counterexample "roundTiesToAway" (p roundTiesToAway)+ , counterexample "roundTiesTowardZero" (p roundTiesTowardZero)+ , counterexample "roundTowardPositive" (p roundTowardPositive)+ , counterexample "roundTowardNegative" (p roundTowardNegative)+ , counterexample "roundTowardZero" (p roundTowardZero)+ , counterexample "roundToOdd" (p roundToOdd)+ ]++testUnary :: RealFloat b => (a -> b) -> [(String, a, b)] -> Property+testUnary f cases = conjoin+ [ counterexample t $ f a `sameFloatP` result+ | (t,a,result) <- cases+ ]++{-# NOINLINE spec #-}+spec :: Spec+spec = do+ describe "Double" $ do+ let proxy :: Proxy Double+ proxy = Proxy+ prop "fromIntegerR vs fromIntegralR" $ eachStrategy (prop_fromIntegerR_vs_fromIntegralR proxy (Proxy :: Proxy Int))+ prop "fromIntegerR vs fromIntegralR" $ eachStrategy (prop_fromIntegerR_vs_fromIntegralR proxy (Proxy :: Proxy Word64))+ prop "fromIntegerR vs fromRationalR" $ eachStrategy (prop_fromIntegerR_vs_fromRationalR proxy)+ prop "fromIntegerR vs encodeFloatR" $ eachStrategy (prop_fromIntegerR_vs_encodeFloatR proxy)+ prop "fromRationalR vs encodeFloatR" $ eachStrategy (prop_fromRationalR_vs_encodeFloatR proxy)+ prop "fromRationalR vs fromRational" $ prop_fromRationalR_vs_fromRational proxy+ prop "scaleFloatR vs fromRationalR" $ eachStrategy (prop_scaleFloatR_vs_fromRationalR proxy)+ prop "scaleFloatR vs encodeFloatR" $ eachStrategy (prop_scaleFloatR_vs_encodeFloatR proxy)+ prop "result of fromIntegerR" $ \x -> prop_order proxy (fromIntegerR x)+ prop "result of fromRationalR" $ \x -> prop_order proxy (fromRationalR x)+ prop "result of encodeFloatR" $ \m k -> prop_order proxy (encodeFloatR m k)+ prop "encodeFloatR/decodeFloat" $ forAllFloats $ \x -> eachStrategy (prop_encodeFloatR_roundtrip proxy x)+ prop "addToOdd" $ forAllFloats2 $ prop_addToOdd proxy+ it "fromIntegralR/(maxBound :: Int64)" $ fromIntegralTiesToEven (maxBound :: Int64) `sameFloatP` (0x1p63 :: Double)+ it "fromIntegralR/(maxBound :: Word64)" $ fromIntegralTiesToEven (maxBound :: Word64) `sameFloatP` (0x1p64 :: Double)+ it "fromIntegralR/(0xffff_ffff_ffff_fc00 :: Word64)" $ fromIntegralTiesToEven (0xffff_ffff_ffff_fc00 :: Word64) `sameFloatP` (0x1p64 :: Double)++ do let cases :: [(String, Rational, Double)]+ cases = [ let t = 11435996997111233 % 1660860084017817297360368008619370227400073727045418226348482155039064904553019973177107435363660614816513137110180404061646380785658477636245443559462428597275694780106044074992747404797486457853074429979899122551795724461450521406238742712434733270295344316890429535153317233021396948961884411359194146958100478088711873454042107514097515809485603670823814576138204139943337375836756405167181947093525325738801370702465460537395969617160395178613194019276299200610817420725783045671692771793360418111369105879747924354309959938911042057102540038489527102833880604228417018090258140649799612644290906038462100262234760641844967425501906703079079531111883261520094019262965803907605528809355522427428605283171700681998722400652411744851193007546978988038363226440325125816593274436339451950472293881264365176866099134907912252035904613356400091473040550399623768278773198402959131216609632370028659088546103031543716668650443675061896807069455112892464207615075528889823150217287305246018046657536654015550308954692439217754082060020956581265580805928178408368880094563736441111304424147055967579092700683418565515720301167266647150173895623838705449444022652355565392171702345881427096566633769494957447420015296687812138177576466001557317056675111027221005969582058022899529333118501380166134607864676483828739116173461269178580186257490266486677839206143742952162243351494227378653938710593503436239164822135914417914190306326552366654989657047816161866088059657348484650208804648917587381647596311004763609009433923628807524614747087370907674848755682961586688315674280522685036343187379852233394640325214899081294504832057011229815959420037782873168548447320460610743719611348921807017679017481761450571271353121504538616599488981500090579800223920074190259243090373197975821900780700994983554220578443939059789455319157185586612934190457927556018513712845810999355955231047286405348577289698269949529315641676747401507179920683528906096656269865346604825245447613068900673228624597839983919623198678889246997823457303425538250074268963180449541718530763258905302809600007299944411273014987193501682320824514373693134713866527503191580073279143086275003713746690240664855814859039455876481938038239569220725678019039050480876390746831297254406921453270519267262507843820232264191737352673268925464180832643899338691638282218257385606257475776691059059255302303114445822454278600219173720763694867106875068457113502491500388073504656799152135200251581518256215104764399515301707283274754264151543807825364161450197108879883727387093477427770004354318482968886709591946818257266432018518668134005188950818936559490195651342132066807456183872268255020846456930757669626740368630473237568715189840731662896998327481598778409201158765383448914364093919235518500273991995313625439096723754872384506907411868540620101022260019920486730850164320257564380330469491975531388141021789624314602105976973474026654086478535953344727481858929880747213733511596028875230104172769919771254751076195477658238344543363620834339799493240979523682870604654974849807411458413970564431884272290785767041903645182376449237883070663106400054251118347592277048642471665850924191188071391188795617326279324544211665128645710824853683627205877921300176381646070686087465015189344127757236896514029243563980383479813573936253276755173117350734421524872428449939741005549450504235910411855579757233304417120352975265436957913138078770206426593236938077956476982936814123774536338991098653758589247087558267603817517200390767537146533680698510122118199916051754470078537238491169553359792229740918073741384817330552101229726860709591659018519799482781149265985004923079601834415995143876094479546972166462535851542643215260243141498224867577987788423766186348687317679115018525104558716345706749942890553933642650461618674671154556006755314390616549147093936804564986443463961450438220362152406762190515061823854149008437459939334182355104574856378203866544401000626238988568308977840150220171310744124565624246900478266970170867838195019777202546995092582751359519995005632488038116545366585729919917509021256056617930088346818246573722242278351202844467835535078947626466439417333092836098812554627989117607752545931702872303942308409649609541986879146218441452717032910609434831215306455063339901653706420056993908069607050862479753834786944380384807128177688568002157423912412284060326246610084680338789899668589451070097117651298167403754077408452603311106679269461516981669288627428528214985766284440659354036167812199161489266566736683801438390018297720643002232031866138861219487931264851019071593248506045777980832084764662336685649221969889059160428833116253588012280798203184065757408956940520408997184425057879282238950799253433771440870506862193781343867894277617920304811869690899227908547152726181311361021942187101849272547858549820527191290014454746676308089316055111376988866151778299477802926255655650697478276694050540148953139848340830296268047498950+ in ('(' : shows t ")", t, 0.0)+ ]+ prop "roundTiesToEven" $ testUnary (roundTiesToEven . fromRationalR) cases++ let cases :: [(String, Rational, Double)]+ cases = [("0x1.ffff_ffff_ffff_f8p1023", 0x1.ffff_ffff_ffff_f8p1023, maxFinite)+ ,("(0x1.ffff_ffff_ffff_f8p1023 + 1/723)", 0x1.ffff_ffff_ffff_f8p1023 + 1/723, 1/0)+ ,("(0x1.ffff_ffff_ffff_f8p1023 - 1/255)", 0x1.ffff_ffff_ffff_f8p1023 - 1/255, maxFinite)+ ,("0xdead_beef.8p-1074", 0xdead_beef.8p-1074, 0xdead_beefp-1074)+ ,("0xdead_beef.9p-1074", 0xdead_beef.9p-1074, 0xdead_bef0p-1074)+ ,("-0xdead_beef.7p-1074", -0xdead_beef.7p-1074, -0xdead_beefp-1074)+ ,("-0x0.8p-1074", -0x0.8p-1074, -0)+ ,("-0x0.80007p-1074", -0x0.80007p-1074, -0x1p-1074)+ ]+ prop "roundTiesTowardZero" $ testUnary (roundTiesTowardZero . fromRationalR) cases++ describe "Float" $ do+ let proxy :: Proxy Float+ proxy = Proxy+ prop "fromIntegerR vs fromRationalR" $ eachStrategy (prop_fromIntegerR_vs_fromRationalR proxy)+ prop "fromIntegerR vs encodeFloatR" $ eachStrategy (prop_fromIntegerR_vs_encodeFloatR proxy)+ prop "fromRationalR vs encodeFloatR" $ eachStrategy (prop_fromRationalR_vs_encodeFloatR proxy)+ prop "fromRationalR vs fromRational" $ prop_fromRationalR_vs_fromRational proxy+ prop "scaleFloatR vs fromRationalR" $ eachStrategy (prop_scaleFloatR_vs_fromRationalR proxy)+ prop "scaleFloatR vs encodeFloatR" $ eachStrategy (prop_scaleFloatR_vs_encodeFloatR proxy)+ prop "result of fromIntegerR" $ \x -> prop_order proxy (fromIntegerR x)+ prop "result of fromRationalR" $ \x -> prop_order proxy (fromRationalR x)+ prop "result of encodeFloatR" $ \m k -> prop_order proxy (encodeFloatR m k)+ prop "encodeFloatR/decodeFloat" $ forAllFloats $ \x -> eachStrategy (prop_encodeFloatR_roundtrip proxy x)+ prop "addToOdd" $ forAllFloats2 $ prop_addToOdd proxy+ it "fromIntegralR/(maxBound :: Int32)" $ fromIntegralTiesToEven (maxBound :: Int32) `sameFloatP` (0x1p31 :: Float)+ it "fromIntegralR/(maxBound :: Word32)" $ fromIntegralTiesToEven (maxBound :: Word32) `sameFloatP` (0x1p32 :: Float)+ it "fromIntegralR/(maxBound :: Int64)" $ fromIntegralTiesToEven (maxBound :: Int64) `sameFloatP` (0x1p63 :: Float)+ it "fromIntegralR/(maxBound :: Word64)" $ fromIntegralTiesToEven (maxBound :: Word64) `sameFloatP` (0x1p64 :: Float)+ it "fromIntegralR/(0xffff_ff80_0000_0000 :: Word64)" $ fromIntegralTiesToEven (0xffff_ff80_0000_0000 :: Word64) `sameFloatP` (0x1p64 :: Float)++ do let cases :: [(String, Rational, Float)]+ cases = [ let t = 20113311130255 % 822127761653273855988822146978202976557090789271144163906483851513046701868339517444102604474616762490976436939594169664101896669409817473587913461546435532885567073887954501607977104895740769882295378286300234464764201845440572849224022844453347299057834829757872072616746710668820893729486742297607776797874+ in ('(' : shows t ")", t, 0.0)+ ]+ prop "roundTiesToEven" $ testUnary (roundTiesToEven . fromRationalR) cases++ do let cases :: [(String, Rational, Float)]+ cases = [ ("0x1.ffff_ffp127", 0x1.ffff_ffp127, maxFinite)+ , ("(0x1.ffff_ffp127 + 1/723)", 0x1.ffff_ffp127 + 1/723, 1/0)+ , ("(0x1.ffff_ffp127 - 1/255)", 0x1.ffff_ffp127 - 1/255, maxFinite)+ , ("0xbeef.8p-149", 0xbeef.8p-149, 0xbeefp-149)+ , ("0xbeef.9p-149", 0xbeef.9p-149, 0xbef0p-149)+ , ("-0xbeef.7p-149", -0xbeef.7p-149, -0xbeefp-149)+ , ("-0x0.8p-149", -0x0.8p-149, -0)+ , ("-0x0.80007p-149", -0x0.80007p-149, -0x1p-149)+ ]+ prop "roundTiesTowardZero" $ testUnary (roundTiesTowardZero . fromRationalR) cases
+ test/Spec.hs view
@@ -0,0 +1,54 @@+{-# LANGUAGE CPP #-}+import qualified AugmentedArithSpec+import qualified ClassificationSpec+import qualified FMASpec+import qualified IntegerInternalsSpec+import qualified MinMaxSpec+import qualified NaNSpec+import qualified NextFloatSpec+import qualified RoundingSpec+import qualified RoundToIntegralSpec+import System.Environment (getArgs, withArgs)+import Test.Hspec hiding (hspec)+import Test.Hspec.Core.Runner hiding (hspec)+import qualified TwoSumSpec+#if defined(USE_HALF)+import qualified HalfSpec+#endif+#if defined(USE_FLOAT128)+import qualified Float128Spec+#endif++-- "Extra" tests are not run by default; set --skip "***" to run them.+myFilter :: Path -> Bool+myFilter (groups, _description) = "Extra" `elem` groups++withDefaultFilter :: Config -> Config+withDefaultFilter config@(Config { configSkipPredicate = Nothing }) = config { configSkipPredicate = Just myFilter }+withDefaultFilter config = config++hspec :: Spec -> IO ()+hspec spec =+ getArgs+ >>= readConfig defaultConfig+ >>= withArgs [] . runSpec spec . withDefaultFilter+ >>= evaluateSummary++main :: IO ()+main = hspec $ do+ describe "Classification" ClassificationSpec.spec+ describe "TwoSum" TwoSumSpec.spec+ describe "FMA" FMASpec.spec+ describe "IntegerInternals" IntegerInternalsSpec.spec+ describe "NextFloat" NextFloatSpec.spec+ describe "AugmentedArith" AugmentedArithSpec.spec+ describe "Rounding" RoundingSpec.spec+ describe "RoundToIntegral" RoundToIntegralSpec.spec+ describe "NaN" NaNSpec.spec+ describe "MinMax" MinMaxSpec.spec+#if defined(USE_HALF)+ describe "Half" HalfSpec.spec+#endif+#if defined(USE_FLOAT128)+ describe "Float128" Float128Spec.spec+#endif
+ test/TwoSumSpec.hs view
@@ -0,0 +1,39 @@+module TwoSumSpec where+import Data.Coerce+import Data.Functor.Identity+import Data.Proxy+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck+import Util (forAllFloats2, sameFloatP)++twoProduct_generic :: RealFloat a => a -> a -> (a, a)+twoProduct_generic x y = coerce (twoProduct (Identity x) (Identity y))++prop_twoSum :: (RealFloat a, Show a) => Proxy a -> a -> a -> Property+prop_twoSum _ x y = exponent x < expMax && exponent y < expMax ==> case twoSum x y of+ (s, t) -> x + y `sameFloatP` s .&&. (isFinite x && isFinite y && isFinite s ==> isFinite t .&&. toRational x + toRational y === toRational s + toRational t)+ where (_,expMax) = floatRange x++prop_twoProduct :: (RealFloat a, Show a) => Proxy a -> (a -> a -> (a, a)) -> a -> a -> Property+prop_twoProduct _ tp x y = case tp x y of+ (s, t) -> x * y `sameFloatP` s .&&. (isFinite x && isFinite y && isFinite s ==> isFinite t .&&. fromRational (toRational x * toRational y - toRational s) === t) -- The result of twoProduct is not exact if the product underflows++{-# NOINLINE spec #-}+spec :: Spec+spec = modifyMaxSuccess (* 100) $ do+ describe "Double" $ do+ let proxy :: Proxy Double+ proxy = Proxy+ prop "twoSum" $ forAllFloats2 $ prop_twoSum proxy+ prop "twoProduct" $ forAllFloats2 $ prop_twoProduct proxy twoProduct+ prop "twoProduct_generic" $ forAllFloats2 $ prop_twoProduct proxy twoProduct_generic+ describe "Float" $ do+ let proxy :: Proxy Float+ proxy = Proxy+ prop "twoSum" $ forAllFloats2 $ prop_twoSum proxy+ prop "twoProduct" $ forAllFloats2 $ prop_twoProduct proxy twoProduct+ prop "twoProduct_generic" $ forAllFloats2 $ prop_twoProduct proxy twoProduct_generic+ prop "twoProductFloat_viaDouble" $ forAllFloats2 $ prop_twoProduct proxy twoProductFloat_viaDouble