Computers and calculating pi: A zillion-bit word?

http://www.straightdope.com/classics/a3_357.html

Cecil talked about the mathematics of computing pi, but not the computer science.

To calculate a more pedestrian value, say 1/2, you load a 1 into one word and a 2 into another word and then execute a divide instruction (more or less). This is oversimplified but the point is that the registers hold a word and the ALU does arithmetic on words.

So to calculate pi to a greater precision than that available in whatever word length is provided in hardware would seem to me to require a method like one of these:

  1. Calculate the next digit, and keep the result as a string of characters rather than a numerical value. But none of the math series I’ve seen would allow you to do that, AFAIK. This seems unlikely.

  2. Do some very low-level programming that combines multiple words to serve as a virtual very long word. This seems like quite a bit of gymnastics and would involve arithmetic on eventually hundreds of millions of words at once. Seems unlikely.

  3. Calculate the next chunk of significant digits, and save off the most significant digits from the running result as a set of characters. So at any one time you are just working with the group of least significant digits that are not yet fully determined. This seems at least possible.

But I have no idea what they really do.

Anyone know?

I’m not sure what the difference is between your options 1 and 3, but lots of math software has support for arbitrary precision numbers stored internally as strings of fixed-width numbers. (Not necessarily strings of ASCII characters, which may strike you as ridiculously inefficient and thus be why you dismiss the possibility, but I don’t see what’s so unlikely about storing arbitrary-precision real numbers as strings of fixed-width integers).

Seems I maybe misunderstood your objection; it was not that software wouldn’t let you manipulate arbitrarily large integers/arbitrary-precision reals as long lists of digits, but that the series, somehow, wouldn’t let you do it in this way. But I don’t understand the root of this objection. What makes this approach problematic for working with, say, the Leibniz series from the article?

In trying to use option 1, we find it is not possible to simply calculate one additional digit. Using the Liebniz algorithm, you are calculating like this:

1 - 1/3 + 1/5 - 1/7 + … converges on pi/4

so I presume that

4 - 4/3 + 4/5 - 4/7 + . . . . converges to pi

so we start with

4.0

the next term is -1.3… for an intermediate result of

3.66…

This presents a problem right off since you have an infinite number of digits already. You cannot represent this result exactly in decimal form so I’m not sure how to proceed. But it is clear that you can’t simply generate pi one digit at a time. You have to store an arbitrarily precise intermediate result after each operation. This also seems to blow my option #3 out of the water as well.

Storing them isn’t what seems unlikely, I just don’t know how you do arithmetic on them. But I am not knowledgable of the type of math software you mention; maybe that’s the key.

Well, note that all your partial sums are rational numbers; you can store them exactly by storing the numerator and denominator exactly. And the numerator and denominator can be stored exactly as just long strings of digits (presumably, for speed reasons, in a base 2^k system, where your word size is k bits, so that they’re stored as just long strings of words).

Now, as for how to get an approximation to a desired level of precision, the useful series generally come with known asymptotically vanishing error bounds: rules that tell you an upper bound for your distance from the actual limit after you’ve computed the Nth partial sum, with the upper bound approaching 0 as N approaches infinity. Keep going till the error bound is small enough, and you have enough precision. Alternatively (and essentially equivalently), get two series that both converge to pi (or whatever value it is you’re calculating); one approaching it from below (always) and one approaching it from above (always). As you calculate the two partial sums, they give you tighter and tighter bounds on the value, and you can stop when you achieve the amount of precision you want. [This is what Cecil was getting at with using inscribed and circumscribed polygons of many sides to calculate the circumference of a circle].

Now, I’m afraid I can’t recall off the top of my head a good error bound for the Leibniz series but I’m pretty sure you could derive one using Taylor’s theorem (at any rate, as Cecil remarked, it converges extremely slowly). I can also tell you that, while the Leibniz series neither stays below nor above pi, the series 1/1^2 + 1/2^2 + 1/3^2 + … converges to pi^2/6 from below, from which you can extract an approximation which converges to pi from below. Similarly, there’s also Wallis’s formula: the infinite product (2^2)/(2^2 - 1) * (4^2)/(4^2-1) * (6^2)/(6^2-1) * (8^2)/(8^2-1) … converges to pi/2 from below. So using either of those two and some formula which converges to pi from above, you could (slowly) calculate pi to any desired precision. Also, another way to get series which converge to pi from below or from above is to use lower or upper Riemann sums to calculate the integral of sqrt(1-x^2) from x = 0 to x = 1, which is equal to pi/4. (Thus, from this, if you know how to compute square roots to arbitrary precision, you get an easy method for computing pi to arbitrary precision).

But all of those are fairly slow convergers. Real pi hunters often use a variation on Leibniz’s formula called Machin’s formula, which converges rather faster. Or they use various other algorithm’s (e.g., the Gauss-Legendre algorithm, or Borwein’s algorithm). I’m no numerical analyst, so I can’t tell you much about the convergence rates, but the point is, they’re out there, useful asymptotically vanishing error bounds are known, so you can eventually tell with any of these methods that you’ve reached the desired level of precision, even if it’s not something as simple as “Each step gives you precisely one more digit”. And like I said, it’s easy enough to store numbers exactly enough to do the calculations: if the calculations only involve rationals, you can store all those perfectly; similarly for if, say, the calculations only involve sums of K-th roots of rationals, you could store those all perfectly and do arithmetic with them, with a properly designed encoding system. If your calculations are more hairy than would allow for any of that, you can do interval arithmetic (store perfect rational lower and upper bounds on all your numbers, treating them as fuzzy intervals rather than precise values) or, essentially equivalently, you can do some computation to figure out what level of precision you need to work with in your calculations to get what level of precision out. Like I said, I unfortunately don’t know all the formulas telling you exactly how much you have to do to get out how much precision, but people who do this sort of thing do know such formulas.
[Since I’m no expert in this area, anyone who is a numerical analyst, or has any other reason to want to correct me, feel free to go ahead and do so on any point. But I think I’ve represented the situation quite accurately.]

One might also add that, historically, not all computers have been binary, or register-oriented. The most popular computer of the early 1960’s, for example, the IBM 1401, performed decimal arithmetic in a storage-to-storage mode that was limited only by the amount of memory.

There are in existance thousands of infinite series (and products) that converge on pi. Some of these converge very slowly while others quickly. Some are difficult calculations, while others are straightforward. Which of course makes some more useful for calculating zillions of decimal places than others.

I forget the exact details – they are in a book I have at home – “The book of curious and interesting numbers” bu David Wells. He mentions a series where the denomenators are powers of 2 (in repeating blocks of four terms iirc). This is a fantastic series to work with since it enables one to calculate an arbitrary hexadecimal digit of pi. You do not need to calculate the billionth digit before calculating the billion and first. This is quite an advantage. And it completely does away with a whole lot of problems of storing and representing the expansion once calculated.

Of course converting from hexadecimal to base ten does resurrect some of these difficulties, but there are well established algorithms for doing this and it is rightly treated as a fairly routine task.

You are almost certainly referring to the Bailey-Borwein-Plouffe formula. It is indeed a very nifty way to extract isolated binary/hexadecimal digits of π, without having to compute/store any of the previous digits. But I think you’ve been a bit glib in passing over the difficulties in converting to decimal; so far as I know, there’s no easy way to use this to get the value of an isolated decimal digit without first computing the whole string of corresponding hexadecimal digits up to the corresponding one.

Furthermore, so far as I can see, the BBP algorithm still does require some use of arbitrary-precision arithmetic. And, if your goal isn’t to grab isolated digits but to compute many, many digits (e.g., π up to a desired level of precision), I don’t think there’s much advantage to the BBP algorithm over other, older methods. But it still is a very nice discovery.

You are right. I have been a bit glib in the base 16 to base 10 conversion. I think it is routine, but certainly not trivial. If that makes any sense.

As has often been pointed out, calculating pi to a particular precision is actually not that useful from a practical standpoint. But people pursue the task for other reasons. One is to investigate the normality of the digits, which may as well be done in hexadecimal as base ten for any statistical analysis. The algorithm (and we are probably talking about the same one) allows one to compare a strings of digits from different places in the expansion without having to compute all digits in between.

Of course, if one is attempting to investigate pi’s normality or otherwise in any base (or all bases) then base conversions will be required. And not just to base 10. This definitely is not a routine task and it will require a full expansion of pi.

Heh; I understand and agree completely.

Incidentally, for some reason, in this text editing box (created by clicking Reply), there are two spaces between your sentences, even though your post displayed with only one space between them. Is the board software actually bothering to override your preferences for how to space sentences, or what’s going on? Total minutiae, but my curiosity burns strong…

It’s not the board, it’s your browser (and everyone else’s). The default behaviour for HTML is to compress whitespace. You have to make a special effort if you want whitespace to remain uncompressed. But that doesn’t stop those of us long in the habit from leaving two spaces after a period, anyway.

Ah, I see. I somehow assumed the board software would go to great pains to shield our meticulously thought out and formatted text from being put through the HTML interpretation wringer (after all, I can’t just type <a href =“www.google.com”>an HTML-style link</a> and have it work out), but apparently not quite as such.

The book Numerical recipes in C contains math routines that do exactly that, achieving arbitrarly high precision. Calculation of pi to as many digits as you wish is provided as an example of thier utility.

A thought just occurred to me, and now I’m deeply ashamed to have been so stupid as to not realize just how very easy it is to extract asymptotically vanishing error bounds for the Leibniz series. The terms alternate in sign while decreasing in magnitude, which is precisely to say that the partial sums alternate between an increasing series and a decreasing series, and thus between approaching the limit from above and below. So to calculate π/4 to a desired level of precision, just calculate successive partial sums of the Leibniz series (using exact rational arithmetic, as explained above) until they differ by less than that level of precision. Simplicity itself.

Of course, Leibniz’s series converges so very slowly, but the same idea extends to all the other faster-converging series with terms of alternating sign and decreasing magnitude. Man, I am quite embarrassed to not have realized and mentioned this earlier.

In fact, someone once derived, on the board, a faster algorithm based on Leibniz’s series. Apparently, if you construct another sequence from the averages of successive partial sums, that sequence alternates above and below pi, too. And if you take the average sequence of that, it alternates as well, and so on. Take it to enough levels, and you can converge on pi much more quickly than the original series ever did (though probably still far short of the best known methods).

As a minor note, I wonder how many digits they got in the original-series Star Trek episode “Wolf in the Fold.” A state-of-the-art mid-22nd-century computer going all out for, I dunno, a hour or so? How long did it take to drive Redjac out of the system?

[sub]glavin![/sub]

Is that one of those awful plots where they have to run the computer all-out so there is no resources left for <something bad> and they concoct some ridiculous scheme to keep the computer busy with something poetic like calculating Pi? You want to keep a computer busy? Here:



int main(void)
{
   int i = 1; 
   while(i--)
   {
        fork(); 
        i++; 
   }
}

Bingo. Fortunately for the good guys, the bad guys’ computers run under single-task operating systems.

That’s good, but you hardly need that i variable, do you?

It does nothing! (But the bad guys won’t know that!) It’s just in case your futuristic compiler is armed with a supercompiling optimizer. :wink: You know, so it doesn’t optimize the fork() out of the loop.

You’d see the opposite happen: the variable i could be optimized away (it does nothing!), but not the function call. In fact, it would be illegal for the compiler to drop it in this case. The function “fork” might have side effects, for all it knows. You couldn’t drop it without changing the meaning of the program.

And of course, as we humans know, “fork” does have some mighty big side effects.