rounded-hw 0.1.0.0 → 0.2.0
raw patch · 24 files changed
+292/−908 lines, 24 filesdep +fp-ieeedep −integer-logarithms
Dependencies added: fp-ieee
Dependencies removed: integer-logarithms
Files
- ChangeLog.md +6/−0
- benchmark/Benchmark.hs +0/−54
- benchmark/Conversion.hs +87/−58
- rounded-hw.cabal +11/−13
- src/Numeric/Rounded/Hardware/Backend/C.hs +10/−46
- src/Numeric/Rounded/Hardware/Backend/Default.hs +13/−12
- src/Numeric/Rounded/Hardware/Backend/FastFFI.hs +1/−0
- src/Numeric/Rounded/Hardware/Backend/Float128.hs +2/−3
- src/Numeric/Rounded/Hardware/Backend/ViaRational.hs +10/−11
- src/Numeric/Rounded/Hardware/Backend/X87LongDouble.hs +2/−3
- src/Numeric/Rounded/Hardware/Internal.hs +0/−1
- src/Numeric/Rounded/Hardware/Internal/Class.hs +6/−0
- src/Numeric/Rounded/Hardware/Internal/Conversion.hs +34/−193
- src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs +1/−307
- src/Numeric/Rounded/Hardware/Internal/RoundedResult.hs +0/−54
- src/Numeric/Rounded/Hardware/Internal/Show.hs +47/−41
- src/Numeric/Rounded/Hardware/Interval.hs +7/−0
- src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs +7/−0
- test/FloatUtilSpec.hs +0/−75
- test/FromIntegerSpec.hs +15/−15
- test/FromRationalSpec.hs +8/−8
- test/ShowFloatSpec.hs +25/−9
- test/Spec.hs +0/−2
- test/X87LongDoubleSpec.hs +0/−3
ChangeLog.md view
@@ -1,5 +1,11 @@ # Changelog for rounded-hw +## 0.2.0 (2020-12-27)++* Some functionality was moved to fp-ieee.+* Fix roundedFusedMultiplyAdd of ViaRational.+* Fix showFFloatRn.+ ## 0.1.0.0 (2020-06-23) * Initial release.
benchmark/Benchmark.hs view
@@ -13,38 +13,26 @@ import qualified Data.Vector.Unboxed as VU import Gauge.Main import IGA-import Numeric import Numeric.Rounded.Hardware.Internal import Numeric.Rounded.Hardware.Interval import Numeric.Rounded.Hardware.Interval.Class (makeInterval) import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as NE import qualified Numeric.Rounded.Hardware.Vector.Unboxed as RVU -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)- main :: IO () main = defaultMain@@ -114,8 +102,6 @@ , bench "mul" $ nf (uncurry (*)) (iv1, iv2) , bench "div" $ nf (uncurry (/)) (iv1, iv2) , bench "sqrt" $ nf sqrt iv1- , bench "fromInteger" $ nf (fromInteger :: Integer -> Interval Double) (2^60 + 1)- , bench "fromIntegral/Int64" $ nf (fromIntegral :: Int64 -> Interval Double) (2^60 + 1) ] , let vec :: V.Vector (Interval Double) vec = V.generate 100000 $ \i -> fromRational (1 % (1 + fromIntegral i))@@ -128,46 +114,6 @@ , bench "NE.exp" $ nf exp (0.3 :: NE.Interval Double) , bench "sin" $ nf sin (7.3 :: Interval Double) , bench "NE.sin" $ nf sin (7.3 :: NE.Interval Double)- ]- , 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 ]- ] ] , bgroup "FMA" [ let arg = (1.0, 2.0, 3.0) :: (Double, Double, Double)
benchmark/Conversion.hs view
@@ -4,11 +4,15 @@ {-# OPTIONS_GHC -Wno-type-defaults #-} module Conversion (benchmark) where import Data.Bits+import Data.Functor.Product import Data.Int import Data.Ratio import Data.Word import Gauge.Benchmark+import Numeric.Floating.IEEE+import qualified Numeric.Floating.IEEE.Internal as IEEE.Internal import Numeric.Rounded.Hardware+import qualified Numeric.Rounded.Hardware.Backend.C as C import Numeric.Rounded.Hardware.Class import Numeric.Rounded.Hardware.Interval @@ -38,62 +42,87 @@ benchmark :: Benchmark benchmark = bgroup "Conversion"- [ let smallInteger = -2^50+2^13+127 :: Integer- mediumInteger = -2^60 + 42 * 2^53 - 137 * 2^24 + 3 :: Integer- largeInteger = -2^100-37*2^80+2^13+127 :: Integer- in bgroup "fromInteger"- [ bench "Double/small" $ nf (fromInteger :: Integer -> Double) smallInteger- , bench "Double/medium" $ nf (fromInteger :: Integer -> Double) mediumInteger- , bench "Double/large" $ nf (fromInteger :: Integer -> Double) largeInteger- , bench "RoundedDouble/ToNearest/small" $ nf (fromInteger :: Integer -> Rounded 'ToNearest Double) smallInteger- , bench "RoundedDouble/ToNearest/medium" $ nf (fromInteger :: Integer -> Rounded 'ToNearest Double) mediumInteger- , bench "RoundedDouble/ToNearest/large" $ nf (fromInteger :: Integer -> Rounded 'ToNearest Double) largeInteger- , bench "RoundedDouble/TowardInf/small" $ nf (fromInteger :: Integer -> Rounded 'TowardInf Double) smallInteger- , bench "RoundedDouble/TowardInf/medium" $ nf (fromInteger :: Integer -> Rounded 'TowardInf Double) mediumInteger- , bench "RoundedDouble/TowardInf/large" $ nf (fromInteger :: Integer -> Rounded 'TowardInf Double) largeInteger- , bench "roundedFromInteger/Double/ToNearest/small" $ nf (roundedFromInteger ToNearest :: Integer -> Double) smallInteger- , bench "roundedFromInteger/Double/ToNearest/medium" $ nf (roundedFromInteger ToNearest :: Integer -> Double) mediumInteger- , bench "roundedFromInteger/Double/ToNearest/large" $ nf (roundedFromInteger ToNearest :: Integer -> Double) largeInteger- , bench "roundedFromInteger/Double/TowardInf/small" $ nf (roundedFromInteger TowardInf :: Integer -> Double) smallInteger- , bench "roundedFromInteger/Double/TowardInf/medium" $ nf (roundedFromInteger TowardInf :: Integer -> Double) mediumInteger- , bench "roundedFromInteger/Double/TowardInf/large" $ nf (roundedFromInteger TowardInf :: Integer -> Double) largeInteger- , bench "IntervalDouble/small" $ nf (fromInteger :: Integer -> Interval Double) smallInteger- , bench "IntervalDouble/medium" $ nf (fromInteger :: Integer -> Interval Double) mediumInteger- , bench "IntervalDouble/large" $ nf (fromInteger :: Integer -> Interval Double) largeInteger- ]- , let smallInteger = -2^50+2^13+127 :: Int64- mediumInteger = -2^60 + 42 * 2^53 - 137 * 2^24 + 3 :: Int64- in bgroup "fromIntegral/Int64"- [ bench "Double/small" $ nf (fromIntegral :: Int64 -> Double) smallInteger- , bench "Double/medium" $ nf (fromIntegral :: Int64 -> Double) mediumInteger- , bench "RoundedDouble/ToNearest/small" $ nf (fromIntegral :: Int64 -> Rounded 'ToNearest Double) smallInteger- , bench "RoundedDouble/ToNearest/medium" $ nf (fromIntegral :: Int64 -> Rounded 'ToNearest Double) mediumInteger- , bench "RoundedDouble/TowardInf/small" $ nf (fromIntegral :: Int64 -> Rounded 'TowardInf Double) smallInteger- , bench "RoundedDouble/TowardInf/medium" $ nf (fromIntegral :: Int64 -> Rounded 'TowardInf Double) mediumInteger- , bench "roundedFromInteger/Double/ToNearest/small" $ nf (roundedFromInteger ToNearest . fromIntegral :: Int64 -> Double) smallInteger- , bench "roundedFromInteger/Double/ToNearest/medium" $ nf (roundedFromInteger ToNearest . fromIntegral :: Int64 -> Double) mediumInteger- , bench "roundedFromInteger/Double/TowardInf/small" $ nf (roundedFromInteger TowardInf . fromIntegral :: Int64 -> Double) smallInteger- , bench "roundedFromInteger/Double/TowardInf/medium" $ nf (roundedFromInteger TowardInf . fromIntegral :: Int64 -> Double) mediumInteger- , bench "int64ToDouble/Double/ToNearest/small" $ nf (int64ToDouble ToNearest :: Int64 -> Double) smallInteger- , bench "int64ToDouble/Double/ToNearest/medium" $ nf (int64ToDouble ToNearest :: Int64 -> Double) mediumInteger- , bench "int64ToDouble/Double/TowardInf/small" $ nf (int64ToDouble TowardInf :: Int64 -> Double) smallInteger- , bench "int64ToDouble/Double/TowardInf/medium" $ nf (int64ToDouble TowardInf :: Int64 -> Double) mediumInteger- ]- , let pi' = 3.14159265358979323846264338327950 :: Rational- smallRational = 22 % 7 :: Rational- largeRational = 78326489123342523452342137498719847192 % 348912374981749170413424213275017 :: Rational- in bgroup "fromRational"- [ bench "Double/decimal" $ nf (fromRational :: Rational -> Double) pi'- , bench "Double/small" $ nf (fromRational :: Rational -> Double) smallRational- , bench "Double/large" $ nf (fromRational :: Rational -> Double) largeRational- , bench "RoundedDouble/ToNearest/decimal" $ nf (fromRational :: Rational -> Rounded 'ToNearest Double) pi'- , bench "RoundedDouble/ToNearest/small" $ nf (fromRational :: Rational -> Rounded 'ToNearest Double) smallRational- , bench "RoundedDouble/ToNearest/large" $ nf (fromRational :: Rational -> Rounded 'ToNearest Double) largeRational- , bench "RoundedDouble/TowardInf/decimal" $ nf (fromRational :: Rational -> Rounded 'TowardInf Double) pi'- , bench "RoundedDouble/TowardInf/small" $ nf (fromRational :: Rational -> Rounded 'TowardInf Double) smallRational- , bench "RoundedDouble/TowardInf/large" $ nf (fromRational :: Rational -> Rounded 'TowardInf Double) largeRational- , bench "IntervalDouble/decimal" $ nf (fromRational :: Rational -> Interval Double) pi'- , bench "IntervalDouble/small" $ nf (fromRational :: Rational -> Interval Double) smallRational- , bench "IntervalDouble/large" $ nf (fromRational :: Rational -> Interval Double) largeRational- ]+ [ bgroup "fromInteger/to Double"+ [ bgroup name $ map ($ value)+ [ bench "plain" . nf (fromInteger :: Integer -> Double)+ , bench "Rounded/ToNearest" . nf (fromInteger :: Integer -> Rounded 'ToNearest Double)+ , bench "Rounded/TowardInf" . nf (fromInteger :: Integer -> Rounded 'TowardInf Double)+ , bench "roundedFromInteger/ToNearest" . nf (roundedFromInteger ToNearest :: Integer -> Double)+ , bench "roundedFromInteger/TowardInf" . nf (roundedFromInteger TowardInf :: Integer -> Double)+ , bench "fp-ieee/ToNearest" . nf (fromIntegerTiesToEven :: Integer -> Double)+ , bench "fp-ieee/TowardInf" . nf (fromIntegerTowardPositive :: Integer -> Double)+ , bench "Interval/default" . nf (fromInteger :: Integer -> Interval Double)+ , bench "Interval/individual" . nf (\n -> (fromIntegerTowardNegative n, fromIntegerTowardPositive n) :: (Double, Double))+ , bench "Interval/fromIntegerR" . nf (\n -> case IEEE.Internal.fromIntegerR n of+ Pair (IEEE.Internal.RoundTowardNegative x) (IEEE.Internal.RoundTowardPositive y) -> (x, y) :: (Double, Double)+ )+ ]+ | (name, value) <- [ ("small", -2^50 + 2^13 + 127)+ , ("medium", -2^60 + 42 * 2^53 - 137 * 2^24 + 3)+ , ("large", -2^100 - 37 * 2^80 + 2^13 + 127)+ ] :: [(String, Integer)] ]+ , bgroup "fromIntegral/Int64->Double"+ [ bgroup name $ map ($ value)+ [ bench "plain" . nf (fromIntegral :: Int64 -> Double)+ , bench "Rounded/ToNearest" . nf (fromIntegral :: Int64 -> Rounded 'ToNearest Double)+ , bench "Rounded/TowardInf" . nf (fromIntegral :: Int64 -> Rounded 'TowardInf Double)+ , bench "roundedFromInteger/ToNearest" . nf (roundedFromInteger ToNearest . fromIntegral :: Int64 -> Double)+ , bench "roundedFromInteger/TowardInf" . nf (roundedFromInteger TowardInf . fromIntegral :: Int64 -> Double)+ , bench "fp-ieee/ToNearest" . nf (fromIntegralTiesToEven :: Int64 -> Double)+ , bench "fp-ieee/TowardInf" . nf (fromIntegralTowardPositive :: Int64 -> Double)+ , bench "int64ToDouble/ToNearest" . nf (int64ToDouble ToNearest :: Int64 -> Double)+ , bench "int64ToDouble/TowardInf" . nf (int64ToDouble TowardInf :: Int64 -> Double)+ , bench "Interval/default" . nf (fromIntegral :: Int64 -> Interval Double)+ , bench "Interval/individual" . nf (\n -> (fromIntegralTowardNegative n, fromIntegralTowardPositive n) :: (Double, Double))+ , bench "Interval/fromIntegralR" . nf (\n -> case IEEE.Internal.fromIntegralR n of+ Pair (IEEE.Internal.RoundTowardNegative x) (IEEE.Internal.RoundTowardPositive y) -> (x, y) :: (Double, Double)+ )+ , bench "Interval/individual/C" . nf (\n -> (C.roundedDoubleFromInt64 TowardNegInf n, C.roundedDoubleFromInt64 TowardInf n))+ ]+ | (name, value) <- [ ("small", -2^50 + 2^13 + 127)+ , ("medium", -2^60 + 42 * 2^53 - 137 * 2^24 + 3)+ ] :: [(String, Int64)]+ ]+ , bgroup "fromIntegral/Word64->Double"+ [ bgroup name $ map ($ value)+ [ bench "plain" . nf (fromIntegral :: Word64 -> Double)+ , bench "Rounded/ToNearest" . nf (fromIntegral :: Word64 -> Rounded 'ToNearest Double)+ , bench "Rounded/TowardInf" . nf (fromIntegral :: Word64 -> Rounded 'TowardInf Double)+ , bench "roundedFromInteger/ToNearest" . nf (roundedFromInteger ToNearest . fromIntegral :: Word64 -> Double)+ , bench "roundedFromInteger/TowardInf" . nf (roundedFromInteger TowardInf . fromIntegral :: Word64 -> Double)+ , bench "fp-ieee/ToNearest" . nf (fromIntegralTiesToEven :: Word64 -> Double)+ , bench "fp-ieee/TowardInf" . nf (fromIntegralTowardPositive :: Word64 -> Double)+ , bench "word64ToDouble/ToNearest" . nf (word64ToDouble ToNearest :: Word64 -> Double)+ , bench "word64ToDouble/TowardInf" . nf (word64ToDouble TowardInf :: Word64 -> Double)+ , bench "Interval/default" . nf (fromIntegral :: Word64 -> Interval Double)+ , bench "Interval/individual" . nf (\n -> (fromIntegralTowardNegative n, fromIntegralTowardPositive n) :: (Double, Double))+ , bench "Interval/fromIntegralR" . nf (\n -> case IEEE.Internal.fromIntegralR n of+ Pair (IEEE.Internal.RoundTowardNegative x) (IEEE.Internal.RoundTowardPositive y) -> (x, y) :: (Double, Double)+ )+ , bench "Interval/individual/C" . nf (\n -> (C.roundedDoubleFromWord64 TowardNegInf n, C.roundedDoubleFromWord64 TowardInf n))+ ]+ | (name, value) <- [ ("small", 2^50 + 2^13 + 127)+ , ("medium", 2^63 + 42 * 2^53 - 137 * 2^24 + 3)+ ] :: [(String, Word64)]+ ]+ , bgroup "fromRational/to Double"+ [ bgroup name $ map ($ value)+ [ bench "plain" . nf (fromRational :: Rational -> Double)+ , bench "Rounded/ToNearest" . nf (fromRational :: Rational -> Rounded 'ToNearest Double)+ , bench "Rounded/TowardInf" . nf (fromRational :: Rational -> Rounded 'TowardInf Double)+ , bench "fp-ieee/ToNearest" . nf (fromRationalTiesToEven :: Rational -> Double)+ , bench "fp-ieee/TowardInf" . nf (fromRationalTowardPositive :: Rational -> Double)+ , bench "Interval/default" . nf (fromRational :: Rational -> Interval Double)+ , bench "Interval/individual" . nf (\x -> (fromRationalTowardNegative x :: Double, fromRationalTowardPositive x :: Double))+ , bench "Interval/fromRationalR" . nf (\x -> case IEEE.Internal.fromRationalR x of+ Pair (IEEE.Internal.RoundTowardNegative a) (IEEE.Internal.RoundTowardPositive b) -> (a, b) :: (Double, Double)+ )+ ]+ | (name, value) <- [ ("decimal", 3.14159265358979323846264338327950)+ , ("binary", 0xcafec0ffeecafec0ffeep-177)+ , ("small", 22 % 7)+ , ("large", 78326489123342523452342137498719847192 % 348912374981749170413424213275017)+ ] :: [(String, Rational)]+ ]+ ]
rounded-hw.cabal view
@@ -4,21 +4,21 @@ -- -- see: https://github.com/sol/hpack ----- hash: 5ac460c3766d27889d7b32a2cbe50446f113a5403dec5125affc436d900f452d+-- hash: c52684ece684d2a3e1a2bde6e2919d961a59bd69724c32a8e7ab1ac0b9230685 name: rounded-hw-version: 0.1.0.0+version: 0.2.0 synopsis: Directed rounding for built-in floating types-description: Please see the README on GitHub at <https://github.com/minoki/rounded-hw#readme>+description: Please see the README on GitHub at <https://github.com/minoki/haskell-floating-point/tree/master/rounded-hw#readme> category: Numeric, Math-homepage: https://github.com/minoki/rounded-hw#readme-bug-reports: https://github.com/minoki/rounded-hw/issues+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-tested-with: GHC == 8.6.5, GHC == 8.8.3+tested-with: GHC == 8.6.5, GHC == 8.8.4, GHC == 8.10.2 build-type: Custom extra-source-files: README.md@@ -29,7 +29,7 @@ source-repository head type: git- location: https://github.com/minoki/rounded-hw+ location: https://github.com/minoki/haskell-floating-point custom-setup setup-depends:@@ -85,7 +85,6 @@ Numeric.Rounded.Hardware.Internal.Constants Numeric.Rounded.Hardware.Internal.Conversion Numeric.Rounded.Hardware.Internal.FloatUtil- Numeric.Rounded.Hardware.Internal.RoundedResult Numeric.Rounded.Hardware.Internal.Show Numeric.Rounded.Hardware.Backend.Default Numeric.Rounded.Hardware.Interval.ElementaryFunctions@@ -95,7 +94,7 @@ array , base >=4.12 && <5 , deepseq- , integer-logarithms+ , fp-ieee ==0.1.* , primitive , tagged , vector@@ -150,7 +149,7 @@ , base >=4.12 && <5 , deepseq , doctest >=0.8- , integer-logarithms+ , fp-ieee ==0.1.* , primitive , vector default-language: Haskell2010@@ -160,7 +159,6 @@ main-is: Spec.hs other-modules: ConstantsSpec- FloatUtilSpec FromIntegerSpec FromRationalSpec IntervalArithmeticSpec@@ -177,8 +175,8 @@ , array , base >=4.12 && <5 , deepseq+ , fp-ieee ==0.1.* , hspec- , integer-logarithms , primitive , random , rounded-hw@@ -210,8 +208,8 @@ array , base >=4.12 && <5 , deepseq+ , fp-ieee ==0.1.* , gauge- , integer-logarithms , primitive , rounded-hw , vector
src/Numeric/Rounded/Hardware/Backend/C.hs view
@@ -32,6 +32,10 @@ , CDouble(..) , VUM.MVector(..) , VU.Vector(..)+ , roundedFloatFromInt64+ , roundedFloatFromWord64+ , roundedDoubleFromInt64+ , roundedDoubleFromWord64 ) where import Control.DeepSeq (NFData (..)) import Data.Bifunctor@@ -39,7 +43,6 @@ import Data.Int (Int64) import Data.Primitive (Prim) import Data.Primitive.ByteArray-import Data.Ratio import Data.Tagged import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM@@ -87,23 +90,6 @@ (F.roundedFromWord64 r x) {-# INLINE roundedFloatFromWord64 #-} -roundedFloatFromInteger :: RoundingMode -> Integer -> Float-roundedFloatFromInteger r x- | -0x1000000 <= x && x <= 0x1000000 {- abs x <= 2^24 -} = fromInteger x- | otherwise = fromInt r x-{-# NOINLINE [1] roundedFloatFromInteger #-}--{-# RULES-"roundeFloatFromInteger/Int" forall r (x :: Int).- roundedFloatFromInteger r (fromIntegral x) = roundedFloatFromInt64 r (fromIntegral x)-"roundeFloatFromInteger/Int64" forall r (x :: Int64).- roundedFloatFromInteger r (fromIntegral x) = roundedFloatFromInt64 r x-"roundeFloatFromInteger/Word" forall r (x :: Word).- roundedFloatFromInteger r (fromIntegral x) = roundedFloatFromWord64 r (fromIntegral x)-"roundeFloatFromInteger/Word64" forall r (x :: Word64).- roundedFloatFromInteger r (fromIntegral x) = roundedFloatFromWord64 r x- #-}- intervalFloatFromInteger :: Integer -> (Rounded 'TowardNegInf Float, Rounded 'TowardInf Float) intervalFloatFromInteger x | -0x1000000 <= x && x <= 0x1000000 {- abs x <= 2^24 -} = (Rounded (fromInteger x), Rounded (fromInteger x))@@ -127,8 +113,8 @@ roundedFusedMultiplyAdd = coerce F.roundedFMA intervalMul x x' y y' = (coerce F.intervalMul_down x x' y y', coerce F.intervalMul_up x x' y y') intervalMulAdd x x' y y' z z' = (coerce F.intervalMulAdd_down x x' y y' z, coerce F.intervalMulAdd_up x x' y y' z')- roundedFromInteger r x = CFloat (roundedFloatFromInteger r x)- intervalFromInteger = coerce intervalFloatFromInteger+ roundedFromInteger r x = CFloat (roundedFromInteger_default r x)+ intervalFromInteger = (coerce `asTypeOf` (bimap (CFloat <$>) (CFloat <$>) .)) intervalFromInteger_default backendNameT = Tagged cBackendName {-# INLINE roundedAdd #-} {-# INLINE roundedSub #-}@@ -142,7 +128,7 @@ roundedDiv = coerce F.roundedDiv intervalDiv x x' y y' = (coerce F.intervalDiv_down x x' y y', coerce F.intervalDiv_up x x' y y') intervalDivAdd x x' y y' z z' = (coerce F.intervalDivAdd_down x x' y y' z, coerce F.intervalDivAdd_up x x' y y' z')- roundedFromRational r x = CFloat $ fromRatio r (numerator x) (denominator x)+ roundedFromRational r x = CFloat (roundedFromRational_default r x) intervalFromRational = (coerce `asTypeOf` (bimap (CFloat <$>) (CFloat <$>) .)) intervalFromRational_default roundedFromRealFloat r x = coerce (roundedFloatFromRealFloat r x) {-# INLINE roundedDiv #-}@@ -209,28 +195,6 @@ (D.roundedFromWord64 r x) {-# INLINE roundedDoubleFromWord64 #-} -roundedDoubleFromInteger :: RoundingMode -> Integer -> Double-roundedDoubleFromInteger r x- | -0x20000000000000 <= x && x <= 0x20000000000000 {- abs x <= 2^53 -} = fromInteger x- | otherwise = fromInt r x-{-# NOINLINE [1] roundedDoubleFromInteger #-}--{-# RULES-"roundedDoubleFromInteger/Int" forall r (x :: Int).- roundedDoubleFromInteger r (fromIntegral x) = roundedDoubleFromInt64 r (fromIntegral x)-"roundedDoubleFromInteger/Int64" forall r (x :: Int64).- roundedDoubleFromInteger r (fromIntegral x) = roundedDoubleFromInt64 r x-"roundedDoubleFromInteger/Word" forall r (x :: Word).- roundedDoubleFromInteger r (fromIntegral x) = roundedDoubleFromWord64 r (fromIntegral x)-"roundedDoubleFromInteger/Word64" forall r (x :: Word64).- roundedDoubleFromInteger r (fromIntegral x) = roundedDoubleFromWord64 r x- #-}--intervalDoubleFromInteger :: Integer -> (Rounded 'TowardNegInf Double, Rounded 'TowardInf Double)-intervalDoubleFromInteger x- | -0x20000000000000 <= x && x <= 0x20000000000000 {- abs x <= 2^53 -} = (Rounded (fromInteger x), Rounded (fromInteger x))- | otherwise = intervalFromInteger_default x- roundedDoubleFromRealFloat :: RealFloat a => RoundingMode -> a -> Double roundedDoubleFromRealFloat r x | isNaN x = 0/0 | isInfinite x = if x > 0 then 1/0 else -1/0@@ -251,8 +215,8 @@ roundedFusedMultiplyAdd = coerce D.roundedFMA intervalMul x x' y y' = (coerce D.intervalMul_down x x' y y', coerce D.intervalMul_up x x' y y') intervalMulAdd x x' y y' z z' = (coerce D.intervalMulAdd_down x x' y y' z, coerce D.intervalMulAdd_up x x' y y' z')- roundedFromInteger = coerce roundedDoubleFromInteger- intervalFromInteger = coerce intervalDoubleFromInteger+ roundedFromInteger r x = CDouble (roundedFromInteger_default r x)+ intervalFromInteger = (coerce `asTypeOf` (bimap (CDouble <$>) (CDouble <$>) .)) intervalFromInteger_default backendNameT = Tagged cBackendName {-# INLINE roundedAdd #-} {-# INLINE roundedSub #-}@@ -266,7 +230,7 @@ roundedDiv = coerce D.roundedDiv intervalDiv x x' y y' = (coerce D.intervalDiv_down x x' y y', coerce D.intervalDiv_up x x' y y') intervalDivAdd x x' y y' z z' = (coerce D.intervalDivAdd_down x x' y y' z, coerce D.intervalDivAdd_up x x' y y' z')- roundedFromRational r x = CDouble $ fromRatio r (numerator x) (denominator x)+ roundedFromRational r x = CDouble (roundedFromRational_default r x) intervalFromRational = (coerce `asTypeOf` (bimap (CDouble <$>) (CDouble <$>) .)) intervalFromRational_default -- TODO: Specialize small case in ***FromRational? roundedFromRealFloat r x = coerce (roundedDoubleFromRealFloat r x)
src/Numeric/Rounded/Hardware/Backend/Default.hs view
@@ -1,13 +1,13 @@ {-# LANGUAGE CPP #-}-{-# LANGUAGE StandaloneDeriving #-}-{-# LANGUAGE DerivingVia #-} {-# LANGUAGE DataKinds #-}+{-# LANGUAGE DerivingVia #-} {-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE StandaloneDeriving #-} {-# OPTIONS_GHC -Wno-orphans -Wno-unused-imports #-} module Numeric.Rounded.Hardware.Backend.Default () where-import Numeric.Rounded.Hardware.Internal.Class import qualified Numeric.Rounded.Hardware.Backend.ViaRational as VR+import Numeric.Rounded.Hardware.Internal.Class #ifdef USE_FFI import qualified Numeric.Rounded.Hardware.Backend.C as C #ifdef USE_GHC_PRIM@@ -20,10 +20,11 @@ import Numeric.Rounded.Hardware.Backend.Float128 () #endif #endif+import Data.Coerce import qualified Data.Vector.Storable as VS import qualified Data.Vector.Unboxed as VU+import Numeric.Floating.IEEE import Unsafe.Coerce-import Data.Coerce #ifdef USE_FFI #ifdef USE_GHC_PRIM@@ -91,19 +92,19 @@ -- orphaned rules {-# RULES "fromIntegral/a->Rounded ToNearest Float"- forall x. fromIntegral x = Rounded (roundedFromInteger ToNearest (fromIntegral x)) :: Rounded 'ToNearest Float+ fromIntegral = \x -> (Rounded (fromIntegralTiesToEven x) :: Rounded 'ToNearest Float) "fromIntegral/a->Rounded TowardInf Float"- forall x. fromIntegral x = Rounded (roundedFromInteger TowardInf (fromIntegral x)) :: Rounded 'TowardInf Float+ fromIntegral = \x -> (Rounded (fromIntegralTowardPositive x) :: Rounded 'TowardInf Float) "fromIntegral/a->Rounded TowardNegInf Float"- forall x. fromIntegral x = Rounded (roundedFromInteger TowardNegInf (fromIntegral x)) :: Rounded 'TowardNegInf Float+ fromIntegral = \x -> (Rounded (fromIntegralTowardNegative x) :: Rounded 'TowardNegInf Float) "fromIntegral/a->Rounded TowardZero Float"- forall x. fromIntegral x = Rounded (roundedFromInteger TowardZero (fromIntegral x)) :: Rounded 'TowardZero Float+ fromIntegral = \x -> (Rounded (fromIntegralTowardZero x) :: Rounded 'TowardZero Float) "fromIntegral/a->Rounded ToNearest Double"- forall x. fromIntegral x = Rounded (roundedFromInteger ToNearest (fromIntegral x)) :: Rounded 'ToNearest Double+ fromIntegral = \x -> (Rounded (fromIntegralTiesToEven x) :: Rounded 'ToNearest Double) "fromIntegral/a->Rounded TowardInf Double"- forall x. fromIntegral x = Rounded (roundedFromInteger TowardInf (fromIntegral x)) :: Rounded 'TowardInf Double+ fromIntegral = \x -> (Rounded (fromIntegralTowardPositive x) :: Rounded 'TowardInf Double) "fromIntegral/a->Rounded TowardNegInf Double"- forall x. fromIntegral x = Rounded (roundedFromInteger TowardNegInf (fromIntegral x)) :: Rounded 'TowardNegInf Double+ fromIntegral = \x -> (Rounded (fromIntegralTowardNegative x) :: Rounded 'TowardNegInf Double) "fromIntegral/a->Rounded TowardZero Double"- forall x. fromIntegral x = Rounded (roundedFromInteger TowardZero (fromIntegral x)) :: Rounded 'TowardZero Double+ fromIntegral = \x -> (Rounded (fromIntegralTowardZero x) :: Rounded 'TowardZero Double) #-}
src/Numeric/Rounded/Hardware/Backend/FastFFI.hs view
@@ -27,6 +27,7 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UnboxedTuples #-} {-# LANGUAGE UnliftedFFITypes #-}+{-# OPTIONS_GHC -fobject-code #-} module Numeric.Rounded.Hardware.Backend.FastFFI ( CDouble(..) , fastIntervalAdd
src/Numeric/Rounded/Hardware/Backend/Float128.hs view
@@ -4,7 +4,6 @@ module Numeric.Rounded.Hardware.Backend.Float128 ( ) where-import Data.Ratio import Data.Tagged import Foreign.C.String (CString, peekCString) import Foreign.Marshal (alloca, with)@@ -125,7 +124,7 @@ roundedSub = roundedSub_f128 roundedMul = roundedMul_f128 roundedFusedMultiplyAdd = roundedFMA_f128- roundedFromInteger = fromInt+ roundedFromInteger = roundedFromInteger_default intervalFromInteger = intervalFromInteger_default backendNameT = Tagged cBackendName {-# INLINE roundedAdd #-}@@ -137,7 +136,7 @@ instance RoundedFractional Float128 where roundedDiv = roundedDiv_f128- roundedFromRational r x = fromRatio r (numerator x) (denominator x)+ roundedFromRational = roundedFromRational_default intervalFromRational = intervalFromRational_default {-# INLINE roundedDiv #-} {-# INLINE roundedFromRational #-}
src/Numeric/Rounded/Hardware/Backend/ViaRational.hs view
@@ -8,8 +8,6 @@ import Control.DeepSeq (NFData (..)) import Control.Exception (assert) import Data.Coerce-import Data.Functor.Product-import Data.Ratio import Data.Tagged import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM@@ -21,7 +19,7 @@ import Numeric.Rounded.Hardware.Internal.Class import Numeric.Rounded.Hardware.Internal.Constants import Numeric.Rounded.Hardware.Internal.Conversion-import Numeric.Rounded.Hardware.Internal.FloatUtil (nextDown, nextUp)+import Numeric.Floating.IEEE (isFinite, nextDown, nextUp) newtype ViaRational a = ViaRational a deriving (Eq,Ord,Show,Generic,Num,Storable)@@ -59,20 +57,21 @@ | isNaN x || isNaN y || isInfinite x || isInfinite y || isNegativeZero x || isNegativeZero y = ViaRational (x * y) | otherwise = roundedFromRational r (toRational x * toRational y) roundedFusedMultiplyAdd r (ViaRational x) (ViaRational y) (ViaRational z)- | isNaN x || isNaN y || isNaN z || isInfinite x || isInfinite y || isInfinite z = ViaRational (x * y + z)- | otherwise = case toRational x * toRational y + toRational z of+ | isFinite x && isFinite y && isFinite z = case toRational x * toRational y + toRational z of 0 -> if z == 0 && isNegativeZero (x * y) == isNegativeZero z then ViaRational z else ViaRational roundedZero w -> roundedFromRational r w+ | isFinite x && isFinite y = ViaRational z -- Infinity or NaN+ | otherwise = ViaRational (x * y + z) where roundedZero = case r of ToNearest -> 0 TowardNegInf -> -0 TowardInf -> 0 TowardZero -> 0- roundedFromInteger r x = ViaRational (fromInt r x)- intervalFromInteger x = case fromIntF x :: Product (Rounded 'TowardNegInf) (Rounded 'TowardInf) a of- Pair a b -> (ViaRational <$> a, ViaRational <$> b)+ roundedFromInteger r x = ViaRational (roundedFromInteger_default r x)+ intervalFromInteger x = case intervalFromInteger_default x of+ (a, b) -> (ViaRational <$> a, ViaRational <$> b) backendNameT = Tagged "via Rational" {-# INLINE roundedFromInteger #-} {-# INLINE intervalFromInteger #-}@@ -83,13 +82,13 @@ roundedDiv r (ViaRational x) (ViaRational y) | isNaN x || isNaN y || isInfinite x || isInfinite y || x == 0 || y == 0 = ViaRational (x / y) | otherwise = roundedFromRational r (toRational x / toRational y)- roundedFromRational r x = ViaRational $ fromRatio r (numerator x) (denominator x)+ roundedFromRational r x = ViaRational $ roundedFromRational_default r x roundedFromRealFloat r x | isNaN x = ViaRational (0/0) | isInfinite x = ViaRational (if x > 0 then 1/0 else -1/0) | isNegativeZero x = ViaRational (-0) | otherwise = roundedFromRational r (toRational x)- intervalFromRational x = case fromRatioF (numerator x) (denominator x) :: Product (Rounded 'TowardNegInf) (Rounded 'TowardInf) a of- Pair a b -> (ViaRational <$> a, ViaRational <$> b)+ intervalFromRational x = case intervalFromRational_default x of+ (a, b) -> (ViaRational <$> a, ViaRational <$> b) {-# INLINE roundedFromRational #-} {-# INLINE intervalFromRational #-} {-# SPECIALIZE instance RoundedFractional (ViaRational Float) #-}
src/Numeric/Rounded/Hardware/Backend/X87LongDouble.hs view
@@ -3,7 +3,6 @@ module Numeric.Rounded.Hardware.Backend.X87LongDouble ( ) where-import Data.Ratio import Data.Tagged import Foreign.C.String (CString, peekCString) import Foreign.Marshal (alloca, with)@@ -124,7 +123,7 @@ roundedSub = roundedSub_ld roundedMul = roundedMul_ld roundedFusedMultiplyAdd = roundedFMA_ld- roundedFromInteger = fromInt+ roundedFromInteger = roundedFromInteger_default intervalFromInteger = intervalFromInteger_default backendNameT = Tagged cBackendName {-# INLINE roundedAdd #-}@@ -138,7 +137,7 @@ -- Note that 'LongDouble' may not work correctly on Win64. instance RoundedFractional LongDouble where roundedDiv = roundedDiv_ld- roundedFromRational r x = fromRatio r (numerator x) (denominator x)+ roundedFromRational = roundedFromRational_default intervalFromRational = intervalFromRational_default {-# INLINE roundedDiv #-} {-# INLINE roundedFromRational #-}
src/Numeric/Rounded/Hardware/Internal.hs view
@@ -8,6 +8,5 @@ import Numeric.Rounded.Hardware.Internal.Constants as Internal import Numeric.Rounded.Hardware.Internal.Conversion as Internal import Numeric.Rounded.Hardware.Internal.FloatUtil as Internal-import Numeric.Rounded.Hardware.Internal.RoundedResult as Internal import Numeric.Rounded.Hardware.Internal.Rounding as Internal import Numeric.Rounded.Hardware.Internal.Show as Internal
src/Numeric/Rounded/Hardware/Internal/Class.hs view
@@ -18,6 +18,8 @@ import Data.Ratio import Data.Tagged import qualified Data.Vector.Generic as VG+import Numeric.Floating.IEEE+import Numeric.Rounded.Hardware.Internal.Conversion import Numeric.Rounded.Hardware.Internal.Rounding import Prelude hiding (fromInteger, fromRational, recip, sqrt, (*), (+), (-), (/))@@ -30,6 +32,8 @@ roundedMul :: RoundingMode -> a -> a -> a roundedFusedMultiplyAdd :: RoundingMode -> a -> a -> a -> a roundedFromInteger :: RoundingMode -> Integer -> a+ default roundedFromInteger :: RealFloat a => RoundingMode -> Integer -> a+ roundedFromInteger = roundedFromInteger_default -- roundedToFloat :: RoundingMode -> a -> Float -- roundedToDouble :: RoundingMode -> a -> Double @@ -95,6 +99,8 @@ default roundedRecip :: Num a => RoundingMode -> a -> a roundedRecip r = roundedDiv r 1 roundedFromRational :: RoundingMode -> Rational -> a+ default roundedFromRational :: RealFloat a => RoundingMode -> Rational -> a+ roundedFromRational = roundedFromRational_default roundedFromRealFloat :: RealFloat b => RoundingMode -> b -> a default roundedFromRealFloat :: (Fractional a, RealFloat b) => RoundingMode -> b -> a roundedFromRealFloat r x | isNaN x = 0 Prelude./ 0
src/Numeric/Rounded/Hardware/Internal/Conversion.hs view
@@ -1,203 +1,44 @@-{-# LANGUAGE HexFloatLiterals #-}-{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-}-{-# LANGUAGE ScopedTypeVariables #-} module Numeric.Rounded.Hardware.Internal.Conversion- ( fromInt- , fromIntF+ ( roundedFromInteger_default+ , roundedFromRational_default , intervalFromInteger_default- , fromRatio- , fromRatioF+ , intervalFromIntegral , intervalFromRational_default ) where-import Numeric.Rounded.Hardware.Internal.Rounding-import Numeric.Rounded.Hardware.Internal.RoundedResult-import Numeric.Rounded.Hardware.Internal.FloatUtil-import Data.Bits-import Data.Functor.Product-import Math.NumberTheory.Logarithms (integerLog2')-import Data.Ratio-import Control.Exception (assert)--- import GHC.Integer.Logarithms.Internals (integerLog2IsPowerOf2#)--- integerLog2IsPowerOf2# :: Integer -> (# Int#, Int# #)--intervalFromInteger_default :: RealFloat a => Integer -> (Rounded 'TowardNegInf a, Rounded 'TowardInf a)-intervalFromInteger_default x = case fromIntF x of Pair a b -> (a, b)-{-# SPECIALIZE intervalFromInteger_default :: Integer -> (Rounded 'TowardNegInf Float, Rounded 'TowardInf Float) #-}-{-# SPECIALIZE intervalFromInteger_default :: Integer -> (Rounded 'TowardNegInf Double, Rounded 'TowardInf Double) #-}--intervalFromRational_default :: RealFloat a => Rational -> (Rounded 'TowardNegInf a, Rounded 'TowardInf a)-intervalFromRational_default x = case fromRatioF (numerator x) (denominator x) of Pair a b -> (a, b)-{-# SPECIALIZE intervalFromRational_default :: Rational -> (Rounded 'TowardNegInf Float, Rounded 'TowardInf Float) #-}-{-# SPECIALIZE intervalFromRational_default :: Rational -> (Rounded 'TowardNegInf Double, Rounded 'TowardInf Double) #-}--fromInt :: RealFloat a => RoundingMode -> Integer -> a-fromInt r n = withRoundingMode (fromIntF n) r-{-# SPECIALIZE fromInt :: RoundingMode -> Integer -> Float #-}-{-# SPECIALIZE fromInt :: RoundingMode -> Integer -> Double #-}+import Data.Functor.Product+import Numeric.Floating.IEEE+import Numeric.Floating.IEEE.Internal (fromIntegerR, fromIntegralR,+ fromRationalR,+ roundTowardNegative,+ roundTowardPositive)+import Numeric.Rounded.Hardware.Internal.Rounding -fromIntF :: forall a f. (RealFloat a, Result f) => Integer -> f a-fromIntF !_ | floatRadix (undefined :: a) /= 2 = error "radix other than 2 is not supported"-fromIntF 0 = exact 0-fromIntF n | n < 0 = negate <$> withOppositeRoundingMode (fromPositiveIntF (- n))- | otherwise = fromPositiveIntF n-{-# INLINE fromIntF #-}+roundedFromInteger_default :: RealFloat a => RoundingMode -> Integer -> a+roundedFromInteger_default ToNearest = fromIntegerTiesToEven+roundedFromInteger_default TowardZero = fromIntegerTowardZero+roundedFromInteger_default TowardInf = fromIntegerTowardPositive+roundedFromInteger_default TowardNegInf = fromIntegerTowardNegative+{-# INLINE roundedFromInteger_default #-} --- n > 0-fromPositiveIntF :: forall a f. (RealFloat a, Result f) => Integer -> f a-fromPositiveIntF !n- = let !k = integerLog2' n -- floor (log2 n)- -- 2^k <= n < 2^(k+1)- !fDigits = floatDigits (undefined :: a) -- 53 for Double- in if k < fDigits- then exact (fromInteger n)- else let e = k - (fDigits - 1)- -- (!q, !r) = n `quotRem` (1 `unsafeShiftL` e)- q = n `unsafeShiftR` e- r = n .&. ((1 `unsafeShiftL` e) - 1)- -- 2^52 <= q < 2^53, 0 <= r < 2^(k-52)- (_expMin, !expMax) = floatRange (undefined :: a) -- (-1021, 1024) for Double- in if k >= expMax- then- -- infinity- inexact (1 / 0) -- ToNearest- (1 / 0) -- TowardInf- maxFinite_ieee -- TowardNegInf- maxFinite_ieee -- TowardZero- else- if r == 0- then exact $ encodeFloat q e -- exact- else- -- inexact- let down = encodeFloat q e- up = encodeFloat (q + 1) e- toNearest = case compare r (1 `unsafeShiftL` (e-1)) of- LT -> down- EQ | even q -> down- | otherwise -> up- GT -> up- in inexact toNearest up down down-{-# SPECIALIZE fromPositiveIntF :: Integer -> DynamicRoundingMode Float #-}-{-# SPECIALIZE fromPositiveIntF :: Integer -> OppositeRoundingMode DynamicRoundingMode Float #-}-{-# SPECIALIZE fromPositiveIntF :: Rounding r => Integer -> Rounded r Float #-}-{-# SPECIALIZE fromPositiveIntF :: Rounding r => Integer -> OppositeRoundingMode (Rounded r) Float #-}-{-# SPECIALIZE fromPositiveIntF :: Integer -> Product (Rounded 'TowardNegInf) (Rounded 'TowardInf) Float #-}-{-# SPECIALIZE fromPositiveIntF :: Integer -> OppositeRoundingMode (Product (Rounded 'TowardNegInf) (Rounded 'TowardInf)) Float #-}-{-# SPECIALIZE fromPositiveIntF :: Integer -> DynamicRoundingMode Double #-}-{-# SPECIALIZE fromPositiveIntF :: Integer -> OppositeRoundingMode DynamicRoundingMode Double #-}-{-# SPECIALIZE fromPositiveIntF :: Rounding r => Integer -> Rounded r Double #-}-{-# SPECIALIZE fromPositiveIntF :: Rounding r => Integer -> OppositeRoundingMode (Rounded r) Double #-}-{-# SPECIALIZE fromPositiveIntF :: Integer -> Product (Rounded 'TowardNegInf) (Rounded 'TowardInf) Double #-}-{-# SPECIALIZE fromPositiveIntF :: Integer -> OppositeRoundingMode (Product (Rounded 'TowardNegInf) (Rounded 'TowardInf)) Double #-}+roundedFromRational_default :: RealFloat a => RoundingMode -> Rational -> a+roundedFromRational_default ToNearest = fromRationalTiesToEven+roundedFromRational_default TowardZero = fromRationalTowardZero+roundedFromRational_default TowardInf = fromRationalTowardPositive+roundedFromRational_default TowardNegInf = fromRationalTowardNegative+{-# INLINE roundedFromRational_default #-} -fromRatio :: (RealFloat a)- => RoundingMode- -> Integer -- ^ numerator- -> Integer -- ^ denominator- -> a-fromRatio r n d = withRoundingMode (fromRatioF n d) r-{-# SPECIALIZE fromRatio :: RoundingMode -> Integer -> Integer -> Float #-}-{-# SPECIALIZE fromRatio :: RoundingMode -> Integer -> Integer -> Double #-}+intervalFromInteger_default :: RealFloat a => Integer -> (Rounded 'TowardNegInf a, Rounded 'TowardInf a)+intervalFromInteger_default x = case fromIntegerR x of+ Pair a b -> (Rounded (roundTowardNegative a), Rounded (roundTowardPositive b))+{-# INLINE intervalFromInteger_default #-} -fromRatioF :: forall a f. (RealFloat a, Result f)- => Integer -- ^ numerator- -> Integer -- ^ denominator- -> f a-fromRatioF !_ !_ | floatRadix (undefined :: a) /= 2 = error "radix other than 2 is not supported"-fromRatioF 0 _ = exact 0-fromRatioF n 0 | n > 0 = exact (1 / 0) -- positive infinity- | otherwise = exact (- 1 / 0) -- negative infinity-fromRatioF n d | d < 0 = error "fromRatio: negative denominator"- | n < 0 = negate <$> withOppositeRoundingMode (fromPositiveRatioF (- n) d)- | otherwise = fromPositiveRatioF n d-{-# INLINE fromRatioF #-}+intervalFromRational_default :: RealFloat a => Rational -> (Rounded 'TowardNegInf a, Rounded 'TowardInf a)+intervalFromRational_default x = case fromRationalR x of+ Pair a b -> (Rounded (roundTowardNegative a), Rounded (roundTowardPositive b))+{-# INLINE intervalFromRational_default #-} --- n > 0, d > 0-fromPositiveRatioF :: forall a f. (RealFloat a, Result f)- => Integer -> Integer -> f a-fromPositiveRatioF !n !d- = let ln, ld, e :: Int- ln = integerLog2' n- ld = integerLog2' d- e = ln - ld - fDigits- q, r, d_ :: Integer- d_ | e >= 0 = d `unsafeShiftL` e- | otherwise = d- (!q, !r) | e >= 0 = n `quotRem` d_- | otherwise = (n `unsafeShiftL` (-e)) `quotRem` d- -- e >= 0: n = q * (d * 2^e) + r, 0 <= r < d * 2^e- -- e <= 0: n * 2^(-e) = q * d + r, 0 <= r < d- -- n / d * 2^^(-e) = q + r / d_- -- 52 <= log2 q < 54- q', r', d' :: Integer- e' :: Int- (!q', !r', !d', !e') | q < (1 `unsafeShiftL` fDigits) = (q, r, d_, e)- | otherwise = let (q'', r'') = q `quotRem` 2- in (q'', r'' * d_ + r, 2 * d_, e + 1)- -- n / d * 2^^(-e') = q' + r' / d', 2^52 <= q' < 2^53, 0 <= r' < d'- -- q' * 2^^e' <= n/d < (q'+1) * 2^^e', 2^52 <= q' < 2^53- -- (q'/2^53) * 2^^(e'+53) <= n/d < (q'+1)/2^53 * 2^^(e'+53), 1/2 <= q'/2^53 < 1- -- normal: 0x1p-1022 <= x <= 0x1.fffffffffffffp+1023- in assert (n % d * 2^^(-e) == fromInteger q + r % d_) $- assert (n % d * 2^^(-e') == fromInteger q' + r' % d') $- if expMin <= e' + fDigits && e' + fDigits <= expMax- then- -- normal- if r' == 0- then- exact $ encodeFloat q' e' -- exact- else- -- inexact- let down = encodeFloat q' e'- up = encodeFloat (q' + 1) e' -- may be infinity- toNearest = case compare (2 * r') d' of- LT -> down- EQ | even q' -> down- | otherwise -> up -- q' + 1 is even- GT -> up- in inexact toNearest up down down- else- -- infinity or subnormal- if expMax <= e' + fDigits- then- -- infinity- inexact (1 / 0) -- ToNearest- (1 / 0) -- TowardInf- maxFinite_ieee -- TowardNegInf- maxFinite_ieee -- TowardZero- else- -- subnormal- -- e' + fDigits < expMin (or, e' < expMin - fDigits = -1074)- -- 0 <= rounded(n/d) <= 2^(expMin - 1) = 0x1p-1022, minimum (positive) subnormal: 0x1p-1074- let (!q'', !r'') = q' `quotRem` (1 `unsafeShiftL` (expMin - fDigits - e'))- -- q' = q'' * 2^(expMin - fDigits - e') + r'', 0 <= r'' < 2^(expMin - fDigits - e')- -- 2^(fDigits-1) <= q' = q'' * 2^(expMin - fDigits - e') + r'' < 2^fDigits- -- n / d * 2^^(-e') = q' + r' / d' = q'' * 2^(expMin - fDigits - e') + r'' + r' / d'- -- n / d = q'' * 2^^(expMin - fDigits) + (r'' + r' / d') * 2^^e'- -- 0 <= r'' < 2^(expMin - fDigits - e')- in if r' == 0 && r'' == 0- then exact $ encodeFloat q'' (expMin - fDigits) -- exact- else let down = encodeFloat q'' (expMin - fDigits)- up = encodeFloat (q'' + 1) (expMin - fDigits)- toNearest = case compare r'' (1 `unsafeShiftL` (expMin - fDigits - e' - 1)) of- LT -> down- GT -> up- EQ | r' /= 0 -> up- | even q' -> down- | otherwise -> up- in inexact toNearest up down down- where- !fDigits = floatDigits (undefined :: a) -- 53 for Double- (!expMin, !expMax) = floatRange (undefined :: a) -- (-1021, 1024) for Double-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> DynamicRoundingMode Float #-}-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> OppositeRoundingMode DynamicRoundingMode Float #-}-{-# SPECIALIZE fromPositiveRatioF :: Rounding r => Integer -> Integer -> Rounded r Float #-}-{-# SPECIALIZE fromPositiveRatioF :: Rounding r => Integer -> Integer -> OppositeRoundingMode (Rounded r) Float #-}-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> Product (Rounded 'TowardNegInf) (Rounded 'TowardInf) Float #-}-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> OppositeRoundingMode (Product (Rounded 'TowardNegInf) (Rounded 'TowardInf)) Float #-}-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> DynamicRoundingMode Double #-}-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> OppositeRoundingMode DynamicRoundingMode Double #-}-{-# SPECIALIZE fromPositiveRatioF :: Rounding r => Integer -> Integer -> Rounded r Double #-}-{-# SPECIALIZE fromPositiveRatioF :: Rounding r => Integer -> Integer -> OppositeRoundingMode (Rounded r) Double #-}-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> Product (Rounded 'TowardNegInf) (Rounded 'TowardInf) Double #-}-{-# SPECIALIZE fromPositiveRatioF :: Integer -> Integer -> OppositeRoundingMode (Product (Rounded 'TowardNegInf) (Rounded 'TowardInf)) Double #-}+intervalFromIntegral :: (Integral i, RealFloat a) => i -> (Rounded 'TowardNegInf a, Rounded 'TowardInf a)+intervalFromIntegral x = case fromIntegralR x of+ Pair a b -> (Rounded (roundTowardNegative a), Rounded (roundTowardPositive b))+{-# INLINE intervalFromIntegral #-}
src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs view
@@ -4,317 +4,11 @@ ( nextUp , nextDown , nextTowardZero- , minPositive_ieee- , maxFinite_ieee , distanceUlp , fusedMultiplyAdd ) where-import Data.Bits import Data.Ratio-import GHC.Float (castDoubleToWord64, castFloatToWord32,- castWord32ToFloat, castWord64ToDouble)---- $setup--- >>> :set -XHexFloatLiterals -XNumericUnderscores---- |--- prop> (minPositive_ieee :: Double) == 0x1p-1074--- prop> (minPositive_ieee :: Float) == 0x1p-149-minPositive_ieee :: RealFloat a => a-minPositive_ieee = let d = floatDigits x- (expMin,_expMax) = floatRange x- x = encodeFloat 1 (expMin - d)- in x-{-# SPECIALIZE minPositive_ieee :: Double #-}-{-# SPECIALIZE minPositive_ieee :: Float #-}---- |--- prop> (maxFinite_ieee :: Double) == 0x1.ffff_ffff_ffff_fp+1023--- prop> (maxFinite_ieee :: Float) == 0x1.fffffep+127-maxFinite_ieee :: RealFloat a => a-maxFinite_ieee = let d = floatDigits x- (_expMin,expMax) = floatRange x- r = floatRadix x- x = encodeFloat (r ^! d - 1) (expMax - d)- in x-{-# SPECIALIZE maxFinite_ieee :: Double #-}-{-# SPECIALIZE maxFinite_ieee :: Float #-}---- A variant of (^) allowing constant folding for base = 2-infixr 8 ^!-(^!) :: Integer -> Int -> Integer-(^!) = (^)-{-# INLINE [2] (^!) #-}-{-# RULES-"2^!" forall y. 2 ^! y = staticIf (y >= 0) (1 `shiftL` y) (2 ^ y)- #-}--staticIf :: Bool -> a -> a -> a-staticIf _ _ x = x-{-# INLINE [0] staticIf #-}-{-# RULES-"staticIf/True" forall x y. staticIf True x y = x-"staticIf/False" forall x y. staticIf False x y = y- #-}---- |--- prop> nextUp 1 == (0x1.0000_0000_0000_1p0 :: Double)--- prop> nextUp 1 == (0x1.000002p0 :: Float)--- prop> nextUp (1/0) == (1/0 :: Double)--- prop> nextUp (-1/0) == (- maxFinite_ieee :: Double)--- prop> nextUp 0 == (0x1p-1074 :: Double)--- prop> nextUp (-0) == (0x1p-1074 :: Double)--- prop> nextUp (-0x1p-1074) == (-0 :: Double)--- prop> isNegativeZero (nextUp (-0x1p-1074) :: Double)-nextUp :: RealFloat a => a -> a-nextUp x | not (isIEEE x) = error "non-IEEE numbers are not supported"- | floatRadix x /= 2 = error "non-binary types are not supported"- | isNaN x || (isInfinite x && x > 0) = x -- NaN or positive infinity- | x >= 0 = nextUp_ieee_positive x- | otherwise = - nextDown_ieee_positive (- x)-{-# INLINE [1] nextUp #-}---- |--- prop> nextDown 1 == (0x1.ffff_ffff_ffff_fp-1 :: Double)--- prop> nextDown 1 == (0x1.fffffep-1 :: Float)--- prop> nextDown (1/0) == (maxFinite_ieee :: Double)--- prop> nextDown (-1/0) == (-1/0 :: Double)--- prop> nextDown 0 == (-0x1p-1074 :: Double)--- prop> nextDown (-0) == (-0x1p-1074 :: Double)--- prop> nextDown 0x1p-1074 == (0 :: Double)-nextDown :: RealFloat a => a -> a-nextDown x | not (isIEEE x) = error "non-IEEE numbers are not supported"- | floatRadix x /= 2 = error "non-binary types are not supported"- | isNaN x || (isInfinite x && x < 0) = x -- NaN or negative infinity- | x >= 0 = nextDown_ieee_positive x- | otherwise = - nextUp_ieee_positive (- x)-{-# INLINE [1] nextDown #-}---- |--- prop> nextTowardZero 1 == (0x1.ffff_ffff_ffff_fp-1 :: Double)--- prop> nextTowardZero 1 == (0x1.fffffep-1 :: Float)--- prop> nextTowardZero (1/0) == (maxFinite_ieee :: Double)--- prop> nextTowardZero (-1/0) == (-maxFinite_ieee :: Double)--- prop> nextTowardZero 0 == (0 :: Double)--- prop> isNegativeZero (nextTowardZero (-0 :: Double))--- prop> nextTowardZero 0x1p-1074 == (0 :: Double)-nextTowardZero :: RealFloat a => a -> a-nextTowardZero x | not (isIEEE x) = error "non-IEEE numbers are not supported"- | floatRadix x /= 2 = error "non-binary types are not supported "- | isNaN x || x == 0 = x -- NaN or zero- | x >= 0 = nextDown_ieee_positive x- | otherwise = - nextDown_ieee_positive (- x)-{-# INLINE [1] nextTowardZero #-}--nextUp_ieee_positive :: RealFloat a => a -> a-nextUp_ieee_positive x- | isNaN x || x < 0 = error "nextUp_ieee_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 encodeFloat (m + 1) e -- normal- else let m' = m `shiftR` (expMin - d - e)- in encodeFloat (m' + 1) (expMin - d) -- subnormal- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,_expMax) = floatRange x -- (-1021,1024) for Double-{-# INLINE nextUp_ieee_positive #-}--nextDown_ieee_positive :: RealFloat a => a -> a-nextDown_ieee_positive x- | isNaN x || x < 0 = error "nextDown_ieee_positive"- | isInfinite x = encodeFloat ((1 `unsafeShiftL` d) - 1) (expMax - d) -- max finite- | 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 .&. m1 == 0- then encodeFloat (2 * m - 1) (e - 1)- else encodeFloat m1 e- else -- subnormal- let m' = m `shiftR` (expMin - d - e)- in encodeFloat (m' - 1) (expMin - d)- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,expMax) = floatRange x -- (-1021,1024) for Double-{-# INLINE nextDown_ieee_positive #-}--{-# RULES-"nextUp/Float" [~1] nextUp = nextUpFloat-"nextUp/Double" [~1] nextUp = nextUpDouble-"nextDown/Float" [~1] nextDown = nextDownFloat-"nextDown/Double" [~1] nextDown = nextDownDouble-"nextTowardZero/Float" [~1] nextTowardZero = nextTowardZeroFloat-"nextTowardZero/Double" [~1] nextTowardZero = nextTowardZeroDouble- #-}---- |--- prop> nextUpFloat 1 == 0x1.000002p0--- prop> nextUpFloat (1/0) == 1/0--- prop> nextUpFloat (-1/0) == - maxFinite_ieee--- prop> nextUpFloat 0 == 0x1p-149--- prop> nextUpFloat (-0) == 0x1p-149--- prop> isNegativeZero (nextUpFloat (-0x1p-149))-nextUpFloat :: Float -> Float-nextUpFloat x- | not (isIEEE x) || floatRadix x /= 2 || d /= 24 || expMin /= -125 || expMax /= 128 = error "rounded-hw assumes Float is IEEE binary32"- | isNaN x = x -- NaN -> itself- | isNegativeZero x = encodeFloat 1 (expMin - d) -- -0 -> min positive- | x < 0 = castWord32ToFloat (castFloatToWord32 x - 1) -- negative- | otherwise = case castFloatToWord32 x of- 0x7f80_0000 -> x -- positive infinity -> itself- w -> castWord32ToFloat (w + 1) -- positive- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,expMax) = floatRange x -- (-1021,1024) for Double- -- Note: castFloatToWord32 is buggy on GHC <= 8.8 on x86_64, so we can't use it to test for NaN or negative number- -- https://gitlab.haskell.org/ghc/ghc/issues/16617---- |--- prop> nextUpDouble 1 == 0x1.0000_0000_0000_1p0--- prop> nextUpDouble (1/0) == 1/0--- prop> nextUpDouble (-1/0) == - maxFinite_ieee--- prop> nextUpDouble 0 == 0x1p-1074--- prop> nextUpDouble (-0) == 0x1p-1074--- prop> isNegativeZero (nextUpDouble (-0x1p-1074))-nextUpDouble :: Double -> Double-nextUpDouble x- | not (isIEEE x) || floatRadix x /= 2 || d /= 53 || expMin /= -1021 || expMax /= 1024 = error "rounded-hw assumes Double is IEEE binary64"- | otherwise = case castDoubleToWord64 x of- w | w .&. 0x7ff0_0000_0000_0000 == 0x7ff0_0000_0000_0000- , w /= 0xfff0_0000_0000_0000 -> x -- NaN or positive infinity -> itself- 0x8000_0000_0000_0000 -> encodeFloat 1 (expMin - d) -- -0 -> min positive- w | testBit w 63 -> castWord64ToDouble (w - 1) -- negative- | otherwise -> castWord64ToDouble (w + 1) -- positive- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,expMax) = floatRange x -- (-1021,1024) for Double---- |--- prop> nextDownFloat 1 == 0x1.fffffep-1--- prop> nextDownFloat (1/0) == maxFinite_ieee--- 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- | not (isIEEE x) || floatRadix x /= 2 || d /= 24 || expMin /= -125 || expMax /= 128 = error "rounded-hw assumes Float is IEEE binary32"- | isNaN x || (isInfinite x && x < 0) = x -- NaN or negative infinity -> itself- | isNegativeZero x || x < 0 = castWord32ToFloat (castFloatToWord32 x + 1) -- negative- | x == 0 = encodeFloat (-1) (expMin - d) -- +0 -> max negative- | otherwise = castWord32ToFloat (castFloatToWord32 x - 1) -- positive- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,expMax) = floatRange x -- (-1021,1024) for Double- -- Note: castFloatToWord32 is buggy on GHC <= 8.8 on x86_64, so we can't use it to test for NaN or negative number- -- https://gitlab.haskell.org/ghc/ghc/issues/16617---- |--- prop> nextDownDouble 1 == 0x1.ffff_ffff_ffff_fp-1--- prop> nextDownDouble (1/0) == maxFinite_ieee--- 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- | not (isIEEE x) || floatRadix x /= 2 || d /= 53 || expMin /= -1021 || expMax /= 1024 = error "rounded-hw assumes Double is IEEE binary64"- | otherwise = case castDoubleToWord64 x of- w | w .&. 0x7ff0_0000_0000_0000 == 0x7ff0_0000_0000_0000- , w /= 0x7ff0_0000_0000_0000 -> x -- NaN or negative infinity -> itself- 0x0000_0000_0000_0000 -> encodeFloat (-1) (expMin - d) -- +0 -> max negative- w | testBit w 63 -> castWord64ToDouble (w + 1) -- negative- | otherwise -> castWord64ToDouble (w - 1) -- positive- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,expMax) = floatRange x -- (-1021,1024) for Double---- |--- prop> nextTowardZeroFloat 1 == 0x1.fffffep-1--- prop> nextTowardZeroFloat (-1) == -0x1.fffffep-1--- prop> nextTowardZeroFloat (1/0) == maxFinite_ieee--- prop> nextTowardZeroFloat (-1/0) == -maxFinite_ieee--- prop> nextTowardZeroFloat 0 == 0--- prop> isNegativeZero (nextTowardZeroFloat (-0))--- prop> nextTowardZeroFloat 0x1p-149 == 0-nextTowardZeroFloat :: Float -> Float-nextTowardZeroFloat x- | not (isIEEE x) || floatRadix x /= 2 || d /= 24 || expMin /= -125 || expMax /= 128 = error "rounded-hw assumes Float is IEEE binary32"- | isNaN x || x == 0 = x -- NaN or zero -> itself- | otherwise = castWord32ToFloat (castFloatToWord32 x - 1) -- positive / negative- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,expMax) = floatRange x -- (-1021,1024) for Double- -- Note: castFloatToWord32 is buggy on GHC <= 8.8 on x86_64, so we can't use it to test for NaN or negative number- -- https://gitlab.haskell.org/ghc/ghc/issues/16617---- |--- prop> nextTowardZeroDouble 1 == 0x1.ffff_ffff_ffff_fp-1--- prop> nextTowardZeroDouble (-1) == -0x1.ffff_ffff_ffff_fp-1--- prop> nextTowardZeroDouble (1/0) == maxFinite_ieee--- prop> nextTowardZeroDouble (-1/0) == -maxFinite_ieee--- prop> nextTowardZeroDouble 0 == 0--- prop> isNegativeZero (nextTowardZeroDouble (-0))--- prop> nextTowardZeroDouble 0x1p-1074 == 0-nextTowardZeroDouble :: Double -> Double-nextTowardZeroDouble x- | not (isIEEE x) || floatRadix x /= 2 || d /= 53 || expMin /= -1021 || expMax /= 1024 = error "rounded-hw assumes Double is IEEE binary64"- | otherwise = case castDoubleToWord64 x of- w | w .&. 0x7ff0_0000_0000_0000 == 0x7ff0_0000_0000_0000- , w .&. 0x000f_ffff_ffff_ffff /= 0 -> x -- NaN -> itself- 0x8000_0000_0000_0000 -> x -- -0 -> itself- 0x0000_0000_0000_0000 -> x -- +0 -> itself- w -> castWord64ToDouble (w - 1) -- positive / negative- where- d, expMin :: Int- d = floatDigits x -- 53 for Double- (expMin,expMax) = floatRange x -- (-1021,1024) for Double--fusedMultiplyAdd :: RealFloat a => a -> a -> a -> a-fusedMultiplyAdd x y z- | isNaN x || isNaN y || isNaN z || isInfinite x || isInfinite y || isInfinite z = x * y + z- | otherwise = case toRational x * toRational y + toRational z of- 0 | isNegativeZero (x * y + z) -> -0- r -> fromRational r-{-# NOINLINE [1] fusedMultiplyAdd #-}--#ifdef USE_FFI--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- #-}--#endif+import Numeric.Floating.IEEE distanceUlp :: RealFloat a => a -> a -> Maybe Integer distanceUlp x y
− src/Numeric/Rounded/Hardware/Internal/RoundedResult.hs
@@ -1,54 +0,0 @@-{-# LANGUAGE DeriveFunctor #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE ScopedTypeVariables #-}-module Numeric.Rounded.Hardware.Internal.RoundedResult where-import Data.Proxy-import Data.Functor.Product-import Numeric.Rounded.Hardware.Internal.Rounding--class Functor f => Result f where- exact :: a -> f a- inexact :: a -- toward nearest- -> a -- toward inf- -> a -- toward neg inf- -> a -- toward zero- -> f a--newtype Exactness a = Exactness { getExactness :: Bool }- deriving (Eq, Ord, Show, Functor)--instance Rounding r => Result (Rounded r) where- exact x = Rounded x- inexact n inf ninf z = case rounding (Proxy :: Proxy r) of- ToNearest -> Rounded n- TowardInf -> Rounded inf- TowardNegInf -> Rounded ninf- TowardZero -> Rounded z--newtype DynamicRoundingMode a = DynamicRoundingMode { withRoundingMode :: RoundingMode -> a }- deriving (Functor)-instance Result DynamicRoundingMode where- exact x = DynamicRoundingMode (\_ -> x)- inexact n inf ninf z = DynamicRoundingMode $ \r ->- case r of- ToNearest -> n- TowardInf -> inf- TowardNegInf -> ninf- TowardZero -> z--instance Result Exactness where- exact _ = Exactness True- inexact _ _ _ _ = Exactness False---- Usage: Product (Rounded TowardNegInf) (Rounded TowardInf)-instance (Result f, Result g) => Result (Product f g) where- exact x = Pair (exact x) (exact x)- inexact n inf ninf z = Pair (inexact n inf ninf z) (inexact n inf ninf z)--newtype OppositeRoundingMode f a = OppositeRoundingMode { withOppositeRoundingMode :: f a }- deriving (Eq, Ord, Show, Functor)--instance Result f => Result (OppositeRoundingMode f) where- exact x = OppositeRoundingMode (exact x)- inexact n inf ninf z = OppositeRoundingMode (inexact n ninf inf z)
src/Numeric/Rounded/Hardware/Internal/Show.hs view
@@ -1,26 +1,15 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE ScopedTypeVariables #-} module Numeric.Rounded.Hardware.Internal.Show where-import Numeric.Rounded.Hardware.Internal.Rounding-import Data.Char (intToDigit)-import Data.Bifunctor (first)-import Data.Bits-import Math.NumberTheory.Logarithms+import Data.Bifunctor (first)+import Data.Bits+import Data.Char (intToDigit)+import Numeric.Floating.IEEE.Internal (countTrailingZerosInteger)+import Numeric.Rounded.Hardware.Internal.Rounding -- $setup -- >>> import Data.Int --- |--- prop> \x -> x == 0 || countTrailingZerosInteger (fromIntegral x) == countTrailingZeros (x :: Int64)--- >>> countTrailingZerosInteger 7--- 0--- >>> countTrailingZerosInteger 8--- 3-countTrailingZerosInteger :: Integer -> Int-countTrailingZerosInteger x- | x == 0 = error "countTrailingZerosInteger: zero"- | otherwise = integerLog2 (x `xor` (x - 1))- -- ratToDigitsRn :: RoundingMode -> Int -> Int -> Rational -> ([Int], Int) -- binaryFloatToDecimalDigitsRn _ prec x = ([d1,d2,...,dn], e)@@ -221,8 +210,8 @@ (d:ds) -> showString $ (intToDigit d : '.' : map intToDigit ds) ++ ('e' : show e') where padRight0 :: Int -> [Int] -> [Int]- padRight0 0 ys = ys- padRight0 !n [] = replicate n 0+ padRight0 0 ys = ys+ padRight0 !n [] = replicate n 0 padRight0 !n (y:ys) = y : padRight0 (n - 1) ys {-# SPECIALIZE showEFloatRn :: RoundingMode -> Maybe Int -> Double -> ShowS #-} @@ -243,32 +232,49 @@ | otherwise = case mprec of Nothing -> let (xs,e) = binaryFloatToDecimalDigits x l = length xs- in if e >= l- then if null xs- then showString "0.0"- else showString (map intToDigit xs ++ replicate (e - l) '0' ++ ".0")- else if e > 0 -- 0 < e < l- then if l == e -- null zs- then showString (map intToDigit xs ++ ".0")- else let (ys,zs) = splitAt (l - e) xs- ys' | null ys = [0]- | otherwise = ys- in showString (map intToDigit ys' ++ "." ++ map intToDigit zs)- else -- e < 0- showString ("0." ++ replicate (-e) '0' ++ map intToDigit xs)+ -- binaryFloatToDecimalDigits x = ([d1,d2,...,dl], e)+ -- x = 0.d1d2...dl * (10^^e)+ -- 0 <= di < 10+ in if e >= l then+ -- d1d2...dl<replicate (e-l) '0'>.0+ if null xs then+ showString "0.0"+ else+ showString (map intToDigit xs ++ replicate (e - l) '0' ++ ".0")+ else+ if e > 0 then -- 0 < e < l+ -- d1d2...d<e>.d<e+1>...dl+ if l == e then-- null zs+ showString (map intToDigit xs ++ ".0")+ else+ let (ys,zs) = splitAt e xs+ ys' = if null ys then [0] else ys+ in showString (map intToDigit ys' ++ "." ++ map intToDigit zs)+ else -- e < 0+ -- 0.<replicate (-e) '0'>d1d2...dl+ showString ("0." ++ replicate (-e) '0' ++ map intToDigit xs) Just prec -> let prec' = max prec 0 xs = binaryFloatToFixedDecimalDigitsRn r prec' x l = length xs- in if prec' == 0- then if null xs- then showString "0"- else showString $ map intToDigit xs- else if l <= prec'- then showString $ "0." ++ replicate (prec' - l) '0' ++ map intToDigit xs- else let (ys,zs) = splitAt (l - prec') xs- ys' | null ys = [0]- | otherwise = ys- in showString $ map intToDigit ys' ++ "." ++ map intToDigit zs+ -- binaryFloatToFixedDecimalDigitsRn _ prec' x = [d1,d2,...,dl]+ -- x = d1d2...dl * (10^^(-prec')) up to rounding+ -- 0 <= di < 10+ in if prec' == 0 then+ -- d1d2...dl or "0"+ if null xs then+ showString "0"+ else+ showString $ map intToDigit xs+ else+ if l <= prec' then+ -- 0.<replicate (prec'-l) '0'>d1d2...dl+ showString $ "0." ++ replicate (prec' - l) '0' ++ map intToDigit xs+ else+ -- l > prec'+ -- d1d2...d<l-prec'>.d<l-prec'+1>...dl+ let (ys,zs) = splitAt (l - prec') xs+ ys' = if null ys then [0] else ys+ in showString $ map intToDigit ys' ++ "." ++ map intToDigit zs {-# SPECIALIZE showFFloatRn :: RoundingMode -> Maybe Int -> Double -> ShowS #-} showGFloatRn :: RealFloat a => RoundingMode -> Maybe Int -> a -> ShowS
src/Numeric/Rounded/Hardware/Interval.hs view
@@ -315,3 +315,10 @@ y = indexByteArray (ByteArray byteArr) (2 * i + 1) in pairToInterval (x, y) -- unsafeReplace, unsafeAccum, unsafeAccumArray: Use default++{-# RULES+"fromIntegral/a->Interval Float"+ fromIntegral = \x -> case intervalFromIntegral x of (l, u) -> I l u :: Interval Float+"fromIntegral/a->Interval Double"+ fromIntegral = \x -> case intervalFromIntegral x of (l, u) -> I l u :: Interval Double+ #-}
src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs view
@@ -354,3 +354,10 @@ y = indexByteArray (ByteArray byteArr) (2 * i + 1) in pairToInterval (x, y) -- unsafeReplace, unsafeAccum, unsafeAccumArray: Use default++{-# RULES+"fromIntegral/a->Interval Float"+ fromIntegral = \x -> case intervalFromIntegral x of (l, u) -> I l u :: Interval Float+"fromIntegral/a->Interval Double"+ fromIntegral = \x -> case intervalFromIntegral x of (l, u) -> I l u :: Interval Double+ #-}
− test/FloatUtilSpec.hs
@@ -1,75 +0,0 @@-module FloatUtilSpec where-import Numeric.Rounded.Hardware.Internal-import Test.Hspec-import Test.Hspec.QuickCheck (prop)-import Test.QuickCheck-import Util (sameFloatP, variousFloats)--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, c_nextTowardZero :: (RealFloat a, CFloat a) => a -> a-c_nextUp x = c_nextafter x (1/0)-c_nextDown x = c_nextafter x (-1/0)-c_nextTowardZero x | isNegativeZero x = x- | otherwise = c_nextafter x 0--prop_nextUp_match :: (RealFloat a, CFloat a, Show a) => a -> Property-prop_nextUp_match x = nextUp x `sameFloatP` c_nextUp x--prop_nextDown_match :: (RealFloat a, CFloat a, Show a) => a -> Property-prop_nextDown_match x = nextDown x `sameFloatP` c_nextDown x--prop_nextTowardZero_match :: (RealFloat a, CFloat a, Show a) => a -> Property-prop_nextTowardZero_match x = nextTowardZero x `sameFloatP` c_nextTowardZero x--prop_fma_match :: (RealFloat a, CFloat a, Show a) => a -> a -> a -> Property-prop_fma_match x y z = fusedMultiplyAdd x y z `sameFloatP` c_fma x y z--isPositiveZero :: RealFloat a => a -> Bool-isPositiveZero x = x == 0 && not (isNegativeZero x)--prop_nextUp_nextDown :: (RealFloat a, Show 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) => a -> Property-prop_nextDown_nextUp x = x /= (1/0) ==>- let x' = nextDown (nextUp x)- in x' `sameFloatP` x .||. (isNegativeZero x .&&. isPositiveZero x')--spec :: Spec-spec = do- describe "Double" $ do- prop "nextUp vs C nextafter" $ forAll variousFloats (prop_nextUp_match :: Double -> Property)- prop "nextDown vs C nextafter" $ forAll variousFloats (prop_nextDown_match :: Double -> Property)- prop "nextTowardZero vs C nextafter" $ forAll variousFloats (prop_nextTowardZero_match :: Double -> Property)- prop "nextUp . nextDown == id (unless -inf)" $ forAll variousFloats (prop_nextUp_nextDown :: Double -> Property)- prop "nextDown . nextUp == id (unless inf)" $ forAll variousFloats (prop_nextDown_nextUp :: Double -> Property)- prop "fusedMultiplyAdd vs C fma" $ forAll variousFloats (prop_fma_match :: Double -> Double -> Double -> Property)- describe "Float" $ do- prop "nextUp vs C nextafter" $ forAll variousFloats (prop_nextUp_match :: Float -> Property)- prop "nextDown vs C nextafter" $ forAll variousFloats (prop_nextDown_match :: Float -> Property)- prop "nextTowardZero vs C nextafter" $ forAll variousFloats (prop_nextTowardZero_match :: Float -> Property)- prop "nextUp . nextDown == id (unless -inf)" $ forAll variousFloats (prop_nextUp_nextDown :: Float -> Property)- prop "nextDown . nextUp == id (unless inf)" $ forAll variousFloats (prop_nextDown_nextUp :: Float -> Property)- prop "fusedMultiplyAdd vs C fma" $ forAll variousFloats (prop_fma_match :: Float -> Float -> Float -> Property)
test/FromIntegerSpec.hs view
@@ -19,32 +19,32 @@ prop_roundedFromInteger_check :: forall a. (RealFloat a, RoundedRing a) => Proxy a -> RoundingMode -> Integer -> Property prop_roundedFromInteger_check _proxy r x = (roundedFromInteger r x :: a)- `sameFloatP` (fromInt r x :: a)+ `sameFloatP` (roundedFromInteger_default r x :: a) prop_roundedFromInt64_check :: forall a. (RealFloat a, RoundedRing a) => Proxy a -> RoundingMode -> Int64 -> Property prop_roundedFromInt64_check _proxy r x = (roundedFromInteger r (fromIntegral x) :: a)- `sameFloatP` (fromInt r (fromIntegral x) :: a)+ `sameFloatP` (roundedFromInteger_default r (fromIntegral x) :: a) prop_roundedFromWord64_check :: forall a. (RealFloat a, RoundedRing a) => Proxy a -> RoundingMode -> Word64 -> Property prop_roundedFromWord64_check _proxy r x = (roundedFromInteger r (fromIntegral x) :: a)- `sameFloatP` (fromInt r (fromIntegral x) :: a)+ `sameFloatP` (roundedFromInteger_default r (fromIntegral x) :: a) -prop_fromInt_order :: forall a. RealFloat a => Proxy a -> Integer -> Property-prop_fromInt_order _proxy x- = let ne = fromInt ToNearest x :: a- ze = fromInt TowardZero x :: a- inf = fromInt TowardInf x :: a- ninf = fromInt TowardNegInf x :: a+prop_roundedFromInteger_order :: forall a. (RealFloat a, RoundedRing a) => Proxy a -> Integer -> Property+prop_roundedFromInteger_order _proxy x+ = let ne = roundedFromInteger ToNearest x :: a+ ze = roundedFromInteger TowardZero x :: a+ inf = roundedFromInteger TowardInf x :: a+ ninf = roundedFromInteger TowardNegInf x :: a in ninf <= inf .&&. (ne == ninf || ne == inf) .&&. (if x < 0 then ze == inf else ze == ninf) -prop_fromInt_exact :: forall a. RealFloat a => Proxy a -> Integer -> Property-prop_fromInt_exact _proxy x- = let inf = fromInt TowardInf x :: a- ninf = fromInt TowardNegInf x :: a+prop_roundedFromInteger_exact :: forall a. (RealFloat a, RoundedRing a) => Proxy a -> Integer -> Property+prop_roundedFromInteger_exact _proxy x+ = let inf = roundedFromInteger TowardInf x :: a+ ninf = roundedFromInteger TowardNegInf x :: a in if ninf == inf then not (isInfinite inf) .&&. toRational inf === fromInteger x else if isInfinite inf@@ -71,9 +71,9 @@ prop "roundedFromInteger/Word64" $ \r -> prop_roundedFromWord64_check proxy r prop "order" $- forAllShrink variousIntegers shrinkIntegral (prop_fromInt_order proxy)+ forAllShrink variousIntegers shrinkIntegral (prop_roundedFromInteger_order proxy) prop "exactness" $- forAllShrink variousIntegers shrinkIntegral (prop_fromInt_exact proxy)+ forAllShrink variousIntegers shrinkIntegral (prop_roundedFromInteger_exact proxy) spec :: Spec spec = do
test/FromRationalSpec.hs view
@@ -17,23 +17,23 @@ prop_roundedFromRational_check :: forall a. (RealFloat a, RoundedFractional a) => Proxy a -> RoundingMode -> Rational -> Property prop_roundedFromRational_check _proxy r x- = (fromRatio r (numerator x) (denominator x) :: a) -- the standard implementation+ = (roundedFromRational_default r x :: a) -- the standard implementation `sameFloatP` (roundedFromRational r x :: a) -- may be optimized prop_fromRatio_order :: forall a. RealFloat a => Proxy a -> Rational -> Property prop_fromRatio_order _proxy x- = let ne = fromRatio ToNearest (numerator x) (denominator x) :: a- ze = fromRatio TowardZero (numerator x) (denominator x) :: a- inf = fromRatio TowardInf (numerator x) (denominator x) :: a- ninf = fromRatio TowardNegInf (numerator x) (denominator x) :: a+ = let ne = roundedFromRational_default ToNearest x :: a+ ze = roundedFromRational_default TowardZero x :: a+ inf = roundedFromRational_default TowardInf x :: a+ ninf = roundedFromRational_default TowardNegInf x :: a in ninf <= inf .&&. (ne == ninf || ne == inf) .&&. (if x < 0 then ze == inf else ze == ninf) prop_fromRatio_exact :: forall a. RealFloat a => Proxy a -> Rational -> Property prop_fromRatio_exact _proxy x- = let inf = fromRatio TowardInf (numerator x) (denominator x) :: a- ninf = fromRatio TowardNegInf (numerator x) (denominator x) :: a+ = let inf = roundedFromRational_default TowardInf x :: a+ ninf = roundedFromRational_default TowardNegInf x :: a in if ninf == inf then not (isInfinite inf) .&&. toRational inf === x else if isInfinite inf@@ -47,7 +47,7 @@ else toRational inf =/= x .&&. toRational ninf =/= x -specT :: forall a. (RealFloat a, RoundedFractional a) => Proxy a -> Bool -> Spec+specT :: (RealFloat a, RoundedFractional a) => Proxy a -> Bool -> Spec specT proxy checkAgainstStock = do when checkAgainstStock $ do -- Although fromRational for Double/Float correctly round to nearest, other types may not.
test/ShowFloatSpec.hs view
@@ -69,15 +69,31 @@ prop "showGFloat" $ prop_showGFloat proxy -- 0.5 should be exactly representable in the type...- prop "showFFloatRn Nothing 0.5" $ \r -> showFFloatRn r Nothing (0.5 :: a) "" === "0.5"- prop "showFFloatRn (Just 0) 0.5" $ \r -> showFFloatRn r (Just 0) (0.5 :: a) "" === (if r == TowardInf then "1" else "0")- prop "showFFloatRn (Just 3) 0.5" $ \r -> showFFloatRn r (Just 3) (0.5 :: a) "" === "0.500"- prop "showGFloatRn Nothing 0.5" $ \r -> showGFloatRn r Nothing (0.5 :: a) "" === "0.5"- prop "showGFloatRn (Just 0) 0.5" $ \r -> showGFloatRn r (Just 0) (0.5 :: a) "" === (if r == TowardInf then "1" else "0")- prop "showGFloatRn (Just 3) 0.5" $ \r -> showGFloatRn r (Just 3) (0.5 :: a) "" === "0.500"- prop "showEFloatRn Nothing 0.5" $ \r -> showEFloatRn r Nothing (0.5 :: a) "" === "5.0e-1"- prop "showEFloatRn (Just 0) 0.5" $ \r -> showEFloatRn r (Just 0) (0.5 :: a) "" === "5e-1"- prop "showEFloatRn (Just 3) 0.5" $ \r -> showEFloatRn r (Just 3) (0.5 :: a) "" === "5.000e-1"+ do let x = 0.5 `asProxyTypeOf` proxy+ prop "showFFloatRn Nothing 0.5" $ \r -> showFFloatRn r Nothing x "" === "0.5"+ prop "showFFloatRn (Just 0) 0.5" $ \r -> showFFloatRn r (Just 0) x "" === (if r == TowardInf then "1" else "0")+ prop "showFFloatRn (Just 3) 0.5" $ \r -> showFFloatRn r (Just 3) x "" === "0.500"+ prop "showGFloatRn Nothing 0.5" $ \r -> showGFloatRn r Nothing x "" === "0.5"+ prop "showGFloatRn (Just 0) 0.5" $ \r -> showGFloatRn r (Just 0) x "" === (if r == TowardInf then "1" else "0")+ prop "showGFloatRn (Just 3) 0.5" $ \r -> showGFloatRn r (Just 3) x "" === "0.500"+ prop "showEFloatRn Nothing 0.5" $ \r -> showEFloatRn r Nothing x "" === "5.0e-1"+ prop "showEFloatRn (Just 0) 0.5" $ \r -> showEFloatRn r (Just 0) x "" === "5e-1"+ prop "showEFloatRn (Just 3) 0.5" $ \r -> showEFloatRn r (Just 3) x "" === "5.000e-1"++ -- -17.5625 should be exactly representable in the type...+ do let x = (-17.5625) `asProxyTypeOf` proxy+ prop "showFFloatRn Nothing -17.5625" $ \r -> showFFloatRn r Nothing x "" === "-17.5625"+ prop "showFFloatRn (Just 0) -17.5625" $ \r -> showFFloatRn r (Just 0) x "" === (if r == TowardInf || r == TowardZero then "-17" else "-18")+ prop "showFFloatRn (Just 3) -17.5625" $ \r -> showFFloatRn r (Just 3) x "" === (if r == TowardNegInf then "-17.563" else "-17.562")+ prop "showFFloatRn (Just 6) -17.5625" $ \r -> showFFloatRn r (Just 6) x "" === "-17.562500"+ prop "showGFloatRn Nothing -17.5625" $ \r -> showGFloatRn r Nothing x "" === "-17.5625"+ prop "showGFloatRn (Just 0) -17.5625" $ \r -> showGFloatRn r (Just 0) x "" === (if r == TowardInf || r == TowardZero then "-17" else "-18")+ prop "showGFloatRn (Just 3) -17.5625" $ \r -> showGFloatRn r (Just 3) x "" === (if r == TowardNegInf then "-17.563" else "-17.562")+ prop "showGFloatRn (Just 6) -17.5625" $ \r -> showGFloatRn r (Just 6) x "" === "-17.562500"+ prop "showEFloatRn Nothing -17.5625" $ \r -> showEFloatRn r Nothing x "" === "-1.75625e1"+ prop "showEFloatRn (Just 0) -17.5625" $ \r -> showEFloatRn r (Just 0) x "" === (if r == TowardInf || r == TowardZero then "-1e1" else "-2e1")+ prop "showEFloatRn (Just 3) -17.5625" $ \r -> showEFloatRn r (Just 3) x "" === (if r == TowardNegInf then "-1.757e1" else "-1.756e1")+ prop "showEFloatRn (Just 6) -17.5625" $ \r -> showEFloatRn r (Just 6) x "" === "-1.756250e1" spec :: Spec spec = do
test/Spec.hs view
@@ -9,7 +9,6 @@ import Numeric.Rounded.Hardware.Backend (backendName) import qualified RoundedArithmeticSpec import qualified ShowFloatSpec-import qualified FloatUtilSpec import Test.Hspec import qualified VectorSpec #ifdef TEST_X87_LONG_DOUBLE@@ -41,7 +40,6 @@ describe "showFloat" ShowFloatSpec.spec describe "rounded arithmetic" RoundedArithmeticSpec.spec describe "interval arithmetic" IntervalArithmeticSpec.spec- describe "FloatUtil" FloatUtilSpec.spec describe "Vector" VectorSpec.spec describe "Constants" ConstantsSpec.spec #ifdef TEST_X87_LONG_DOUBLE
test/X87LongDoubleSpec.hs view
@@ -3,7 +3,6 @@ import qualified ConstantsSpec import Data.Int import Data.Proxy-import qualified FloatUtilSpec import qualified FromIntegerSpec import qualified FromRationalSpec import qualified IntervalArithmeticSpec@@ -37,8 +36,6 @@ describe "fromRational" $ FromRationalSpec.specT ldProxy False describe "showFloat" $ ShowFloatSpec.specT ldProxy describe "constants" $ ConstantsSpec.specT ldProxy- prop "nextUp . nextDown == id (unless -inf)" $ forAll variousFloats (FloatUtilSpec.prop_nextUp_nextDown :: LongDouble -> Property)- prop "nextDown . nextUp == id (unless inf)" $ forAll variousFloats (FloatUtilSpec.prop_nextDown_nextUp :: LongDouble -> Property) where ldProxy :: Proxy LongDouble ldProxy = Proxy