{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE Safe #-}
module Data.Connection.Float (
f32w08,
f32w16,
f32i08,
f32i16,
f32f32,
ulp32,
near32,
shift32,
f64f64,
f64w08,
f64w16,
f64w32,
f64i08,
f64i16,
f64i32,
f64f32,
ulp64,
near64,
shift64,
until,
) where
import safe Data.Bool
import safe Data.Connection.Conn hiding (ceiling, floor)
import safe Data.Int
import safe Data.Order
import safe Data.Order.Syntax
import safe Data.Word
import safe GHC.Float as F
import safe Prelude hiding (Eq (..), Ord (..), until)
import safe qualified Prelude as P
f32f32 :: Conn k (Float, Float) Float
f32f32 = fxxfxx
f32w08 :: Conn k Float (Extended Word8)
f32w08 = fxxext
f32w16 :: Conn k Float (Extended Word16)
f32w16 = fxxext
f32i08 :: Conn k Float (Extended Int8)
f32i08 = fxxext
f32i16 :: Conn k Float (Extended Int16)
f32i16 = fxxext
ulp32 :: Float -> Float -> Maybe (Ordering, Word32)
ulp32 x y = fmap f $ pcompare x y
where
x' = floatInt32 x
y' = floatInt32 y
f LT = (LT, fromIntegral . abs $ y' - x')
f EQ = (EQ, 0)
f GT = (GT, fromIntegral . abs $ x' - y')
near32 :: Word32 -> Float -> Float -> Bool
near32 tol x y = maybe False ((<= tol) . snd) $ ulp32 x y
shift32 :: Int32 -> Float -> Float
shift32 n x =
if isNaN x == True
then case signum n of
-1 -> -1 / 0
1 -> 1 / 0
_ -> 0 / 0
else int32Float . clamp32 . (+ n) . floatInt32 $ x
f64f64 :: Conn k (Double, Double) Double
f64f64 = fxxfxx
f64w08 :: Conn k Double (Extended Word8)
f64w08 = fxxext
f64w16 :: Conn k Double (Extended Word16)
f64w16 = fxxext
f64w32 :: Conn k Double (Extended Word32)
f64w32 = fxxext
f64i08 :: Conn k Double (Extended Int8)
f64i08 = fxxext
f64i16 :: Conn k Double (Extended Int16)
f64i16 = fxxext
f64i32 :: Conn k Double (Extended Int32)
f64i32 = fxxext
f64f32 :: Conn k Double Float
f64f32 = Conn f g h
where
f x =
let est = double2Float x
in if g est >~ x
then est
else ascend32 est g x
g = float2Double
h x =
let est = double2Float x
in if g est <~ x
then est
else descend32 est g x
ascend32 z g1 y = until (\x -> g1 x >~ y) (<~) (shift32 1) z
descend32 z h1 x = until (\y -> h1 y <~ x) (>~) (shift32 (-1)) z
{-# INLINE f64f32 #-}
ulp64 :: Double -> Double -> Maybe (Ordering, Word64)
ulp64 x y = fmap f $ pcompare x y
where
x' = doubleInt64 x
y' = doubleInt64 y
f LT = (LT, fromIntegral . abs $ y' - x')
f EQ = (EQ, 0)
f GT = (GT, fromIntegral . abs $ x' - y')
near64 :: Word64 -> Double -> Double -> Bool
near64 tol x y = maybe False ((<= tol) . snd) $ ulp64 x y
shift64 :: Int64 -> Double -> Double
shift64 n x =
if isNaN x == True
then case signum n of
-1 -> -1 / 0
1 -> 1 / 0
_ -> 0 / 0
else int64Double . clamp64 . (+ n) . doubleInt64 $ x
{-# INLINE until #-}
until :: (a -> Bool) -> (a -> a -> Bool) -> (a -> a) -> a -> a
until pre rel f seed = go seed
where
go x
| x' `rel` x = x
| pre x = x
| otherwise = go x'
where
x' = f x
signed32 :: Word32 -> Int32
signed32 x
| x < 0x80000000 = fromIntegral x
| otherwise = fromIntegral (maxBound - (x - 0x80000000))
signed64 :: Word64 -> Int64
signed64 x
| x < 0x8000000000000000 = fromIntegral x
| otherwise = fromIntegral (maxBound - (x - 0x8000000000000000))
unsigned32 :: Int32 -> Word32
unsigned32 x
| x >= 0 = fromIntegral x
| otherwise = 0x80000000 + (maxBound - (fromIntegral x))
unsigned64 :: Int64 -> Word64
unsigned64 x
| x >= 0 = fromIntegral x
| otherwise = 0x8000000000000000 + (maxBound - (fromIntegral x))
int32Float :: Int32 -> Float
int32Float = castWord32ToFloat . unsigned32
floatInt32 :: Float -> Int32
floatInt32 = signed32 . (+ 0) . castFloatToWord32
int64Double :: Int64 -> Double
int64Double = castWord64ToDouble . unsigned64
doubleInt64 :: Double -> Int64
doubleInt64 = signed64 . (+ 0) . castDoubleToWord64
clamp32 :: Int32 -> Int32
clamp32 = P.max (-2139095041) . P.min 2139095040
clamp64 :: Int64 -> Int64
clamp64 = P.max (-9218868437227405313) . P.min 9218868437227405312
fxxfxx :: (Total a, Fractional a) => Conn k (a, a) a
fxxfxx = Conn f g h
where
f (x, y) = maybe (1 / 0) (bool y x . (>= EQ)) (pcompare x y)
g x = (x, x)
h (x, y) = maybe (-1 / 0) (bool y x . (<= EQ)) (pcompare x y)
fxxext :: forall a b k. (RealFrac a, Preorder a, Bounded b, Integral b) => Conn k a (Extended b)
fxxext = Conn f g h
where
f = extend (~~ -1 / 0) (\x -> x ~~ 0 / 0 || x > high) $ \x -> if x < low then minBound else ceiling x
g = extended (-1 / 0) (1 / 0) fromIntegral
h = extend (\x -> x ~~ 0 / 0 || x < low) (~~ 1 / 0) $ \x -> if x > high then maxBound else floor x
low = fromIntegral $ minBound @b
high = fromIntegral $ maxBound @b
{-# INLINE fxxext #-}