3

I am trying to write a program in Haskell that returns 'e' (Euler's number) to a given decimal place. Here is my code so far:

factorial 0 = 1
factorial n = n * factorial (n - 1)

calculateE a 
    | a == 0 = 1 
    | otherwise = nextLevel
    where nextLevel = (1 / (factorial a)) + calculateE (a-1)

Whenever I call calculateE I only get back 16 decimal places. Is this a limitation of Haskell/My computer? Is there a way to get back any number of decimal places?

Will Ness
  • 70,110
  • 9
  • 98
  • 181
Jack Bodine
  • 105
  • 7
  • 2
    You can use a type like `Scientific` instead: https://hackage.haskell.org/package/scientific-0.3.6.2/docs/Data-Scientific.html – Willem Van Onsem Sep 13 '20 at 20:00
  • 2
    This is not really Haskell-specific: most programming languages provide as primitives floating point types like `Float` and `Double`, which can store only a fixed number of significant digits. If you need something more, you have to use some library like the one suggested by Willem above. – chi Sep 13 '20 at 20:12
  • 1
    @WillemVanOnsem Actually, after reading some docs, I'm unsure if `Scientific` can be used here. Apparently, it stores all the digits eagerly, and fails on things like `1/3` since it can not store all the digits. I'm not sure if this will be a issue or not here. That `1 / factorial a` could easily fail. – chi Sep 13 '20 at 20:18
  • @chi: hmmm... There was some package if I recall correctly that could calculate up to an arbitrary amount of numbers. Although I can not immediately find it :s – Willem Van Onsem Sep 13 '20 at 20:44
  • IOW `calculateE a = sum [ 1 / product [1..n] | n <- [0..a]]`. then call it as `calculateE (100%1)`, as the answer suggests. – Will Ness Sep 16 '20 at 15:35

1 Answers1

8

This code already works to arbitrary precision. You just need to use an arbitrary precision type and not the standard Float/Double. Haskell's standard library has Rational for this purpose, which represents rational numbers as pairs of integers.

ghci> calculateE 100 :: Rational
4299778907798767752801199122242037634663518280784714275131782813346597523870956720660008227544949996496057758175050906671347686438130409774741771022426508339 % 1581800261761765299689817607733333906622304546853925787603270574495213559207286705236295999595873191292435557980122436580528562896896000000000000000000000000

The issue now is getting a sequence of digits out of it. I'm not aware of anything in the standard library that does it, so here's a stupid simple (might still be buggy!) implementation:

import Data.List(unfoldr)
import Data.List.NonEmpty(NonEmpty((:|)))
import Data.Ratio

-- first element is integral part (+ sign), rest are positive and < 10 and are digits
-- after the decimal point (for negative numbers, these digits should be seen as having negative value)
longDivision :: Integral a => Ratio a -> NonEmpty a
longDivision x = hi :| unfoldr go (abs lo)
  where (hi, lo) = numerator x `quotRem` denominator x
        go 0 = Nothing
        go lo = Just $ (lo * 10) `quotRem` denominator x

printDigits :: Show a => NonEmpty a -> String
printDigits (x :| xs) = show x ++ "." ++ concatMap show xs

So

ghci> take 100 $ printDigits $ longDivision $ calculateE 100
"2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642"

This approximation actually seems to be good to ~160 digits after the decimal point.

HTNW
  • 27,182
  • 1
  • 32
  • 60