Approximating Pi

July 29, 2011

In the two previous exercises we have seen three different methods for computing the prime-counting function π(n), all based on a formula of Legendre. We have also seen the difficulty of the calculation, as three eminent mathematicians got their sums wrong, and the difficulty continues even today with the use of computers, as Xavier Gourdon had to stop an attempt to compute π(1023) because of an error. Given the difficulty of calculation, in some cases it might make sense to calculate an approximate value of π, instead of an exact value. In today’s exercise we look at two methods for approximating the prime-counting function.

The first method dates to the German mathematician Carl Friedrich Gauss, who gives an estimate π(x) ∼ Li(x), the logarithmic integral—this is the celebrated prime number theorem, which Gauss figured out when he was fifteen years old! The logarithmic integral can be calculated by expanding the infinite series

\mathrm{Li}(x) = \int_0^x \frac{d\ t}{\log_e t} = \gamma + \log \log x + \sum_{k=1}^\infty \frac{(\log x)^k}{k!\ k}

to the desired degree of accuracy; somewhere around a hundred terms ought to be sufficient for arguments up to 1012. Beware the singularity at x=0.

In 1859, Bernhard Riemann, whose hypothesis about the prime numbers is one of the greatest open questions in mathematics, described a different, and much more accurate, approximation to the prime counting function:

\mathrm{R}(x) = \sum_{n=1}^\infty \frac{\mu(n)}{n} \mathrm{Li}(x^{1/n})

where μ(n) is the Möbius function that takes the value 1 when n is 1, 0 when the factorization of n contains a duplicated factor, and (−1)k, where k is the number of factors in the factorization of n, otherwise. There is no better way to compute the Möbius function than to compute the factors of n. As with the logarithmic integral, about a hundred terms of the infinite series ought to be sufficient to get a good approximation for arguments up to 1012.

Your task is to write functions to compute the logarithmic integral, the Möbius function, and Riemann’s R function, and to compare the results to the values you calculated for the prime-counting function in the previous exercise. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Advertisement

Pages: 1 2

In the previous exercise we studied Legendre’s prime-counting function. In today’s exercise we will look at two more prime-counting functions, one from Ernst Meissel in the late 1800s and the other from Derrick Lehmer in the mid 1990s.

Both formulas are based on Legendre’s original formula; in both cases, they simply rearrange the terms of Legendre’s formula to eliminate some work. We won’t try to give derivations, as they are complex and chock-full of Greek letters; if you are interested, Hans Riesel’s book Prime Numbers and Computer Methods for Factorization was the source for all three formulas.

Meissel: \pi(x) = \phi(x, c) + \frac{(b+c-2)(b-c+1)}{2} - \sum_{i=c+1}^b \pi\left(\left\lfloor\frac{x}{p_i}\right\rfloor\right), where b = \pi(\lfloor \sqrt{x} \rfloor) and c = \pi(\lfloor \sqrt[3]{x} \rfloor)

Lehmer: \pi(x) = \phi(x, a) + \frac{(b+a-2)(b-a+1)}{2} - \sum_{i=a+1}^b \pi \left( \left\lfloor \frac{x}{p_i} \right\rfloor \right) - \sum_{i=a+1}^c \sum_{j=i}^{b_i} \left( \pi \left( \left\lfloor \frac{x}{p_i p_j} \right\rfloor \right) - (j-1) \right), where a = \pi(\lfloor \sqrt[4]{x} \rfloor), b = \pi(\lfloor \sqrt{x} \rfloor), c = \pi(\lfloor \sqrt[3]{x} \rfloor), and b_i = \pi(\lfloor \sqrt{x / p_i} \rfloor) for a < i <= c.

Your task is to implement the prime-counting functions of Meissel and Lehmer, then compare timings with Legendre's prime-counting function of the previous exercise. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

The prime-counting function π(n) computes the number of primes not greater than n, and has been a matter of fascination to mathematicians for centuries. The Prime Number Theorem tells us that π(n) is approximately n / log(n); Carl Fredrich Gauss proposed the prime number theorem in 1792 when he was only fifteen years old, and Jacques Hadamard and Charles-Jean de la Valle Poussin both proved it, independently, in 1896. The first person to make a serious attempt at the calculation, except for trivial attempts that simply enumerated the primes and counted them, was the French mathematician Adrien-Marie Legendre at the end of the eighteenth century; his formula was correct, but he erroneously reported that π(106) = 78526, whereas the correct value is 78498. We will recreate Legendre’s calculation in today’s exercise.

Legendre’s formula works in two parts. First, he defined a function φ(x, a) that counts all the numbers from 1 to x inclusive that are not stricken by sieving with the first a primes. For instance φ(50, 3) = 14, because the numbers 1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 and 49 less than or equal to 50 are not divisible by the first three primes 2, 3, or 5. The φ function can be computed by the recurrence equation φ(x, a) = φ(x, a−1) − φ(x/pa, a−1), where φ(x, 1) is the number of odd integers less than or equal to x and pa is the ath prime number. The second part of Legendre’s formula uses the φ function to compute the value of π recursively: π(x) = φ(x, a) + a − 1, where a = π(⌊√x⌋).

Your task is to write φ and π functions and recreate Legendre’s calculation of π(1000000). When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Sum Of Two Integers

July 19, 2011

We have today another exercise in our on-going series of interview questions:

Write a program that takes a list of integers and a target number and determines if any two integers in the list sum to the target number. If so, return the two numbers. If not, return an indication that no such integers exist.

Your task is to write the indicated program. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

JSON: Writing Output

July 15, 2011

In the previous exercise we wrote a function to read JSON input and parse it into an object in the native language. In today’s exercise we write the inverse function.

Your task is to write a function that takes a JSON object and writes it in text format. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

JSON: Parsing Input

July 12, 2011

JavaScript Object Notation, JSON, is a system for passing structured data between computers, similar to XML but far simpler. JSON provides Unicode strings with a limited number of escapes, decimal numbers, the booleans true and false, the null value, arrays, and key/value dictionaries with string keys.

Your task is to write a parser that converts JSON data to an object in your favorite computer language. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Vedic Divisibility

July 8, 2011

Vedic arithmetic is a system of arithmetic devised in India for pencil-and-paper calculations. One of the calculations determines if one number is a divisor of another; the divisor must be odd and not divisible by 5; if the divisor is even or divisible by 5, factors of 2 and 5 should be removed before continuing. The calculation is a three-step process:

1) Determine the osculator of the divisor by multiplying it by the smallest number that will make the last digit 9, then strip the last digit 9 and add 1.

2) For each digit in the dividend, strip the last digit, then add that last digit times the osculator to the remaining digits in the dividend, giving a new dividend. Repeat as long as the dividend keeps decreasing.

3) Finally, divide the remaining dividend by the original divisor. If the remainder is 0, then the original divisor evenly divides the original dividend; otherwise not.

Here’s an example. For a divisor of 23, the number 23 is multiplied by 3, giving 69, then the last digit 9 is stripped, giving 6, and 1 is added, giving 7; thus, the osculator for 23 is 7. To test if 13174584 is divisible by 23, strip the 4 from the end of 13174584, multiply the stripped 4 times the osculator 7, giving 28, and add 28 to the stripped dividend 1317458 giving 1317486. Repeat, giving a new dividend of 131748 + 42 = 131790. Repeat again, giving a new dividend of 13179 (you can always drop a trailing 0, because 0 times the osculator is 0). Another iteration gives a new dividend of 1317 + 63 = 1380, from which we drop the trailing 0, giving 138. One final iteration gives 13 + 56 = 69. Now 69 ÷ 23 = 3, so 23 divides 13174584 and we are finished.

It is obviously easier for a computer to just perform the division and check the remainder; this method is intended for pencil-and-paper calculation (or even mental calculation).

Your task is to write a function that implements the Vedic divisibility test. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Big Numbers: Examples

July 5, 2011

For our final exercise in the big-numbers series, we use the library to implement some algorithms: display the factorials to 50, compute the factors of a number using trial division, determine if a number is prime using the Miller-Rabin algorithm, and compute the factors of a prime using Pollard’s rho algorithm.

Your task is to implement these algorithms using the big-number library. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Feet And Inches

July 1, 2011

Drawing programs measure distances in ten-thousandths of an inch, like 73.0185, but carpenters work in feet, inches, and thirty-seconds of an inch, like 6 feet 1 and 1/32 inches.

Your task is to write a program that takes a measurement in decimal notation and returns the measurement in carpenter’s notation; readers unfamiliar with imperial measurements will want to know that there are twelve inches in a foot. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2