{-# OPTIONS_GHC -fno-warn-orphans #-}
module Data.BigFloating
  ( piChudnovsky
  , sqr
  , nthRoot
  )
where

import           Data.BigDecimal
import           Data.List  (find)
import           Data.Maybe (fromMaybe, fromJust)
import           GHC.Real   ((%), Ratio ((:%)))
import           GHC.Natural

-- I'm giving some implementation ideas for approximisations for functions on transcendental numbers.

-- The rest is left as an exercise to the interested reader ;-)

instance Floating BigDecimal where
    pi :: BigDecimal
pi    = RoundingAdvice -> BigDecimal
piChudnovsky RoundingAdvice
defaultRounding
    exp :: BigDecimal -> BigDecimal
exp   = forall a. HasCallStack => a
undefined -- e^x

    log :: BigDecimal -> BigDecimal
log   = forall a. HasCallStack => a
undefined
    sin :: BigDecimal -> BigDecimal
sin   = forall a. HasCallStack => a
undefined
    cos :: BigDecimal -> BigDecimal
cos   = forall a. HasCallStack => a
undefined
    asin :: BigDecimal -> BigDecimal
asin  = forall a. HasCallStack => a
undefined
    acos :: BigDecimal -> BigDecimal
acos  = forall a. HasCallStack => a
undefined
    atan :: BigDecimal -> BigDecimal
atan  = forall a. HasCallStack => a
undefined
    sinh :: BigDecimal -> BigDecimal
sinh  = forall a. HasCallStack => a
undefined
    cosh :: BigDecimal -> BigDecimal
cosh  = forall a. HasCallStack => a
undefined
    asinh :: BigDecimal -> BigDecimal
asinh = forall a. HasCallStack => a
undefined
    acosh :: BigDecimal -> BigDecimal
acosh = forall a. HasCallStack => a
undefined
    atanh :: BigDecimal -> BigDecimal
atanh = forall a. HasCallStack => a
undefined

-- not required for minimal implementation

    sqrt :: BigDecimal -> BigDecimal
sqrt BigDecimal
x = BigDecimal -> RoundingAdvice -> BigDecimal
sqr BigDecimal
x RoundingAdvice
defaultRounding
    BigDecimal
x ** :: BigDecimal -> BigDecimal -> BigDecimal
** BigDecimal
y = BigDecimal -> Natural -> RoundingAdvice -> BigDecimal
nthRoot (BigDecimal
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
b) (forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n) RoundingAdvice
defaultRounding
                where
                  (Integer
b :% Integer
n) = forall a. Real a => a -> Rational
toRational BigDecimal
y

defaultRounding :: RoundingAdvice
defaultRounding :: RoundingAdvice
defaultRounding = (RoundingMode
DOWN, forall a. a -> Maybe a
Just Natural
100)

-- | computes the square root of any non-negative BigDecimal, rounding and precision defined by RoundingAdvice.

--   We are using Newton's algorithm.

sqr :: BigDecimal -> RoundingAdvice -> BigDecimal
sqr :: BigDecimal -> RoundingAdvice -> BigDecimal
sqr BigDecimal
x RoundingAdvice
mc
  | BigDecimal
x forall a. Ord a => a -> a -> Bool
<  BigDecimal
0    = forall a. HasCallStack => [Char] -> a
error [Char]
"can't determine the square root of negative numbers"
  | BigDecimal
x forall a. Eq a => a -> a -> Bool
== BigDecimal
0    = BigDecimal
0
  | Bool
otherwise = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a -> a
fromMaybe (forall a. HasCallStack => [Char] -> a
error [Char]
"did not find a sqrt") forall a b. (a -> b) -> a -> b
$ BigDecimal
-> BigDecimal -> RoundingAdvice -> Maybe (BigDecimal, Natural)
refine BigDecimal
x BigDecimal
1 RoundingAdvice
mc
      where
        refine :: BigDecimal
-> BigDecimal -> RoundingAdvice -> Maybe (BigDecimal, Natural)
refine BigDecimal
_ BigDecimal
_ (RoundingMode
_, Maybe Natural
Nothing)           = forall a. HasCallStack => [Char] -> a
error [Char]
"can't produce square root with unlimited precision"
        refine BigDecimal
r BigDecimal
initial ra :: RoundingAdvice
ra@(RoundingMode
_, Just Natural
scl) = forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find (BigDecimal, Natural) -> Bool
withinPrecision forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate forall {b}. Num b => (BigDecimal, b) -> (BigDecimal, b)
nextGuess (BigDecimal
initial, Natural
0)
          where
            withinPrecision :: (BigDecimal, Natural) -> Bool
withinPrecision (BigDecimal
guess, Natural
count) = forall a. Num a => a -> a
abs (BigDecimal
guessforall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int) forall a. Num a => a -> a -> a
- BigDecimal
r) forall a. Ord a => a -> a -> Bool
< Integer -> Natural -> BigDecimal
BigDecimal Integer
10 Natural
scl Bool -> Bool -> Bool
|| Natural
count forall a. Ord a => a -> a -> Bool
> Natural
10 forall a. Num a => a -> a -> a
* Natural
scl forall a. Num a => a -> a -> a
* BigDecimal -> Natural
precision BigDecimal
r
            nextGuess :: (BigDecimal, b) -> (BigDecimal, b)
nextGuess (BigDecimal
guess, b
count) = (BigDecimal -> BigDecimal
nf forall a b. (a -> b) -> a -> b
$ (BigDecimal, BigDecimal) -> RoundingAdvice -> BigDecimal
divide (BigDecimal
guess forall a. Num a => a -> a -> a
+ (BigDecimal, BigDecimal) -> RoundingAdvice -> BigDecimal
divide (BigDecimal
r, BigDecimal
guess) RoundingAdvice
mc, BigDecimal
2) RoundingAdvice
ra, b
countforall a. Num a => a -> a -> a
+b
1)

nthRoot :: BigDecimal -> Natural -> RoundingAdvice -> BigDecimal
nthRoot :: BigDecimal -> Natural -> RoundingAdvice -> BigDecimal
nthRoot BigDecimal
x Natural
n mc :: RoundingAdvice
mc@(RoundingMode
rm, Maybe Natural
maybeScale)
  | BigDecimal
x forall a. Ord a => a -> a -> Bool
<  BigDecimal
0 Bool -> Bool -> Bool
&& forall a. Integral a => a -> Bool
even Natural
n   = forall a. HasCallStack => [Char] -> a
error [Char]
"can't determine even roots of negative numbers"
  | BigDecimal
x forall a. Ord a => a -> a -> Bool
<  BigDecimal
0 Bool -> Bool -> Bool
&& forall a. Integral a => a -> Bool
odd  Natural
n   = - BigDecimal -> Natural -> RoundingAdvice -> BigDecimal
nthRoot BigDecimal
x (-Natural
n) RoundingAdvice
mc
  | BigDecimal
x forall a. Eq a => a -> a -> Bool
== BigDecimal
0    = BigDecimal
0
  | Bool
otherwise = BigDecimal -> RoundingAdvice -> BigDecimal
roundBD (forall a b. (a, b) -> a
fst (forall a. a -> Maybe a -> a
fromMaybe (forall a. HasCallStack => [Char] -> a
error [Char]
"did not find a sqrt") forall a b. (a -> b) -> a -> b
$ BigDecimal
-> BigDecimal -> RoundingAdvice -> Maybe (BigDecimal, Natural)
refine BigDecimal
x BigDecimal
1 (RoundingMode
rm, forall a. a -> Maybe a
Just (Natural
sforall a. Num a => a -> a -> a
+Natural
4)))) RoundingAdvice
mc
      where
        s :: Natural
s  = forall a. HasCallStack => Maybe a -> a
fromJust Maybe Natural
maybeScale
        refine :: BigDecimal
-> BigDecimal -> RoundingAdvice -> Maybe (BigDecimal, Natural)
refine BigDecimal
_ BigDecimal
_ (RoundingMode
_, Maybe Natural
Nothing)             = forall a. HasCallStack => [Char] -> a
error [Char]
"can't produce nth root with unlimited precision"
        refine BigDecimal
r BigDecimal
initial ra :: RoundingAdvice
ra@(RoundingMode
_, Just Natural
scl) = forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find (BigDecimal, Natural) -> Bool
withinPrecision forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate forall {b}. Num b => (BigDecimal, b) -> (BigDecimal, b)
nextGuess (BigDecimal
initial, Natural
0)
          where
            withinPrecision :: (BigDecimal, Natural) -> Bool
withinPrecision (BigDecimal
guess, Natural
count) = forall a. Num a => a -> a
abs (BigDecimal
guessforall a b. (Num a, Integral b) => a -> b -> a
^Natural
n forall a. Num a => a -> a -> a
- BigDecimal
r) forall a. Ord a => a -> a -> Bool
< Integer -> Natural -> BigDecimal
BigDecimal (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ Natural
nforall a. Num a => a -> a -> a
*Natural
10) Natural
scl Bool -> Bool -> Bool
|| Natural
count forall a. Ord a => a -> a -> Bool
> Natural
10 forall a. Num a => a -> a -> a
* Natural
scl forall a. Num a => a -> a -> a
* BigDecimal -> Natural
precision BigDecimal
r
            nextGuess :: (BigDecimal, b) -> (BigDecimal, b)
nextGuess (BigDecimal
guess, b
count) =
              (BigDecimal -> BigDecimal
nf forall a b. (a -> b) -> a -> b
$ (BigDecimal, BigDecimal) -> RoundingAdvice -> BigDecimal
divide ((BigDecimal
guess forall a. Num a => a -> a -> a
* Integer -> Natural -> BigDecimal
BigDecimal (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ Natural
nforall a. Num a => a -> a -> a
-Natural
1) Natural
0) forall a. Num a => a -> a -> a
+ (BigDecimal, BigDecimal) -> RoundingAdvice -> BigDecimal
divide (BigDecimal
r, BigDecimal
guessforall a b. (Num a, Integral b) => a -> b -> a
^(Natural
nforall a. Num a => a -> a -> a
-Natural
1)) RoundingAdvice
ra, Integer -> Natural -> BigDecimal
BigDecimal (forall a b. (Integral a, Num b) => a -> b
fromIntegral Natural
n) Natural
0) RoundingAdvice
ra, b
countforall a. Num a => a -> a -> a
+b
1)


-- | Compute pi using rounding mode and scale of the specified RoundingAdvice

--   Sources: https://wiki.haskell.org/Integers_too_big_for_floats & https://github.com/eobermuhlner/big-math

piChudnovsky :: RoundingAdvice -> BigDecimal
piChudnovsky :: RoundingAdvice -> BigDecimal
piChudnovsky (RoundingMode
_, Maybe Natural
Nothing)           = forall a. HasCallStack => [Char] -> a
error [Char]
"can't compute pi with umlimited precision"
piChudnovsky mc :: RoundingAdvice
mc@(RoundingMode
rMode, Just Natural
scl) = (BigDecimal, BigDecimal) -> RoundingAdvice -> BigDecimal
divide (BigDecimal
1, BigDecimal
12 forall a. Num a => a -> a -> a
* (BigDecimal, BigDecimal) -> RoundingAdvice -> BigDecimal
divide (Rational -> RoundingAdvice -> BigDecimal
fromRatio Rational
s RoundingAdvice
mc,BigDecimal
f) RoundingAdvice
mc') RoundingAdvice
mc
    where
      mc' :: RoundingAdvice
mc'   = (RoundingMode
rMode, forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ Natural
scl forall a. Num a => a -> a -> a
+ Natural
3) -- increase precision to avoid propagation of rounding errors

      steps :: Natural
steps = Natural
1 forall a. Num a => a -> a -> a
+ forall a. Integral a => a -> a -> a
div Natural
scl  Natural
14         -- taken from github.com/eobermuhlner/big-math

      s :: Rational
s = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Integer -> Rational
chudnovsky (forall a b. (Integral a, Num b) => a -> b
fromIntegral Natural
n) | Natural
n <- [Natural
0..Natural
steps]] :: Rational
      f :: BigDecimal
f = BigDecimal -> RoundingAdvice -> BigDecimal
sqr (forall a. Num a => Integer -> a
fromInteger Integer
cforall a b. (Num a, Integral b) => a -> b -> a
^(Int
3::Int)) RoundingAdvice
mc      -- Common factor in the sum


      -- k-th term of the Chudnovsky series

      chudnovsky :: Integer -> Rational
      chudnovsky :: Integer -> Rational
chudnovsky Integer
k
          | forall a. Integral a => a -> Bool
even Integer
k    =  Rational
quotient
          | Bool
otherwise = -Rational
quotient
          where
            quotient :: Rational
quotient = Integer
num forall a. Integral a => a -> a -> Ratio a
% Integer
den
            num :: Integer
num  = Integer -> Integer -> Integer
facDiv (Integer
6 forall a. Num a => a -> a -> a
* Integer
k) (Integer
3 forall a. Num a => a -> a -> a
* Integer
k) forall a. Num a => a -> a -> a
* (Integer
a forall a. Num a => a -> a -> a
+ Integer
b forall a. Num a => a -> a -> a
* Integer
k)
            den :: Integer
den  = forall a. (Enum a, Num a) => a -> a
fac Integer
k forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3::Int) forall a. Num a => a -> a -> a
* (Integer
c forall a b. (Num a, Integral b) => a -> b -> a
^ (Integer
3 forall a. Num a => a -> a -> a
* Integer
k))

      -- Compute n!

      fac :: (Enum a, Num a) => a -> a
      fac :: forall a. (Enum a, Num a) => a -> a
fac a
n = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [a
1..a
n]

      -- Compute n! / m! efficiently

      facDiv :: Integer -> Integer -> Integer
      facDiv :: Integer -> Integer -> Integer
facDiv Integer
n Integer
m
          | Integer
n forall a. Ord a => a -> a -> Bool
> Integer
m     = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
n, Integer
n forall a. Num a => a -> a -> a
- Integer
1 .. Integer
m forall a. Num a => a -> a -> a
+ Integer
1]
          | Integer
n forall a. Eq a => a -> a -> Bool
== Integer
m    = Integer
1
          | Bool
otherwise = Integer -> Integer -> Integer
facDiv Integer
m Integer
n

      a :: Integer
a = Integer
13591409
      b :: Integer
b = Integer
545140134
      c :: Integer
c = Integer
640320