Posts Tagged ‘integral’

Programming Praxis – Two Integrals

January 11, 2011

In today’s Programming Praxis exercise, our task is to write functions to calculate the exponential and logarithmic integrals. Let’s get started, shall we?

The exponential integral is just executable math, with the only addition being a limit on the infinite sum.

ei :: Double -> Double
ei x = 0.5772156649015328606065 + log x +
       sum (takeWhile (> 1e-17) [x**k / k / product [1..k] | k <- [1..]])

The two logarithmic integrals can be expressed in terms of the exponential one, so they become trivial to implement.

li :: Double -> Double
li = ei . log

liOffset :: Double -> Double
liOffset x = li x - li 2

A test to see if everything works properly:

main :: IO ()
main = print . round $ liOffset 1e21

We get a result of 21127269486616088576, which starts to deviate from the mathematically correct solution at the 15th digit; close enough.