Today is Pi Day, but to hex with it
Did you know it’s possible to directly compute individual digits of pi? The BaileyBorweinPlouffe formula below can be used to compute the nth digit of pi… in hexadecimal:
The basic idea is that the nth term of the above summation corresponds to the nth digit of pi when written in hexadecimal — that’s what the (1/16)^{n} essentially means. Of course, since the rest of the term is hardly guaranteed to be an integer between 0 and 15 inclusive, each term “spills over” into the following decimal hexadecimal places, so you can’t just pluck out the nth term and expect it to give you what you want.
The trick, explained in more detail on the Wikipedia article on the formula, is to evaluate only as many terms as you need to include all the ones that could affect the nth digit, which basically means terms 0 through n+1. But since we don’t care what comes before that digit, we can multiply everything by 16^{n}, which shifts the decimal radix point to the position immediately before the digit we want, and use modular arithmetic to throw away the whole part, leaving just the fractional part we’re interested in.
Here’s some Haskell code that implements the basic idea:
{
 Simple implementation of using the BaileyBorweinPlouffe formula
 for computing individual hexadecimal digits of pi.

 Copyright (C) 2009 Paul Kuliniewicz <paul@kuliniewicz.org>

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2, or (at your option)
 any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 021111301, USA.
}
import Data.Ratio
import Numeric
import System.Environment
{
 Compute b^e mod m, using the method of exponentiation by squaring.
}
mod_exp :: Integral a => a > a > a > a
mod_exp modulus base = by_squaring 1 base
where by_squaring acc mult 0 = acc
by_squaring acc mult expo = let next_mult = mult * mult `mod` modulus
next_expo = expo `div` 2
in if odd expo
then by_squaring (acc * mult `mod` modulus) next_mult next_expo
else by_squaring acc next_mult next_expo
{
 Compute the "whole" part of one of the BBP subterms, where "b" is the
 particular constant used in the subterm being computed.

 n 16^(nk) mod (8k + b)
 E 
 k=0 8k + b
}
bbp_subterm_whole :: Integral a => a > a > Ratio a
bbp_subterm_whole n b = sum $ map kth_term [0 .. n]
where kth_term k = let modulus = 8 * k + b
in mod_exp modulus 16 (n  k) % modulus
{
 Compute the "fractional" part of one of the BBP subterms, where "b" is
 the particular constant used in the subterm being computed.

 oo 16^(nk)
 E 
 k=n+1 8k + b

 Actually, assume that only the first term is going to have any
 influence on the first fractional hex digit, so what's *really*
 computed is:

 n+1 16^(nk) 1
 E  = 
 k=n+1 8k + b 16 (8k + b)
}
bbp_subterm_frac :: Integral a => a > a > Ratio a
bbp_subterm_frac n b = 1 % (16 * (8 * k + b))
where k = n + 1
{
 Compute one of the four BBP formula terms. Since they only
 differ by one of the constants used, the same code can be used
 to compute each of them. (They also have different coefficients
 in the BBP formula itself, but we don't worry about that here.)
}
bbp_term :: Integral a => a > a > Ratio a
bbp_term n b = bbp_subterm_whole n b + bbp_subterm_frac n b
{
 Given a fraction, remove the whole part and return the
 fractional part.
}
frac_part :: Integral a => Ratio a > Ratio a
frac_part r = (numerator r `mod` denominator r) % denominator r
{
 Given a fraction, remove the fractional part and return
 the whole part.
}
whole_part :: Integral a => Ratio a > a
whole_part r = numerator r `div` denominator r
{
 Compute the nth hex digit of pi, where digit 0 is the first
 digit *after* the radix point. Uses the BBP digitextraction
 algorithm; see the Wikipedia article for a description of what
 this means and why it works.
}
nth_hex_digit_of_pi :: Integral a => a > a
nth_hex_digit_of_pi n =
let subsum = 4 * bbp_term n 1  2 * bbp_term n 4  bbp_term n 5  bbp_term n 6
shifted = 16 * frac_part subsum
in whole_part shifted
{
 Return a string containing hex digits "begin" through "end" of pi,
 where digit 0 is digit *before* the radix point. Note how the
 digit indexes are adjusted before calling nth_hex_digit_of_pi.
}
hex_substring_of_pi :: Integral a => a > a > String
hex_substring_of_pi begin end
 begin == 0 = "3." ++ hex_substring_of_pi 1 end
 begin > end = ""
 otherwise = concat $ map (hex . nth_hex_digit_of_pi) [begin  1 .. end  1]
{
 Helper function to convert an integer into a hex string.
}
hex :: Integral a => a > String
hex value = showHex value ""
{
 Get the range from the user and print it out.
}
main :: IO ()
main = do argv < getArgs
case map read argv of
[digit] > putStrLn $ hex_substring_of_pi digit digit
[begin, end] > putStrLn $ hex_substring_of_pi begin end
_ > mapM_ putStrLn ["usage: bbp digit",
" bbp first_digit last_digit"]
It will let you compute an individual digit:
$ ./bbp 15
d
Or a range of digits:
$ ./bbp 10 20
5a308d31319
$ ./bbp 0 50
3.243f6a8885a308d313198a2e03707344a4093822299f31d008
While the last example shows you can easily use the program to compute the first n digits of pi, calculating each digit individually is slower than other approximation techniques.
And since someone’s bound to be wondering:
$ ./bbp 700
6
So close. But given that nobody writes jokes in base 13, it shouldn’t be surprising that nobody writes songs about pi in base 16 either.
One Response
Thank you for that Paul. This is my first pi day as a number theorist, which is a lot more dissapointing than you’d expect. Thanks for being my ghost of pi day present.