{-# 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
instance Floating BigDecimal where
pi :: BigDecimal
pi = RoundingAdvice -> BigDecimal
piChudnovsky RoundingAdvice
defaultRounding
exp :: BigDecimal -> BigDecimal
exp = forall a. HasCallStack => a
undefined
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
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)
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)
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)
steps :: Natural
steps = Natural
1 forall a. Num a => a -> a -> a
+ forall a. Integral a => a -> a -> a
div Natural
scl Natural
14
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
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))
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]
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