Happy Pi Day!

Pi Day is celebrated on March 14, i.e. 3-14. Indiana, however, once came close to celebrating it on 3-20 exactly.

Actually, that’s not quite true. The text of the 1897 bill doesn’t come out and say it directly, and it gives several different derivations of pi, none of which are consistent with each other. Values include pi = 3.2, pi = 4, and pi = 16√2 / 7 ≈ 3.23. Who would have thought the mathematical crank who came up with it didn’t check his work?

The intro of the bill is also interesting, and tends to get overlooked in discussions about it:

A Bill for an act introducing a new mathematical truth and offered as a contribution to education to be used only by the State of Indiana free of cost by paying any royalties whatever on the same, provided it is accepted and adopted by the official action of the Legislature of 1897.

The crank behind the bill apparently planned on trying to collect royalties on the “correct” value(s) of pi, after giving Indiana a royalty-free license. Of course, this is nonsense, since you can’t copyright, trademark, or patent a fact, so there’s no way you can compel anyone to pay you royalties for it. Passing a bill to grant a state royalty-free access to a fact that isn’t even true is completely stupid in two independent ways.

Needless to say, it passed the Indiana House unanimously before getting shelved in the Indiana Senate, demonstrating the cluelessness of Indiana’s elected officials in both mathematics and intellectual property law. If not for the intervention of a Purdue mathematics professor, Indiana students today might be making funny-shaped circles to comply with state law.

The whole ordeal made the Indiana government a laughingstock, and since then all state governments have been careful to consult with experts before taking action to reject basic facts. Oh, wait.

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
- 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

Or a range of digits:

$ ./bbp 10 20
$ ./bbp 0 50

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

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.

Does your computer appreciate ? Day?

It’s that time of year again, when we all reflect on what the ratio of a circle’s circumference to its diameter means to us. But have you ever stopped to ponder what it means… to your computer?

For C programmers at least, the answer might first seem to be found in the the standard library’s math.h header:

# define M_PI           3.14159265358979323846  /* pi */

The M_ prefix denotes it as a constant provided by the math part of the standard library, though you can pronounce it as “mmmm, pie” if you like. Either way, though, that certainly seems like a perfectly reasonable approximation to the transcendental number.

But is it?

To find out, first we need to consider how your computer represents real numbers internally. (OK, technically only rational numbers, and even then only a subset of them.) Most likely, your computer uses the standard floating point representation, which works pretty much like a binary version of scientific notation.

As you’ll recall, since scientists often work with really large or really small numbers that are obnoxious to write in normal notation, they like to move the decimal point around to make the number between 1 and 10, and then multiply the result by the appropriate power of 10 to compensate. For example, instead of writing the number of atoms there are in a gram of 12C as about 602,214,790,000,000,000,000,000, they’d write 6.0221479 × 1023.

Floating point works pretty much the same way, but in binary instead of decimal. For example, the number 13.40625 would be 1101.011012 (using the subscript 2 to denote a binary number), or 1.101011012 × 23.

Of course, floating point stores both the fractional part and the exponent as binary integers, but the precise details aren’t important for this discussion. What is important is how many bits are used to store the fraction and the exponent, and in C we have a few options:

Type Sign bit Exponent bits Fraction bits Total bits
float 1 8 23 + 1 32
double 1 11 52 + 1 64
long double 1 15 64 80

The sign bit denotes whether the number is positive or negative. The more exponent bits we have, the larger and smaller (in the “distance from zero” sense) the numbers we can represent, since we can move the binary point (not a decimal point!) farther to the left or right. The more fraction bits we have, the more precisely we can represent a number, just like 2.718281828459045 is more precise than 2.718.

If you’re paying attention, you might protest that the rows in that table don’t add up. How can we have 1 sign bit, 8 exponent bits, and 24 fraction bits in a 32-bit value, where clearly that’s 33 bits? Easy. In scientific notation, remember how we move the decimal point until the fraction part is between 1 (inclusive) and 10 (exclusive)? In floating point, we move the binary point until the fraction part is between 12 (inclusive) and 102 (exclusive). That means that the fraction part always looks like 1.something2. Since that first bit is always a 1, we don’t need to bother storing it. Woo, free bit! That’s why I wrote 23 + 1 in the table and not 24: 23 stored bits plus 1 implied bit.

The long double type is the oddball, since it’s much more platform-dependent than the other two. On x86 machines, it’s the 80-bit floating point value that the FPU uses internally. Other architectures will probably use something different. Also, the long double format explicitly stores that guaranteed-to-be-a-1 bit at the front of the fraction, presumably because that makes the FPU’s work easier.

Now, with that background out of the way, back to our question: just how accurately can our computer represent ?, without rolling our own data type? Here’s how the three types fare, compared with the true value:

Type Value
float 1.100100100001111110110112 × 21
double 1.10010010000111111011010101000100010000101101000110002 × 21
long double 1.1001001000011111101101010100010001000010110100011000010001101012 × 21
? 1.10010010000111111011010101000100010000101101000110000100011010011…2 × 21

(If you want to check for yourself, here’s the source code I used.)

Now, I may be a master of computer science, but even I have trouble looking at those numbers and seeing how they compare to our buddy M_PI, a.k.a. 3.14159265358979323846. Fortunately, it’s not too hard to compute the maximum rounding error in our floating point values. In general, the magnitude of the rounding error in base b is the number with digit b/2 in the position immediately after the last digit in the rounded number, and zeroes everywhere else. If it were otherwise, the number would’ve been rounded to something else. For a decimal example to make that clearer, the error in rounding ? to 3.14 is at most ±0.005: the true value of ? must be between 3.135 and 3.145, otherwise it wouldn’t've been rounded to 3.14.

Crunching the numbers for our three data types, we get:

Type Maximum error in ?
float ±2-23 ±10-6.9237 ±0.0000001192092895507811309
double ±2-52 ±10-15.6536 ±0.0000000000000002220446049
long double ±2-63 ±10-18.9649 ±0.0000000000000000001084202

To make it clearer still, let’s look at what the decimal values of ? in those types look like, both written out to 21 digits and truncated before we run into any errors:

Type Value Truncated Accurate Digits
float 3.14159274101257324219 3.141592 7
double 3.14159265358979311600 3.141592653589793 16
long double 3.14159265358979323851 3.141592653589793238 19

Moving from float to double roughly doubles the precision we get before losing accuracy, and moving on to long double adds a few more.

For comparison:

Constant Value Accurate Digits
M_PI 3.14159265358979323846 21
M_PIl 3.1415926535897932384626433832795029L 35

It turns out that M_PI is more precise than any of our options for storing it in a floating point value! That makes the GNU extension M_PIl, specifying a long double version of the constant, a bit silly, since we can’t use those extra 14 digits anyway! The only functional difference between the two is the L suffix at the end, telling the compiler to treat the constant as a long double instead of an ordinary double. Without that suffix, assigning M_PI into a long double would still only give us the precision of a double, with those extra fraction bits set to 0.

Of course, in C, storage limits of the different data types is machine-dependent. It’s certainly possible for those same definitions to be used on a machine with larger floating point types that can make use of those extra digits. But for poor old kryten here, it is not meant to be.

For extra ? Day fun, consider that many programs don’t want to depend on having a math.h header recent enough to define M_PI for them, and so define their own version of the constant as a backup. A Google Code search for C programs with “#define M_PI” in them turns up no shortage of hits. Amusingly, many programs have their own ideas as to how many digits M_PI should be defined with.

One particularly fun one is this hit for a 97-digit version of M_PI, if only for the comment in front of it:

 * Ok, so some systems require magic incantations for M_PI to be defined.
 * It's not important enough to worry about.
#ifndef M_PI
#define M_PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117

Also amusing are the “you fail at math” examples, like this one:

#ifndef M_PI
#define M_PI 3.14159263

OK, it’s not as bad as the Indiana state legislature’s infamous attempt to implicitly define ? as 3.2. To say nothing of defining PI to be 12. But still.

So now you know: the great thing about ? is that there’s so many different values for it lurking in your computer. Of course, all the x86 assembly programmers are laughing about this whole mess, since they’ll just use the FLDPI instruction to have the FPU use its value for ? and be done with it.

Comments Off