Today is Pi Day, but to hex with it

Did you know it’s possible to directly compute individual digits of pi? The Bailey-Borwein-Plouffe formula below can be used to compute the nth digit of pi… in hexadecimal:

Bailey–Borwein–Plouffe formula

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 16n, 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 Bailey-Borwein-Plouffe 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 02111-1301, 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^(n-k) 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^(n-k)
-     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^(n-k)          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 digit-extraction
- 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.