Two Integrals

January 11, 2011

Our functions are straight forward translations of the mathematical formulas, summing the infinite series until the individual terms are less than the desired tolerance:

(define (exp-integral x)
  (let* ((gamma 0.5772156649015328606065))
    (let loop ((k 1) (fact 1) (z x) (sum (+ gamma (log x))))
      (let ((term (/ z k fact)))
        (if (< term 1e-17) sum
          (loop (+ k 1) (* fact (+ k 1)) (* z x) (+ sum term)))))))

(define (log-integral x)
  (let* ((gamma 0.5772156649015328606065)
         (logx (log x)) (sum (+ gamma (log logx))))
    (let loop ((k 1) (fact 1) (logx^k logx) (sum sum))
      (let* ((term (/ logx^k fact k)) (sum (+ sum term))
             (newk (+ k 1)) (newfact (* fact newk)))
        (if (< term 1e-17) sum
          (loop newk newfact (* logx^k logx) sum))))))

(define (offset-log-integral x)
  (- (log-integral x) 1.04516378011749278))

Note that we lose some precision when x is large; the proper value of Li(1021) is 21127269486616126182, but we compute 21127269486616100864. We computed Li(1021) = 21127269486616129536 using adaptive quadrature in a previous exercise.

You may be curious why we are interested in a function that approximates the prime counting function when we gave an exact version of the function in the previous exercise. The reason is that calculating the logarithmic integral is much faster than enumerating and counting the primes, as we did, and can be done for larger values of x. And by the way, use of the logarithmic integral as an approximation for the prime counting function was discovered by Carl Gauss as a fourteen-year-old boy!

You can run the program at http://programmingpraxis.codepad.org/9hg2s20m.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
5 Comments »

5 Responses to “Two Integrals”

  1. Programming Praxis – Two Integrals « Bonsai Code said
    January 11, 2011 at 9:48 AM

    […] today’s Programming Praxis exercise, our task is to write functions to calculate the exponential and […]

  2. Remco Niemeijer said
    January 11, 2011 at 9:48 AM

    My Haskell solution (see http://bonsaicode.wordpress.com/2011/01/11/programming-praxis-two-integrals/ for a version with comments):

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

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

liOffset :: Double -> Double
liOffset x = li x - li 2
  3. Graham said
    January 11, 2011 at 3:39 PM

    My Python solution. This blog and my free time studies are drawing me more and
    more towards Scheme and Haskell, but since there are two great solutions in
    those languages already I felt I should offer a solution in a different
    language. I’ve moved towards the newer “format” instead of the older printf
    style string formatting.

    #!/usr/bin/env python

from __future__ import division
import math

GAMMA = 0.5772156649015328606065


def exp_int(x):
    s = GAMMA + math.log(x)
    term, k, f = x, 1, 1
    while term > 1e-17:
        s += term
        k += 1
        f *= k
        term = pow(x, k) / (k * f)
    return s


def log_int(x):
    return exp_int(math.log(x))


def offset_log_int(x):
    return log_int(x) - 1.04516378011749278


if __name__ == "__main__":
    print "Li_offset(1e6) = {0:d}".format(int(round(offset_log_int(1e6))))
    print "Li_offset(1e21) = {0:d}".format(int(round(offset_log_int(1e21))))
# Output:
# Li_offset(1e6) = 78627
# Li_offset(1e21) = 21127269486616088576
  4. bubo said
    January 14, 2011 at 6:29 AM

    this is the best blog on the whole internet. thank you!

  5. GameSeven said
    January 19, 2011 at 6:36 AM

    Why 1e-17 is chosen as the stopping point in all of your code?

Leave a comment