{-|
This module provides various combinations of explicit floating point operations.
The are three supported rounding variants:
* ceil (rounding towards positive infinity)
* floor (rounding towards negative infinity)
* truncate (rounding towards zero)
Operators:
* `+` add
* `-` subtract
* `*` multiplicate
* `/` divide
* `sqrt` squareRoot
* `id` converting numbers between formats
The floating point data types:
* `Float` (single precision)
* `Double` (double precision)
The behaviour for `NaN`s is very hardware dependent. In case an operation would result in a `NaN` the corresponding round-to-nearest-even-on-tie variant is used to get the result.
-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ExtendedLiterals #-}
module Rounding
-- * Arithmetic Signum
( i32_arithmic_signum_f64
, i32_arithmic_signum_f64#
, i32_arithmic_signum_f32
, i32_arithmic_signum_f32#
-- * Sign bit
, i32_sign_bit_f64
, i32_sign_bit_f64#
, i32_sign_bit_f32
, i32_sign_bit_f32#
-- * Bit pattern neighbours
, f64_successorIEEE
, f64_successorIEEE#
, f64_predecessorIEEE
, f64_predecessorIEEE#
, f32_successorIEEE
, f32_successorIEEE#
, f32_predecessorIEEE
, f32_predecessorIEEE#
-- * Ceil variants
{-|
Operators of the ceil variant operate as if they calculated the exact mathematic result. Then they take the minimum of all floating point numbers that are greater or equal to the exact result. Special care is taken for the sign of zero.
Should the result from a plus or minus operation be zero the sign is positive.
-}
, f32_squareRoot_ceil
, f32_add_ceil
, f32_subtract_ceil
, f32_multiplicate_ceil
, f32_divide_ceil
, f64_squareRoot_ceil
, f64_add_ceil
, f64_subtract_ceil
, f64_multiplicate_ceil
, f64_divide_ceil
, f32_convert_i32_signed_ceil
, f32_convert_i32_unsigned_ceil
, f32_convert_i64_signed_ceil
, f32_convert_i64_unsigned_ceil
, f32_demote_f64_ceil
, f64_convert_i32_signed_ceil
, f64_convert_i32_unsigned_ceil
, f64_convert_i64_signed_ceil
, f64_convert_i64_unsigned_ceil
, f64_promote_f32_ceil
-- * Floor variants
{-|
The operators operate as if they calculate the exact mathematic result. Then they take the maximum of all floating point numbers that are less or equal to the exact result.
Should the result from a plus or minus operation be zero the sign is negative.
-}
, f32_squareRoot_floor
, f32_add_floor
, f32_subtract_floor
, f32_multiplicate_floor
, f32_divide_floor
, f64_squareRoot_floor
, f64_add_floor
, f64_subtract_floor
, f64_multiplicate_floor
, f64_divide_floor
, f32_convert_i32_signed_floor
, f32_convert_i32_unsigned_floor
, f32_convert_i64_signed_floor
, f32_convert_i64_unsigned_floor
, f32_demote_f64_floor
, f64_convert_i32_signed_floor
, f64_convert_i32_unsigned_floor
, f64_convert_i64_signed_floor
, f64_convert_i64_unsigned_floor
, f64_promote_f32_floor
-- * Truncate variants
{-|
Operators of the truncate variant operate as if they calculated the exact result. Then they take the result from the floor variant in case the exact result is greater than zero. Should the exact result be zero or less the ceil variant is choosen.
-}
, f32_squareRoot_truncate
, f32_add_truncate
, f32_subtract_truncate
, f32_multiplicate_truncate
, f32_divide_truncate
, f64_squareRoot_truncate
, f64_add_truncate
, f64_subtract_truncate
, f64_multiplicate_truncate
, f64_divide_truncate
, f32_convert_i32_signed_truncate
, f32_convert_i32_unsigned_truncate
, f32_convert_i64_signed_truncate
, f32_convert_i64_unsigned_truncate
, f32_demote_f64_truncate
, f64_convert_i32_signed_truncate
, f64_convert_i32_unsigned_truncate
, f64_convert_i64_signed_truncate
, f64_convert_i64_unsigned_truncate
, f64_promote_f32_truncate
) where
import GHC.Exts
import GHC.Float
import GHC.Int
import Data.Word
-- | There are three possible outputs:
--
-- prop> i32_arithmic_signum_f64 (-1.0) = -1
-- prop> i32_arithmic_signum_f64 (-0.0) = 0
-- prop> i32_arithmic_signum_f64 (0.0) = 0
-- prop> i32_arithmic_signum_f64 (1.0) = 1
--
-- @`NaN`@ values are mapped to @`0`@:
--
-- prop> i32_arithmic_signum_f64 (0/0) = 0
{-# INLINABLE i32_arithmic_signum_f64 #-}
i32_arithmic_signum_f64 :: Double -> Int32
i32_arithmic_signum_f64 (D# x) = I32# (intToInt32# (i32_arithmic_signum_f64# x))
{-# INLINABLE i32_arithmic_signum_f32 #-}
i32_arithmic_signum_f32 :: Float -> Int32
i32_arithmic_signum_f32 (F# x) = I32# (intToInt32# (i32_arithmic_signum_f32# x))
-- https://stackoverflow.com/questions/1903954/is-there-a-standard-sign-function-signum-sgn-in-c-c
{-# INLINABLE i32_arithmic_signum_f64# #-}
i32_arithmic_signum_f64# :: Double# -> Int#
i32_arithmic_signum_f64# x = (0.0## <## x) -# (x <## 0.0##)
{-# INLINABLE i32_arithmic_signum_f32# #-}
i32_arithmic_signum_f32# :: Float# -> Int#
i32_arithmic_signum_f32# x = (0.0# `ltFloat#` x) -# (x `ltFloat#` 0.0#)
-- | There are two possible outputs:
--
-- prop> i32_sign_bit_f64 (-1.0) = 1
-- prop> i32_sign_bit_f64 (-0.0) = 1
-- prop> i32_sign_bit_f64 (0.0) = 0
-- prop> i32_sign_bit_f64 (1.0) = 0
--
-- All @`NaN`@s have a proper sign:
--
-- prop> 1 = i32_sign_bit_f64 (0/0)
-- prop> 0 = i32_sign_bit_f64 (castWord64ToDouble 0b0111111111111000000000000000000000000000000000000000000000000000)
{-# INLINABLE i32_sign_bit_f64 #-}
i32_sign_bit_f64 :: Double -> Int32
i32_sign_bit_f64 (D# x) = I32# (i32_sign_bit_f64# x)
{-# INLINABLE i32_sign_bit_f64# #-}
i32_sign_bit_f64# :: Double# -> Int32#
i32_sign_bit_f64# value = result
where
value_as_word64 :: Word64#
value_as_word64 = stgDoubleToWord64 value
sign_bit :: Word64#
sign_bit
= 0b1000000000000000000000000000000000000000000000000000000000000000#Word64
sign_bit_of_value :: Word64#
sign_bit_of_value = and64# sign_bit value_as_word64
result
= case sign_bit_of_value of
0b1000000000000000000000000000000000000000000000000000000000000000#Word64
-> 0b1#Int32
0b0000000000000000000000000000000000000000000000000000000000000000#Word64
-> 0b0#Int32
_ -> 0b0#Int32 -- error "cannot happen"
{-# INLINABLE i32_sign_bit_f32 #-}
i32_sign_bit_f32 :: Float -> Int32
i32_sign_bit_f32 (F# x) = I32# (i32_sign_bit_f32# x)
{-# INLINABLE i32_sign_bit_f32# #-}
i32_sign_bit_f32# :: Float# -> Int32#
i32_sign_bit_f32# value = result
where
value_as_word32 :: Word32#
value_as_word32 = stgFloatToWord32 value
sign_bit :: Word32#
sign_bit
= 0b10000000000000000000000000000000#Word32
sign_bit_of_value :: Word32#
sign_bit_of_value = andWord32# sign_bit value_as_word32
result
= case sign_bit_of_value of
0b10000000000000000000000000000000#Word32
-> 0b1#Int32
0b00000000000000000000000000000000#Word32
-> 0b0#Int32
_ -> 0b0#Int32 -- error "cannot happen"
-----------
{-# INLINABLE f64_successorIEEE# #-}
f64_successorIEEE#
:: Double#
-> Double#
f64_successorIEEE#
value
= result
where
result = resulting_double
value_as_word64 :: Word64#
value_as_word64 = stgDoubleToWord64 value
sign_bit :: Word64#
sign_bit
= 0b1000000000000000000000000000000000000000000000000000000000000000#Word64
exponent_bits :: Word64#
exponent_bits
= 0b0111111111110000000000000000000000000000000000000000000000000000#Word64
_negativ_infinity :: Word64#
_negativ_infinity
= 0b1111111111110000000000000000000000000000000000000000000000000000#Word64
_least_finit_value :: Word64#
_least_finit_value
= 0b1111111111101111111111111111111111111111111111111111111111111111#Word64
-- -1.7976931348623157e308##
_least_positiv_value :: Word64#
_least_positiv_value
= 0b0000000000000000000000000000000000000000000000000000000000000001#Word64
-- 5.0e-324##
_negative_zero = sign_bit
exponent_bits_of_value :: Word64#
exponent_bits_of_value = and64# exponent_bits value_as_word64
sign_bit_of_value :: Word64#
sign_bit_of_value = and64# sign_bit value_as_word64
resulting_double
= case exponent_bits_of_value of
{- not real value or not negativ zero -}
0b0111111111110000000000000000000000000000000000000000000000000000#Word64 {-exponent_bits-}
-> case value_as_word64 of
0b1111111111110000000000000000000000000000000000000000000000000000#Word64 {-negativ_infinity-}
-> -1.7976931348623157e308## {-least_finit_value-}
_
-> value --not real values stay as they are (except for negativ_infinity)
{- real value or negativ zero -}
_
-> case sign_bit_of_value of
{-positive real value-}
0b0000000000000000000000000000000000000000000000000000000000000000#Word64
-> stgWord64ToDouble (plusWord64# 1#Word64 value_as_word64)
{-negative real value-}
_
-> case value_as_word64 of
{-negative_zero-}
0b1000000000000000000000000000000000000000000000000000000000000000#Word64 {-negative_zero-}
-> 5.0e-324## {-least_positiv_value-}
{-proper negativ value-}
_
-> stgWord64ToDouble (subWord64# value_as_word64 1#Word64)
{-|
This function is supposed to be a pure haskell replacement for `Numeric.IEEE.succIEEE`. ([succIEEE](https://hackage.haskell.org/package/ieee754-0.8.0/docs/Numeric-IEEE.html#v:succIEEE))
Floating point numbers of the same sign have canonical ordered bitpatterns. This means that first a float that represents a real number can first be reinterpeted as an integer. Then one can add or subtract `1` inorder to get the two floating point neighbours. Special cases are infinit values and negative and positive zero. @`NaN`@s are returned unchanged.
prop> f64_successorIEEE (-1/0) = -1.7976931348623157e308
prop> f64_successorIEEE (-0.0) = castWord64ToDouble 1
prop> f64_successorIEEE (0.0) = castWord64ToDouble 1
prop> f64_successorIEEE (castWord64ToDouble 1) = castWord64ToDouble 2
prop> f64_successorIEEE (castWord64ToDouble 2) = castWord64ToDouble 3
prop> f64_successorIEEE (1.7976931348623157e308) = 1/0
-}
{-# INLINABLE f64_successorIEEE #-}
f64_successorIEEE
:: Double
-> Double
f64_successorIEEE
(D# value#)
= D# (f64_successorIEEE# value#)
{-# INLINABLE f64_predecessorIEEE# #-}
f64_predecessorIEEE#
:: Double#
-> Double#
f64_predecessorIEEE#
value
= symetric_result
where
symetric_result
= negateDouble#
$# f64_successorIEEE#
$# negateDouble#
$# value
infixr 0 $#
($#) :: (Double# -> Double#) -> Double# -> Double#
f $# x = f x
{-# INLINABLE f64_predecessorIEEE #-}
f64_predecessorIEEE
:: Double
-> Double
f64_predecessorIEEE
(D# value#)
= D# (f64_predecessorIEEE# value#)
------------
{-# INLINABLE f32_successorIEEE# #-}
f32_successorIEEE#
:: Float#
-> Float#
f32_successorIEEE#
value
= result
where
result = resulting_float
value_as_word64 :: Word32# --TODO recent GHC returns Word32#
value_as_word64 = stgFloatToWord32 value
sign_bit :: Word32#
sign_bit
= 0b10000000000000000000000000000000#Word32
exponent_bits :: Word32#
exponent_bits
= 0b01111111100000000000000000000000#Word32
_negativ_infinity :: Word32#
_negativ_infinity
= 0b11111111100000000000000000000000#Word32
_least_finit_value :: Word32#
_least_finit_value
= 0b11111111011111111111111111111111#Word32
-- -3.4028235e38#
_least_positiv_value :: Word32#
_least_positiv_value
= 0b00000000000000000000000000000001#Word32
-- 1.0e-45#
_negative_zero = sign_bit
exponent_bits_of_value :: Word32#
exponent_bits_of_value = andWord32# exponent_bits value_as_word64
sign_bit_of_value :: Word32#
sign_bit_of_value = andWord32# sign_bit value_as_word64
resulting_float
= case exponent_bits_of_value of
{- not real value or not negativ zero -}
0b01111111100000000000000000000000#Word32 {-exponent_bits-}
-> case value_as_word64 of
0b11111111100000000000000000000000#Word32 {-negativ_infinity-}
-> -3.4028235e38# {-least_finit_value-}
_
-> value --not real values stay as they are (except for negativ_infinity)
{- real value or negativ zero -}
_
-> case sign_bit_of_value of
{-positive real value-}
0b00000000000000000000000000000000#Word32
-> stgWord32ToFloat (plusWord32# 1#Word32 value_as_word64)
{-negative real value-}
_
-> case value_as_word64 of
{-negative_zero-}
0b10000000000000000000000000000000#Word32 {-negative_zero-}
-> 1.0e-45# {-least_positiv_value-}
{-proper negativ value-}
_
-> stgWord32ToFloat (minusWord32# value_as_word64 1#Word32)
minusWord32# :: Word32# -> Word32# -> Word32#
minusWord32# x y = wordToWord32# (minusWord# (word32ToWord# x) (word32ToWord# y))
{-# INLINABLE f32_successorIEEE #-}
f32_successorIEEE
:: Float
-> Float
f32_successorIEEE
(F# value#)
= F# (f32_successorIEEE# value#)
{-# INLINABLE f32_predecessorIEEE# #-}
f32_predecessorIEEE#
:: Float#
-> Float#
f32_predecessorIEEE#
value
= result
where
symetric_result
= negateFloat#
$# f32_successorIEEE#
$# negateFloat#
$# value
infixr 0 $#
($#) :: (Float# -> Float#) -> Float# -> Float#
f $# x = f x
result
= if isNegativeZero (F# symetric_result)
then 0.0#
else symetric_result
{-# INLINABLE f32_predecessorIEEE #-}
f32_predecessorIEEE
:: Float
-> Float
f32_predecessorIEEE
(F# value#)
= F# (f32_predecessorIEEE# value#)
-----------
-- |
-- prop> (castWord64ToDouble 1) == f32_demote_f64_ceil (castWord64ToDouble 1)
f32_demote_f64_ceil :: Double -> Float
f32_demote_f64_ceil = sanatizeNaNDemotion $ sanatizeNegtiveZeroDemotionPromotion $ convert RoundUp
f32_demote_f64_floor :: Double -> Float
f32_demote_f64_floor = sanatizeNaNDemotion $ sanatizeNegtiveZeroDemotionPromotion $ convert RoundDown
f32_demote_f64_truncate :: Double -> Float
f32_demote_f64_truncate = sanatizeNaNDemotion $ sanatizeNegtiveZeroDemotionPromotion $ convert Truncate
f64_promote_f32_ceil :: Float -> Double
f64_promote_f32_ceil = sanatizeNaNPromotion $ sanatizeNegtiveZeroDemotionPromotion $ convert RoundUp
f64_promote_f32_floor :: Float -> Double
f64_promote_f32_floor = sanatizeNaNPromotion $ sanatizeNegtiveZeroDemotionPromotion $ convert RoundDown
f64_promote_f32_truncate :: Float -> Double
f64_promote_f32_truncate = sanatizeNaNPromotion $ sanatizeNegtiveZeroDemotionPromotion $ convert Truncate
-- |
-- prop> 0.0 = f32_squareRoot_ceil 0.0
-- prop> 1.0 = f32_squareRoot_ceil 1.0
-- prop> (1/0) = f32_squareRoot_ceil (1/0)
-- prop> True = isNaN (f32_squareRoot_ceil (-1/0))
f32_squareRoot_ceil :: Float -> Float
f32_squareRoot_ceil = squareRoot RoundUp
f32_squareRoot_floor :: Float -> Float
f32_squareRoot_floor = squareRoot RoundDown
-- |
-- prop> f32_squareRoot_truncate = f32_squareRoot_floor
f32_squareRoot_truncate :: Float -> Float
f32_squareRoot_truncate = squareRoot Truncate
f64_squareRoot_ceil :: Double -> Double
f64_squareRoot_ceil = squareRoot RoundUp
f64_squareRoot_floor :: Double -> Double
f64_squareRoot_floor = squareRoot RoundDown
f64_squareRoot_truncate :: Double -> Double
f64_squareRoot_truncate = squareRoot Truncate
f32_add_ceil :: Float -> Float -> Float
f32_add_ceil = binary_operator_rounded RoundUp (+)
-- |
-- prop> (-0.0) = f32_add_floor 1 (-1)
-- prop> (-0.0) = f32_add_floor (-1) 1
f32_add_floor :: Float -> Float -> Float
f32_add_floor = sanitizeZeroForAddition $ binary_operator_rounded RoundDown (+)
f32_add_truncate :: Float -> Float -> Float
f32_add_truncate = binary_operator_rounded Truncate (+)
-- |
-- prop> 0.0 == f32_subtract_ceil 1 1
f32_subtract_ceil :: Float -> Float -> Float
f32_subtract_ceil = binary_operator_rounded RoundUp (-)
f32_subtract_floor :: Float -> Float -> Float
f32_subtract_floor = sanitizeZeroForSubtraction $ binary_operator_rounded RoundDown (-)
f32_subtract_truncate :: Float -> Float -> Float
f32_subtract_truncate = binary_operator_rounded Truncate (-)
f32_multiplicate_ceil :: Float -> Float -> Float
f32_multiplicate_ceil = binary_operator_rounded RoundUp (*)
f32_multiplicate_floor :: Float -> Float -> Float
f32_multiplicate_floor = binary_operator_rounded RoundDown (*)
f32_multiplicate_truncate :: Float -> Float -> Float
f32_multiplicate_truncate = binary_operator_rounded Truncate (*)
f32_divide_ceil :: Float -> Float -> Float
f32_divide_ceil = binary_operator_rounded RoundUp (/)
f32_divide_floor :: Float -> Float -> Float
f32_divide_floor = binary_operator_rounded RoundDown (/)
f32_divide_truncate :: Float -> Float -> Float
f32_divide_truncate = binary_operator_rounded Truncate (/)
-- |
-- prop> 0.0 = f64_add_ceil 1 (-1)
f64_add_ceil :: Double -> Double -> Double
f64_add_ceil = binary_operator_rounded RoundUp (+)
f64_add_floor :: Double -> Double -> Double
f64_add_floor = sanitizeZeroForAddition $ binary_operator_rounded RoundDown (+)
f64_add_truncate :: Double -> Double -> Double
f64_add_truncate = binary_operator_rounded Truncate (+)
f64_subtract_ceil :: Double -> Double -> Double
f64_subtract_ceil = binary_operator_rounded RoundUp (-)
f64_subtract_floor :: Double -> Double -> Double
f64_subtract_floor = sanitizeZeroForSubtraction $ binary_operator_rounded RoundDown (-)
f64_subtract_truncate :: Double -> Double -> Double
f64_subtract_truncate = binary_operator_rounded Truncate (-)
f64_multiplicate_ceil :: Double -> Double -> Double
f64_multiplicate_ceil = binary_operator_rounded RoundUp (*)
f64_multiplicate_floor :: Double -> Double -> Double
f64_multiplicate_floor = binary_operator_rounded RoundDown (*)
f64_multiplicate_truncate :: Double -> Double -> Double
f64_multiplicate_truncate = binary_operator_rounded Truncate (*)
f64_divide_ceil :: Double -> Double -> Double
f64_divide_ceil = binary_operator_rounded RoundUp (/)
f64_divide_floor :: Double -> Double -> Double
f64_divide_floor = binary_operator_rounded RoundDown (/)
f64_divide_truncate :: Double -> Double -> Double
f64_divide_truncate = binary_operator_rounded Truncate (/)
f32_convert_i32_unsigned_ceil :: Word32 -> Float
f32_convert_i32_unsigned_ceil = convert RoundUp
f32_convert_i32_signed_ceil :: Int32 -> Float
f32_convert_i32_signed_ceil = convert RoundUp
f32_convert_i32_unsigned_floor :: Word32 -> Float
f32_convert_i32_unsigned_floor = convert RoundDown
f32_convert_i32_signed_floor :: Int32 -> Float
f32_convert_i32_signed_floor = convert RoundDown
f32_convert_i32_unsigned_truncate :: Word32 -> Float
f32_convert_i32_unsigned_truncate = convert Truncate
f32_convert_i32_signed_truncate :: Int32 -> Float
f32_convert_i32_signed_truncate = convert Truncate
f32_convert_i64_unsigned_ceil :: Word64 -> Float
f32_convert_i64_unsigned_ceil = convert RoundUp
f32_convert_i64_signed_ceil :: Int64 -> Float
f32_convert_i64_signed_ceil = convert RoundUp
f32_convert_i64_unsigned_floor :: Word64 -> Float
f32_convert_i64_unsigned_floor = convert RoundDown
f32_convert_i64_signed_floor :: Int64 -> Float
f32_convert_i64_signed_floor = convert RoundDown
f32_convert_i64_unsigned_truncate :: Word64 -> Float
f32_convert_i64_unsigned_truncate = convert Truncate
f32_convert_i64_signed_truncate :: Int64 -> Float
f32_convert_i64_signed_truncate = convert Truncate
f64_convert_i32_unsigned_ceil :: Word32 -> Double
f64_convert_i32_unsigned_ceil = convert RoundUp
f64_convert_i32_signed_ceil :: Int32 -> Double
f64_convert_i32_signed_ceil = convert RoundUp
f64_convert_i32_unsigned_floor :: Word32 -> Double
f64_convert_i32_unsigned_floor = convert RoundDown
f64_convert_i32_signed_floor :: Int32 -> Double
f64_convert_i32_signed_floor = convert RoundDown
f64_convert_i32_unsigned_truncate :: Word32 -> Double
f64_convert_i32_unsigned_truncate = convert Truncate
f64_convert_i32_signed_truncate :: Int32 -> Double
f64_convert_i32_signed_truncate = convert Truncate
f64_convert_i64_unsigned_ceil :: Word64 -> Double
f64_convert_i64_unsigned_ceil = convert RoundUp
f64_convert_i64_signed_ceil :: Int64 -> Double
f64_convert_i64_signed_ceil = convert RoundUp
f64_convert_i64_unsigned_floor :: Word64 -> Double
f64_convert_i64_unsigned_floor = convert RoundDown
f64_convert_i64_signed_floor :: Int64 -> Double
f64_convert_i64_signed_floor = convert RoundDown
f64_convert_i64_unsigned_truncate :: Word64 -> Double
f64_convert_i64_unsigned_truncate = convert Truncate
f64_convert_i64_signed_truncate :: Int64 -> Double
f64_convert_i64_signed_truncate = convert Truncate
---
-- https://stackoverflow.com/questions/28949774/what-is-0-0-by-ieee-floating-point-standard
sanitizeZeroForAddition
:: (Eq a1, Num a1, RealFloat a2)
=> (a2 -> a2 -> a1)
-> (a2 -> a2 -> a1)
sanitizeZeroForAddition
operator
a
b
= if (0 == x) && (sign_bit a * sign_bit b == -1)
then negate x
else x
where
x = operator a b
sign_bit y
= if 0 == y
then (boolToNumber $ isNegativeZero y)
else signum y
boolToNumber True = -1
boolToNumber False = 1
sanitizeZeroForSubtraction
:: (Eq a1, RealFloat a2, Num a1)
=> (a2 -> a2 -> a1)
-> (a2 -> a2 -> a1)
sanitizeZeroForSubtraction
operator
a
b
= if (0 == x) && (sign_bit a == sign_bit b)
then negate x
else x
where
x = operator a b
sign_bit y
= if 0 == y
then (boolToNumber $ isNegativeZero y)
else signum y
boolToNumber True = -1
boolToNumber False = 1
sanatizeNegtiveZeroDemotionPromotion
:: (RealFloat p, Fractional a)
=> (p -> a)
-> (p -> a)
sanatizeNegtiveZeroDemotionPromotion
conversion
input
= result
where
result
= if isNegativeZero input
then negate 0.0
else conversion input
sanatizeNaNDemotion
:: (Double -> Float)
-> (Double -> Float)
sanatizeNaNDemotion
conversion
input
= result
where
result
= if isNaN input
then double2Float input
else conversion input
sanatizeNaNPromotion
:: (Float -> Double)
-> (Float -> Double)
sanatizeNaNPromotion
conversion
input
= result
where
result
= if isNaN input
then float2Double input
else conversion input
---------
data RationalWithSpecials
= NegativInfinity
| Exisiting Rational
| PositiveInfinity
deriving (Show, Eq)
getExisiting :: Rational -> RationalWithSpecials -> Rational
getExisiting _ (Exisiting x) = x
getExisiting defaultValue _ = defaultValue
-- | left: exact result, right: wannebe result
liftComparator
:: Round
-> (Rational -> Rational -> Bool)
-> RationalWithSpecials
-> RationalWithSpecials
-> Bool
liftComparator _ sorted (Exisiting x) (Exisiting y) = x `sorted` y
liftComparator RoundDown _ PositiveInfinity PositiveInfinity = True
liftComparator RoundDown _ _ PositiveInfinity = False
liftComparator RoundDown _ _ NegativInfinity = True
liftComparator RoundUp _ PositiveInfinity PositiveInfinity = True
liftComparator RoundUp _ _ PositiveInfinity = True
liftComparator RoundUp _ _ NegativInfinity = False
liftComparator Truncate _ xx@(Exisiting x) PositiveInfinity = if 0 < x then liftComparator RoundDown undefined xx PositiveInfinity else undefined --error "off by sign"
liftComparator Truncate _ xx@(Exisiting x) NegativInfinity = if x < 0 then liftComparator RoundUp undefined xx NegativInfinity else undefined --error "off by sign"
liftComparator Truncate _ PositiveInfinity PositiveInfinity = True
liftComparator _roundDirection _sorted _x _y = undefined --error $ show ("liftComparator", roundDirection, x, y)
liftComparatorConvert
:: Round
-> (Rational -> Rational -> Bool)
-> RationalWithSpecials
-> RationalWithSpecials
-> Bool
liftComparatorConvert _ sorted (Exisiting x) (Exisiting y) = x `sorted` y
liftComparatorConvert RoundDown _sorted PositiveInfinity (Exisiting _) = False
liftComparatorConvert RoundDown _ NegativInfinity (Exisiting _) = True
liftComparatorConvert RoundUp _sorted PositiveInfinity (Exisiting _) = True
liftComparatorConvert RoundUp _sorted NegativInfinity (Exisiting _) = False
liftComparatorConvert Truncate _ PositiveInfinity xx@(Exisiting x)
= if 0 < x
then liftComparatorConvert RoundDown undefined PositiveInfinity xx
else undefined --off_by_sign
liftComparatorConvert Truncate _ NegativInfinity xx@(Exisiting x)
= if x < 0
then liftComparatorConvert RoundUp undefined NegativInfinity xx
else undefined --off_by_sign
liftComparatorConvert _ _sorted PositiveInfinity PositiveInfinity = True
liftComparatorConvert _ _sorted NegativInfinity NegativInfinity = True
liftComparatorConvert _roundDirection _sorted _x _y = undefined --error $ show ("liftComparatorConvert", roundDirection, x, y)
liftUnary
:: (Rational -> Rational)
-> (RationalWithSpecials -> RationalWithSpecials)
liftUnary f (Exisiting x) = Exisiting $ f x
liftUnary _ PositiveInfinity = PositiveInfinity
liftUnary _ NegativInfinity = PositiveInfinity
class RationalWithSpecialsLike a where
exactToRational :: {-HasCallStack =>-} a -> RationalWithSpecials
class
(Show a, Num a, RealFloat a, RationalWithSpecialsLike a)
=> IEEE_RoundingFuckUp_Correctable a
where
exactFromRational :: {-HasCallStack =>-} RationalWithSpecials -> a
nudgeDown :: a -> a
nudgeUp :: a -> a
nudgeToZero :: a -> a
nudgeToZero j
= if j < 0
then nudgeUp j
else if 0 < j
then nudgeDown j
else undefined -- error "what now?"
nudgeAwayFromZero :: a -> a
nudgeAwayFromZero j
= if j < 0
then nudgeDown j
else if 0 < j
then nudgeUp j
else undefined -- error "what now?"
instance RationalWithSpecialsLike Word32 where
exactToRational = Exisiting . toRational
instance RationalWithSpecialsLike Word64 where
exactToRational = Exisiting . toRational
instance RationalWithSpecialsLike Int32 where
exactToRational = Exisiting . toRational
instance RationalWithSpecialsLike Int64 where
exactToRational = Exisiting . toRational
withFrozenCallStack :: a -> a
withFrozenCallStack = id
assertNote :: a -> b -> c -> c
assertNote _note _check = id
instance IEEE_RoundingFuckUp_Correctable Double where
exactFromRational = withFrozenCallStack hopefullyExactFromRational
nudgeDown = f64_predecessorIEEE
nudgeUp = f64_successorIEEE
instance RationalWithSpecialsLike Double where
exactToRational = withFrozenCallStack hopeFullyExactToRational
instance IEEE_RoundingFuckUp_Correctable Float where
exactFromRational = withFrozenCallStack hopefullyExactFromRational
nudgeDown = f32_predecessorIEEE
nudgeUp = f32_successorIEEE
instance RationalWithSpecialsLike Float where
exactToRational = withFrozenCallStack hopeFullyExactToRational
hopeFullyExactToRational
:: ({-HasCallStack,-} RealFloat a, Show a)
=> a
-> RationalWithSpecials
hopeFullyExactToRational
value
= id
$ withFrozenCallStack
$ assertNote note check
$ result
where
result
= if isInfinite value
then case signum value of
-1 -> NegativInfinity
1 -> PositiveInfinity
_ -> undefined -- error "infinit has sign 0"
else if not $ isNaN value
then Exisiting $ toRational value
else undefined --error $ "cannot convert to rational: " ++ show value
note = "this double cannot be represented as Rational: " ++ show value
check
= (not $ isNaN $ value)
&& (not $ isNegativeZero $ value)
hopefullyExactFromRational
:: ({-HasCallStack,-} RealFloat a)
=> RationalWithSpecials
-> a
hopefullyExactFromRational (Exisiting rational) = fromRational rational
hopefullyExactFromRational (PositiveInfinity) = 1/0
hopefullyExactFromRational (NegativInfinity) = -1/0
------------
data Round = RoundDown | RoundUp | Truncate
deriving Show
roundingOrder_binaryOperator
:: Round
-> (Rational -> Rational -> Bool)
roundingOrder_binaryOperator RoundDown = (>=)
roundingOrder_binaryOperator RoundUp = (<=)
roundingOrder_binaryOperator Truncate = \exact_result allegedly_rounded_result
-> if exact_result < 0 && allegedly_rounded_result < 0
then (roundingOrder_binaryOperator RoundUp) exact_result allegedly_rounded_result
else if 0 < exact_result && 0 < allegedly_rounded_result
then roundingOrder_binaryOperator RoundDown exact_result allegedly_rounded_result
else if 0 == exact_result && 0 == allegedly_rounded_result
then True
else if 0 == allegedly_rounded_result
then True -- zero is allways rounded ok
else undefined -- error $ "roundingOrder_binaryOperator: unconsiderd case: " ++ show (exact_result, allegedly_rounded_result)
-- exact_result `sorted` rounded_result
-- -2.0 `sorted` -1.0
-- 1.0 `sorted` 2.0
roundingOrder_UnaryOperator
:: Round
-> (Rational -> Rational -> Bool)
roundingOrder_UnaryOperator RoundDown = (<=)
roundingOrder_UnaryOperator RoundUp = (>=)
roundingOrder_UnaryOperator Truncate = \allegedly_rounded_result exact_result
-> if allegedly_rounded_result < 0 && exact_result < 0
then (roundingOrder_UnaryOperator RoundUp) allegedly_rounded_result exact_result
else if 0 < allegedly_rounded_result && 0 < exact_result
then roundingOrder_UnaryOperator RoundDown allegedly_rounded_result exact_result
else if 0 == allegedly_rounded_result && 0 == exact_result
then True
else if 0 == allegedly_rounded_result
then True -- zero is allways rounded ok
else undefined -- error $ "roundingOrder_binaryOperator: unconsiderd case: " ++ show (exact_result, allegedly_rounded_result)
chooseNudging
:: IEEE_RoundingFuckUp_Correctable a
=> Round
-> a
-> a
chooseNudging RoundDown = nudgeDown
chooseNudging RoundUp = nudgeUp
chooseNudging Truncate = nudgeToZero
chooseTheTheOtherNudging
:: IEEE_RoundingFuckUp_Correctable a
=> Round
-> a
-> a
chooseTheTheOtherNudging RoundDown = nudgeUp
chooseTheTheOtherNudging RoundUp = nudgeDown
chooseTheTheOtherNudging Truncate = nudgeAwayFromZero
------------
binary_operator_rounded
:: Round
-> (forall a. (Num a, Fractional a) => a -> a -> a)
-> ((IEEE_RoundingFuckUp_Correctable b) => b -> b -> b)
binary_operator_rounded
upOrDown
operator
left
right
= id
$ result
where
sorted
= liftComparator upOrDown
$ roundingOrder_binaryOperator upOrDown
nudge = chooseNudging upOrDown
the_other_nudge = chooseTheTheOtherNudging upOrDown
badArguments = undefined --error $ show ("case not supported", left, right)
left_rational = getExisiting badArguments $ exactToRational left
right_rational = getExisiting badArguments $ exactToRational right
exact_result :: Rational
exact_result = left_rational `operator` right_rational
somehow_rounded_result = left `operator` right
somehow_rounded_result_as_rational
= exactToRational
$ somehow_rounded_result
was_already_rounded_correctly
= (Exisiting exact_result) `sorted` somehow_rounded_result_as_rational
correctly_rounded_result
= if not (isNaN somehow_rounded_result || 0 == right)
then if was_already_rounded_correctly
then somehow_rounded_result
else nudge somehow_rounded_result
else somehow_rounded_result
note
= "new nudging did not work: "
++ show
( upOrDown
, left
, right
)
check
= is_special_case
||
(True
&& (Exisiting exact_result)
`sorted`
(exactToRational correctly_rounded_result)
&& (exactToRational $ the_other_nudge correctly_rounded_result)
`strict_sorted`
(Exisiting exact_result)
)
strict_sorted x y = (x /= y) && (x `sorted` y)
is_special_case = False
|| isNaN left
|| isNaN right
|| isNegativeZero left
|| isNegativeZero right
|| isInfinite left
|| isInfinite right
result
= if not (isNaN left || isNaN right)
then normal_case
else somehow_rounded_result
normal_case
= if not(isInfinite left || isInfinite right)
then propernumbercase
else somehow_rounded_result
propernumbercase = assertNote note check correctly_rounded_result
------------
squareRoot
:: (Floating a, IEEE_RoundingFuckUp_Correctable a)
=> Round
-> a
-> a
squareRoot
upOrDown
preImage
= result
where
sorted
= liftComparator upOrDown
$ roundingOrder_UnaryOperator upOrDown
-- TODO refactor into one function
nudge = chooseNudging upOrDown
the_other_nudge = chooseTheTheOtherNudging upOrDown
result = normal_case
somehowRoundedResult = sqrt preImage
preImage_Rational = exactToRational preImage
somehowRoundedResult_preImage :: RationalWithSpecials
somehowRoundedResult_preImage
= id
$ liftUnary (\x -> x*x)
$ exactToRational
$ somehowRoundedResult
-- this works because squareroot is monotonic
wasRoundedCorrectly
= somehowRoundedResult_preImage
`sorted` preImage_Rational
correctly_rounded_result
= if (not $ isNaN somehowRoundedResult) || (isNegativeZero preImage)
then
if wasRoundedCorrectly
then somehowRoundedResult
else nudge somehowRoundedResult
else
somehowRoundedResult
normal_case = assertNote note check correctly_rounded_result
note = "something went wrong"
check
= True
&& we_rounded_correctly
&& we_rounded_tightly
we_rounded_correctly = ((liftUnary (\x->x*x) $ exactToRational correctly_rounded_result) `sorted` (exactToRational preImage))
we_rounded_tightly = (exactToRational preImage) `strict_sorted` (liftUnary (\x->x*x) $ exactToRational (the_other_nudge correctly_rounded_result))
where
strict_sorted x y = (x /= y) && (x `sorted` y)
convert
::
( Show a
, Real a
, RationalWithSpecialsLike a
, IEEE_RoundingFuckUp_Correctable b
)
=> Round
-> a
-> b
convert
upOrDown
integral
= id
$ result
where
result = normal_case
sorted
= liftComparatorConvert upOrDown
$ roundingOrder_UnaryOperator upOrDown
nudge = chooseNudging upOrDown
the_other_nudge = chooseTheTheOtherNudging upOrDown
argument_Rational = exactToRational $ integral
somehowRoundedResult = exactFromRational $ argument_Rational
somehow_rounded_result_Rational
= exactToRational somehowRoundedResult
wasRoundedCorrectly
= somehow_rounded_result_Rational
`sorted` argument_Rational
correctly_rounded_result
= if wasRoundedCorrectly
then somehowRoundedResult
else nudge somehowRoundedResult
normal_case = assertNote note check correctly_rounded_result
note = "something went wrong"
check
= True
&& we_rounded_correctly
&& we_rounded_tightly
we_rounded_correctly = (Exisiting $ toRational correctly_rounded_result) `sorted` (Exisiting $ toRational integral)
we_rounded_tightly = (Exisiting $ toRational integral) `strict_sorted` (exactToRational $ the_other_nudge correctly_rounded_result)
where
strict_sorted = liftComparator upOrDown $ \x y -> (x /= y) && (Exisiting x `sorted` Exisiting y)