module Numeric.VariablePrecision.Algorithms
( recodeFloat
, viaDouble
, (=~=)
, genericRecip
, genericSqrt
, genericExp
, genericLog
, genericLog'
, genericLog2
, genericLog''
, genericPi
, genericSin
, genericPositiveZero
, genericNegativeZero
, genericPositiveInfinity
, genericNegativeInfinity
, genericNotANumber
, sameSign
) where
import Data.Bits (bit, shiftR)
import Data.List (foldl')
genericPositiveZero, genericNegativeZero, genericPositiveInfinity, genericNegativeInfinity, genericNotANumber :: RealFloat a => a
genericPositiveZero = 0
genericNegativeZero = 0
genericPositiveInfinity = result
where
result = encodeFloat m e
m = bit (floatDigits (undefined `asTypeOf` result))
e = snd (floatRange (undefined `asTypeOf` result))
genericNegativeInfinity = result
where
result = encodeFloat (negate m) e
m = bit (floatDigits (undefined `asTypeOf` result))
e = snd (floatRange (undefined `asTypeOf` result))
genericNotANumber = genericPositiveInfinity + genericNegativeInfinity
recodeFloat :: (RealFloat a, RealFloat b) => a -> b
recodeFloat !x
| isNaN x = genericNotANumber
| isInfinite x && x > 0 = genericPositiveInfinity
| isInfinite x && x < 0 = genericNegativeInfinity
| isNegativeZero x = genericNegativeZero
| x == 0 = genericPositiveZero
| otherwise = uncurry encodeFloat (decodeFloat x)
sameSign :: (Ord a, Num a) => a -> a -> Bool
sameSign a b = compare 0 a == compare 0 b
(=~=) :: RealFloat a => a -> a -> Int -> Bool
(=~=) !x !y !s
| x == y = True
| isNaN x && isNaN y = True
| isNaN x || isNaN y = False
| isInfinite x || isInfinite y = False
| not (sameSign a b) = False
| otherwise = abs (e f) <= s && abs (x y) <= encodeFloat 1 (s + (e `max` f))
where
(a, e) = decodeFloat x
(b, f) = decodeFloat y
genericRecip :: RealFloat a => Int -> a -> a
genericRecip accuracy y = recip' y
where
recip' f0
| isNaN f0 = f0
| isInfinite f0 && f0 > 0 = genericPositiveZero
| isInfinite f0 && f0 < 0 = genericNegativeZero
| isNegativeZero f0 = genericNegativeInfinity
| f0 == 0 = genericPositiveInfinity
| f0 < 0 = negate . recip' . negate $ f0
| otherwise = scaleFloat sh (go d s0 x0)
where
x0 = k48 k32 * d
d = significand f0
sh = exponent d exponent f0
go !d !s !x
| (x =~= x') accuracy = x'
| s == 0 = x'
| otherwise = go d (s 1) x'
where
x' = scaleFloat 1 x d * x * x
p = floatDigits (undefined `asTypeOf` y)
s0 = ceiling (logBase 2 (fromIntegral (p + 1) / logBase 2 17) :: Double) :: Int
k48 = recodeFloat (48/17 :: Double)
k32 = recodeFloat (32/17 :: Double)
genericSqrt :: RealFloat a => Int -> a -> a
genericSqrt accuracy f0
| f0 < 0 = genericNotANumber
| f0 == 0 = f0
| isNaN f0 = f0
| isInfinite f0 = f0
| otherwise = go (viaDouble sqrt f)
where
e = exponent f0
d = if even e then 2 else 1
s = e d
f = scaleFloat (negate s) f0
go !r =
let r' = scaleFloat (1) (r + f / r)
in if (r =~= r') accuracy then scaleFloat (s `shiftR` 1) r' else go r'
genericExp :: RealFloat a => Int -> a -> a
genericExp accuracy x
| isNaN x = x
| isInfinite x && x < 0 = 0
| isInfinite x = x
| x == 0 = 1
| x < 0 = recip . genericExp accuracy . negate $ x
| otherwise = go 0 1 1
where
go !s !xnnf !n
| (s =~= s') accuracy = s'
| otherwise = go s' (xnnf * x / fromIntegral n) (n + 1 :: Int)
where
s' = s + xnnf
genericLog :: RealFloat a => Int -> a -> a
genericLog accuracy = genericLog' accuracy (genericLog2 accuracy)
genericLog2 :: RealFloat a => Int -> a
genericLog2 accuracy = negate (genericLog'' accuracy 0.5)
genericLog' :: RealFloat a => Int -> a -> a -> a
genericLog' accuracy ln2 x
| isNaN x = x
| x == 0 = genericNegativeInfinity
| x < 0 = genericNotANumber
| isInfinite x = x
| otherwise = mln2 + genericLog'' accuracy s
where
m = exponent x
s = significand x
mln2
| m == 0 = 0
| otherwise = fromIntegral m * ln2
genericLog'' :: RealFloat a => Int -> a -> a
genericLog'' accuracy x = result
where
result = go (1) 1 (encodeFloat 1 m) 0 1 (scaleFloat m x) 0
m2 = accuracy floatDigits (undefined `asTypeOf` result)
m = m2 `shiftR` 1
small y = y == 0 || exponent y <= m2
go !n !a !b !s !c !d !t
| small ds && small dt = recip (1 s') recip (1 t')
| otherwise = go n' a' b' s' c' d' t'
where
a' = scaleFloat (1) (a + b)
c' = scaleFloat (1) (c + d)
b' = sqrt (a * b)
d' = sqrt (c * d)
ds = scaleFloat n (a * a b * b)
dt = scaleFloat n (c * c d * d)
t' = t + dt
s' = s + ds
n' = n + 1
genericPi :: RealFloat a => Int -> a
genericPi accuracy = result
where
sqr x = x * x
result = go 1 (sqrt 0.5) 0.25 0 1
go !a !b !t !k !p
| (p =~= p') accuracy = p'
| otherwise = go a' b' t' k' p'
where
a' = scaleFloat (1) (a + b)
b' = sqrt (a * b)
t' = t scaleFloat k (sqr (a' a))
k' = k + 1
p' = scaleFloat (2) (sqr (a + b) / t)
genericSin :: RealFloat a => Int -> a -> a -> a
genericSin accuracy pi x0 = reduced taylor x0
where
sqr y = y * y
t :: Double
t = fromIntegral (floatDigits x0 accuracy)
k :: Int
k = round (t / 3)
three = 3 ^ k
up 0 !y = y
up n !y = up (n 1) (y * (3 scaleFloat 2 (sqr y)))
reduced f y
| y == 0 = y
| y < 0 = (negate . reduced f . negate) y
| y <= pi / 4 = (up k . f . (/ three)) y
| y <= pi / 2 = (sqrt . (1 ) . sqr . reduced f . (pi / 2 )) y
| y <= pi = (reduced f . (pi )) y
| y <= pi * 2 = (negate . reduced f . (pi * 2 )) y
| otherwise = (reduced f . subtract (pi * 2)) y
taylor y = (sum' . reverse . takeWhile ((> threshold) . abs) . go 1) y
where
threshold = scaleFloat (negate (floatDigits y * 2)) y
x2 = sqr y
go !n !xnf = xnf : go n' xnf'
where
n' = n + 1
xnf' = negate (x2 * xnf / fromInt (2 * n + 1))
fromInt :: Num a => Int -> a
fromInt = fromIntegral
sum' :: Num a => [a] -> a
sum' = foldl' (+) 0
viaDouble :: (RealFloat a, RealFloat b) => (Double -> Double) -> a -> b
viaDouble f = recodeFloat . f . recodeFloat