Pi from Nothin'

Riemann Hypothesis - Numberphile
Stu turned me on to an obscure little math trick for determining if a number is a prime number or not. It's not generally useful as it involves factorials, which get really big really fast. But there is a computer math package that can deal with numbers of any size, as long as you have enough computer memory to accommodate them, so I thought I would take Stu's formula and the GMP math package and see if I could make them play together.
     Dealing with really long numbers is not very difficult, you just have to keep track of overflows. Any half competent programmer should be able to write a set of routines to perform the basic operations, but it's a common enough job that some people got together and wrote a complete library of functions and procedures.
    Like any bit of software there are rules for how to use it. For me the easiest way to learn these rules is find a sample piece of code and start messing with it, so I went looking and I found this: randPi.c by Mitch Richling, which to purports to make Pi out of nothing but random numbers. He uses the GMP library to deal with his really long random numbers, so it was good example to start with. Plus it was weird. Pi from nothin'? How can that be?
     It depends on the observation that between any two random numbers, there is a statistical probability that they are relatively prime, i.e. the only factor they have in common is one. 3 and 4 are relatively prime. 3 and 6 are not, they have a factor of 3 in common. 4 and 6 are not, they have the factor 2 in common. Two numbers by themselves won't tell you much, but if you have big pile of them and start counting the number of pairs that are relatively prime and those that aren't, you end up with a ratio that (if you munge it just a bit) looks just like Pi.
    I got to thinking about this and realized you don't really need random numbers. After all, random numbers are just one number out of a range of numbers. You run long enough and you will cover all numbers in the range. You should be able to get the same effect by just testing the relative primeness of all the numbers in a range. I wrote a small program to test this and it indeed works. Running a cross tabulation of all the numbers from 1 to ten thousand resulted in value of 3.141534, which is not bad considering I made it from nothin'.

    Okay, we have GMP library and we think we know how to use it. Let's test Stu's formula, so I wrote another little program, it works and so Stu's formula seems to be valid. My program only goes up to 31, but that's only because we are starting to fill up the screen with the giant numbers that come from computing factorials. 31! (31 factorial) is like a zillion digits.

P.S. I have decided I don't like github because they don't handle tabs. Google Drive / Docs has a viewer that works okay and the two programs are stored there. It doesn't support highlighting, and there may be some argument over tabs, but it's not bad, and it doesn't require any extra fooling around.

Transcontinental Airmail

United States Contract Air Mail routes 1925–1930. The numbers on the routes are the Contract Air Mail (CAM) route numbers. The original, government-run transcontinental route is shown in a thicker, lighter-coloured, dotted line running East–West between San Francisco and New York.

Comrade Misfit furnished a link and down the rathole I go, wherein I eventually came across this story by John Schamel about how coast-to-coast airmail got started. It's just a little nuts. I found it on a couple of different sites, but because I didn't like the formatting I made my own copy which is more readable.

In the midst of this story I found this paragraph where he talks about beacons.

What resulted was the first ground based civilian navigation system in the world. Beacons were positioned every ten miles along the airway. At the top of a 51-foot steel tower was a 1 million candlepower-rotating beacon. Pilots could see the clear flash of light from a distance of 40 miles. Also at the top of the tower were two color-coded 100,000 candlepower course lights. These pointed up and down the airway. They were colored green, signifying an adjacent airfield, and red, signifying no airfield. The course lights also flashed a Morse code letter. The letter corresponded to the number of the beacon within a 100-mile segment of the airway. To determine their position, a pilot simply had to remember this phrase – “When Undertaking Very Hard Routes, Keep Direction By Good Methods” – and know which 100-mile segment they were on.

 Now why would they choose that particular sequence of letters? Let's see what the Morse Code actually looks like.


Looks like they picked codes that had some kind of rough approximation to the numbers. They also have the advantage of being shorter than numbers. Numbers in Morse code are all five blips, none of these has more than four.

I suspect that the erection of these towers caused several locations to be christened 'Beacon Hill'.

Computing Pi

A while back Stu put up a post that had a newish formula for computing the value of pi. Just for grins, I thought I would see if I could implement the formula in code, that is, write a computer program to compute the value of pi using this formula. I mean, you never know about these things. Sometimes you get formulas that look nice, but when you try and turn them into code you find out there is one piece of psycho-babble in the formula that cannot be implemented in any kind of straight forward manner. Anyway, I tried it and I was successful, and it did indeed compute an accurate value of pi. Bear in mind that accurate is a relative term.

It is claimed that this particular formula is unique in that it can compute any digit of pi without having to compute any of the previous digits. I tried to figure out how this was being done. I failed. The formula computes a series of values and each value is added to the previous total. Each successive value makes the total a little bit more precise. Each successive value is several times smaller than the previous one, so if the first value computes the first digit, then the second value will have no effect on the first digit, but it will have an impact on the second digit, and so on. However, each value that is computed is an insane decimal number in it's own right, and so the string of digits to the right of the decimal point goes on indefinitely. So while computing the Nth value of this series will not effect earlier digits, the Nth digit is still going to be the sum of all the Nth digits from all the previously computed values.

I eventually figured out what they meant. Typically these formulas include a section that needs to be repeated several times in order to achieve the required accuracy. Most of these formulas require that after you have completed the necessary number of steps and have arrived at some total, you now subject this total to another formula for further manipulation. This newish formula does not require this second step.

I took a look at the article on Pi in Wikipedia and it has a whole list of formulas for calculating Pi. And that's when it struck me that although we have these formulas, we have no explanation for how they work. Matter of fact, all we have is someone's word that they do work. And how do you know what the 10,000th digit of Pi is anyway? You can only verify the value empirically (by measuring a physical circle) out to ten or maybe 20 digits. Beyond that we are talking the intellectual equivalent of peacock feathers: very pretty but totally useless. I imagine each one of these formulas was the result of someone's doctoral thesis in mathematics and the proof of each one is probably 400 pages of inscrutable squiggles. Presumably they can all be tested against each other by writing and running a computer program.

One of the formulas was developed by an Indian prodigy about 100 years ago. It is able to compute the value of pi to the limits of my machine in three steps. The new formula takes about a dozen steps.

So now I'm wondering if I can compute the value of pi straight up from nothing. The standard way to do something like this is to divide a circle into a bunch of triangles, like cutting up a pie, and then add the lengths of the bases of all the triangles. The smaller you make the triangles, the more closely the bases of these triangle will approximate the circle, and the closer you will get to the true value of pi. So I figured out how to calculate the base of ever smaller triangles using the pythagorean formula, and put it in a program, and in a couple dozen steps or so I had a pretty good approximation.

Links:

• C Language source code for the demo program
• Output from demo program
• Explanation of geometry/algebra used for tangent computations
• Explanation of geometry/algebra used for chord computations
• An irrelevent Stu story about Srinivasa Ramanujan

I could find Stu's story about BBP.

Dyanmic Two Dimensional Arrays in C

In computer programming, a list of numbers can be stored in an array. The size of the array can be decided in advance when you write the program, or you can create the array at run time to be any size you want. A simple list is called a one dimensional array.

A two dimensional array is more like a spreadsheet. You have rows and columns. The C programming language supports two dimensional arrays that are created when you write the program, but it does not support dynamic creation of two dimensional arrays.

However, a two dimensional array can be viewed as a list of lists, and because C treats pointers and arrays pretty much interchangeably, it is possible to dynamically create a data structure that functions as a two dimensional array. The computation of the actual address in memory of the data field being accessed may not be as elegant as that done for static two dimensional arrays, but it may be faster as it only involves addition and no multiplication. There is some overhead in storage for an extra row of pointers, but that should not be a problem for today's computers.

I did this once before, but the code wandered off. I am playing with matrices these days, so I decided to redo it. It is not particularly complicated, but you do need to get the details right to make it work. I put the essentials in a demo program you can examine should you be so inclined.

More π

So I was thinking about that program for calculating the value of π that I posted yesterday, and I realized that the algorithm was quite capable of calculating the value of pi to an indefinite number of digits. The code would take some work, quit a bit of work actually, but feasible. The algorithm calculates a new value of x with every iteration. This new value of x starts 4 bits to the right of the place of the last value of x. So you could keep adding these new values of x to this string of bits for however long you wanted to. Eventually you would get tired or bored, or run out of disk space, but you could go quite a while. The problem you would run into is that each value of x would be an infinitely repeating pattern, so you would have to keep replicating this pattern as you went, and with each iteration, the pattern would get longer. It might be necessary, or even easier to keep a separate pattern for each new value of x. Then when you get a new value of x, you add the appropriate section of each pattern to it and append that value to your infinitely extending value of pi. There is going to be some point where the addition of the new x is not going to propagate any farther towards the binary point. Then you can say you know the value of pi to that point and can safely store it away. Recognizing how many bits go into each pattern is going to be another trick. I imagine someone has already figured out how to do that.

π, again

Yesterday was pi day and Stu put up a post with a formula for calculating the value of pi, so naturally I had to try it out. This one promised to be a little better than the traditional one. I had written a program to implement that algorithm a while back, and while it worked, it was very slow to close in on the value of pi, requiring hundreds of iterations. This formula closes in on pi in only 8 iterations. Well, it gets as close as it can. It is limited by the format of double precision floating point numbers.

/* Simple program to calculate the value of pi.
  Original formula:

           /   4        2        1        1    \
  16^-K * (  ------ - ------ - ------ - ------  )
           \ 8K + 1   8K + 4   8K + 5   8K + 6 /

  The caret (^) indicates raising to a power: 16^K means 16 to the power of K.
  Multiplying by 16^-K is the same as dividing by 16^K

  A double precision floating point number has 64 bits (8 bytes).
   1 bit for the sign
  11 bits for the exponent
  52 bits for the mantissa

  Three decimal digits require approximately 10 bits. 2^9 = 1 << 10 ~= 1024
  We were able to get 13 decimal digits, which works out to about 43 bits.
  The double angle brackets (<<) indicate a bitwise shift,
  in this case to the left.
*/

#include <stdio.h>
#include <math.h>

int main(void)
{
   int       k;
   int       n;
   double    pi;
   double    x;

   pi = 0;
   for (k=0; k<8; k++)
   {
       n = 8 * k;
       x = ( 4.0/(n + 1)
           - 2.0/(n + 4)
           - 1.0/(n + 5)
           - 1.0/(n + 6))
           / (1 << (4 * k));    // 16^K
       pi += x;
       printf("%16.12f    ", x);
       printf("%16.12f\r\n", pi);
   }
   return 0;
}


Results:
      x                   pi
3.133333333333      3.133333333333
0.008089133089      3.141422466422
0.000164923924      3.141587390347
0.000005067221      3.141592457567
0.000000187893      3.141592645460
0.000000007768      3.141592653228
0.000000000345      3.141592653573
0.000000000016      3.141592653589

No sound Epox motherboard Ubuntu 8.10

I've been fooling around with my Linux box trying to get the sound working. As a practical matter it's basically a waste of time. I could go down to the local computer shop, buy a new sound board for $10, plug it in and it would probably work. But learning the in's and out's of Linux might be a useful skill, so I have taken to investigating the problem, and this means being assimilated into the Linux culture. If I keep it down to an hour or two a day, it is not too bad. I learn a little bit, I make a little progress, I get disgusted with the convoluted way this, that or the other was implemented. Then I stop for the day, and let the disgust dissipate, and the next day I am ready to try a little more. What I wouldn't give for a design document that explains how this all supposed to work. There probably is one somewhere in the thousands (3,715 to be exact) of files in the ASLA subsystem, but I'll be durned if I know where it is. Meanwhile I am sorting it out by stepping through the code for aplay with gdb (Gnu Debugger).

Most of the code looks like it is involved in setting things up, opening the file, figuring out what format it is, what device is going to get this data. Then at some point it flips a switch, something magic happens, and sound comes out of the speakers. Or, as in my case, not.