diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,5 @@ +# Changelog for fp-ieee + +## Version 0.1.0 (2020-12-27) + +Initial release. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/benchmark/Benchmark.hs b/benchmark/Benchmark.hs new file mode 100644 --- /dev/null +++ b/benchmark/Benchmark.hs @@ -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 + ] diff --git a/cbits/canonicalize.c b/cbits/canonicalize.c new file mode 100644 --- /dev/null +++ b/cbits/canonicalize.c @@ -0,0 +1,64 @@ +#include + +#pragma STDC FENV_ACCESS ON + +#if defined(__SSE2__) + +#include + +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 diff --git a/cbits/fma.c b/cbits/fma.c new file mode 100644 --- /dev/null +++ b/cbits/fma.c @@ -0,0 +1,17 @@ +#include + +#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); +} diff --git a/cbits/half.c b/cbits/half.c new file mode 100644 --- /dev/null +++ b/cbits/half.c @@ -0,0 +1,120 @@ +#include // uint16_t +#include + +#if defined(__F16C__) // x86 F16C + +#include + +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 diff --git a/cbits/minmax.c b/cbits/minmax.c new file mode 100644 --- /dev/null +++ b/cbits/minmax.c @@ -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 diff --git a/cbits/roundeven.c b/cbits/roundeven.c new file mode 100644 --- /dev/null +++ b/cbits/roundeven.c @@ -0,0 +1,48 @@ +#include +#include + +#if defined(__SSE4_1__) // SSE 4.1 + +#include + +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 diff --git a/decimal-test/NextFloatSpec.hs b/decimal-test/NextFloatSpec.hs new file mode 100644 --- /dev/null +++ b/decimal-test/NextFloatSpec.hs @@ -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 diff --git a/decimal-test/Spec.hs b/decimal-test/Spec.hs new file mode 100644 --- /dev/null +++ b/decimal-test/Spec.hs @@ -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 diff --git a/doctests.hs b/doctests.hs new file mode 100644 --- /dev/null +++ b/doctests.hs @@ -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" + ] diff --git a/fp-ieee.cabal b/fp-ieee.cabal new file mode 100644 --- /dev/null +++ b/fp-ieee.cabal @@ -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 +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 diff --git a/src/GHC/Float/Compat.hs b/src/GHC/Float/Compat.hs new file mode 100644 --- /dev/null +++ b/src/GHC/Float/Compat.hs @@ -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 diff --git a/src/MyPrelude.hs b/src/MyPrelude.hs new file mode 100644 --- /dev/null +++ b/src/MyPrelude.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE.hs b/src/Numeric/Floating/IEEE.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE.hs @@ -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. diff --git a/src/Numeric/Floating/IEEE/Internal.hs b/src/Numeric/Floating/IEEE/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/Augmented.hs b/src/Numeric/Floating/IEEE/Internal/Augmented.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Augmented.hs @@ -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) #-} diff --git a/src/Numeric/Floating/IEEE/Internal/Base.hs b/src/Numeric/Floating/IEEE/Internal/Base.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Base.hs @@ -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) +-} diff --git a/src/Numeric/Floating/IEEE/Internal/Classify.hs b/src/Numeric/Floating/IEEE/Internal/Classify.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Classify.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/Conversion.hs b/src/Numeric/Floating/IEEE/Internal/Conversion.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Conversion.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/FMA.hs b/src/Numeric/Floating/IEEE/Internal/FMA.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/FMA.hs @@ -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 - \@. +-- +-- 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 diff --git a/src/Numeric/Floating/IEEE/Internal/Float128.hs b/src/Numeric/Floating/IEEE/Internal/Float128.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Float128.hs @@ -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 + #-} diff --git a/src/Numeric/Floating/IEEE/Internal/GenericArith.hs b/src/Numeric/Floating/IEEE/Internal/GenericArith.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/GenericArith.hs @@ -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) diff --git a/src/Numeric/Floating/IEEE/Internal/Half.hs b/src/Numeric/Floating/IEEE/Internal/Half.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Half.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs b/src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/IntegerInternals.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/MinMax.hs b/src/Numeric/Floating/IEEE/Internal/MinMax.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/MinMax.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/NaN.hs b/src/Numeric/Floating/IEEE/Internal/NaN.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/NaN.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/NextFloat.hs b/src/Numeric/Floating/IEEE/Internal/NextFloat.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/NextFloat.hs @@ -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" diff --git a/src/Numeric/Floating/IEEE/Internal/Remainder.hs b/src/Numeric/Floating/IEEE/Internal/Remainder.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Remainder.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs b/src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/RoundToIntegral.hs @@ -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 diff --git a/src/Numeric/Floating/IEEE/Internal/Rounding.hs b/src/Numeric/Floating/IEEE/Internal/Rounding.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Rounding.hs @@ -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 () diff --git a/src/Numeric/Floating/IEEE/Internal/Rounding/Common.hs b/src/Numeric/Floating/IEEE/Internal/Rounding/Common.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Rounding/Common.hs @@ -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 #-} diff --git a/src/Numeric/Floating/IEEE/Internal/Rounding/Encode.hs b/src/Numeric/Floating/IEEE/Internal/Rounding/Encode.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Rounding/Encode.hs @@ -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 + #-} diff --git a/src/Numeric/Floating/IEEE/Internal/Rounding/Integral.hs b/src/Numeric/Floating/IEEE/Internal/Rounding/Integral.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Rounding/Integral.hs @@ -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)) + #-} diff --git a/src/Numeric/Floating/IEEE/Internal/Rounding/Rational.hs b/src/Numeric/Floating/IEEE/Internal/Rounding/Rational.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/Internal/Rounding/Rational.hs @@ -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)) + #-} diff --git a/src/Numeric/Floating/IEEE/NaN.hs b/src/Numeric/Floating/IEEE/NaN.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Floating/IEEE/NaN.hs @@ -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 diff --git a/test/AugmentedArithSpec.hs b/test/AugmentedArithSpec.hs new file mode 100644 --- /dev/null +++ b/test/AugmentedArithSpec.hs @@ -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 diff --git a/test/ClassificationSpec.hs b/test/ClassificationSpec.hs new file mode 100644 --- /dev/null +++ b/test/ClassificationSpec.hs @@ -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) diff --git a/test/FMASpec.hs b/test/FMASpec.hs new file mode 100644 --- /dev/null +++ b/test/FMASpec.hs @@ -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 #-} diff --git a/test/Float128Spec.hs b/test/Float128Spec.hs new file mode 100644 --- /dev/null +++ b/test/Float128Spec.hs @@ -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 diff --git a/test/HalfSpec.hs b/test/HalfSpec.hs new file mode 100644 --- /dev/null +++ b/test/HalfSpec.hs @@ -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." diff --git a/test/IntegerInternalsSpec.hs b/test/IntegerInternalsSpec.hs new file mode 100644 --- /dev/null +++ b/test/IntegerInternalsSpec.hs @@ -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)) diff --git a/test/MinMaxSpec.hs b/test/MinMaxSpec.hs new file mode 100644 --- /dev/null +++ b/test/MinMaxSpec.hs @@ -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 diff --git a/test/NaNSpec.hs b/test/NaNSpec.hs new file mode 100644 --- /dev/null +++ b/test/NaNSpec.hs @@ -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) diff --git a/test/RoundToIntegralSpec.hs b/test/RoundToIntegralSpec.hs new file mode 100644 --- /dev/null +++ b/test/RoundToIntegralSpec.hs @@ -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 diff --git a/test/RoundingSpec.hs b/test/RoundingSpec.hs new file mode 100644 --- /dev/null +++ b/test/RoundingSpec.hs @@ -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 diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -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 diff --git a/test/TwoSumSpec.hs b/test/TwoSumSpec.hs new file mode 100644 --- /dev/null +++ b/test/TwoSumSpec.hs @@ -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