The 10 hex digits of π starting at the 10 millionth place after the decimal point (is it really a decimal point if you are writing numbers base 16?) are 7AF58 63EFF. How do I know? Why, courtesy of my own implementation
of the Bailey-Borwein-Plouffe formula for computing hex digits of pi. This remarkable discovery allows you to compute hex digits of pi without computing all the digits to the left of it. My program is probably accurate out until 10⁷ decimal places and runs reasonably fast. When compiled with -ffast-math and -O3 on my AMD box, it computes the 10 millionth place in about 21.3 seconds. Not bad.

Addendum: Incidently, this program was written largely because I tried downloading this program from the LiteratePrograms wiki, and found that the “Basic Implementation” simply didn’t work for any n > 0.

Addendum²: The final digit is actually off in the above example. It should be 7AF58 63EFE. If you run the program to find the 10000001 digits, you get AF586 3EFEE (again, the final digit is off). We are operating at the limits of precision right around there. (Oddly, I seem to get the proper result on my Intel Xeon box at work. Not sure what’s going on there.) Whether I get the final digit right or not appears to depend on whether I compile with -ffast-math. With, it gets the final digit right. Without, it gets it wrong. Interesting.

Addendum³: You might want to read this paper for more information. The examples seem to differ in one position than the results I got, probably because my program addresses digit 0 as the first digit to the right of the decimal. No biggie.

Addendum⁴: If you replace “double” in the code with “long double”, apparently you can get 128 bit math, at least on my AMD64. It seems to work, I ran it out to 10^8th, and I’m getting the right digits:

bbp (using long double arithmetic)
CB840 E2192
421.034u 2.256s 8:20.93 84.4%   0+0k 0+0io 0pf+0w

I had need for a reminder of a bit twiddling hack that I had forgotten, and found this useful page that included many I don’t think I had seen before. These are mostly “too clever”, and should be used sparingly in your code. You’ve been warned.

Bit Twiddling Hacks

[tags]Programming,Hackery[/tags]

Need something to listen to on a really long trip? Librivox has just completed a huge audiobook version of the sixth edition of Charles Darwin’s classic On the origin of species. It’s a fascinating book, and well worth reading if you have any desire to understand the history of biological science.

LibriVox » The Origin of Species by Charles Darwin

[tags]Public Domain,Origin of Species,Charles Darwin[/tags]

While playing around with the GNU Bignum library and Mersenne primes, I decided to code up a simple program to compute pi to an obscene number of decimal places (I chose 10 million, a much larger number than any of my previous efforts). The GNU Bignum Library supports arbitrary precision floating point arithmetic as well, which made the Brent-Salamin (also called the Gauss-Legendre algorithm or the AGM algorithm) to compute pi. It’s very simple, and has quadratic convergence, which means that only 22 iterations are sufficient to generate 10m digits of accuracy. The resulting program worked the very first time I ran it. It takes 3m 11s or so on my AMD64 box, which makes it much, much faster than any other program I’ve written to date. (It still pales to the gmp-chudnovsky program (which I used to check my generated digits) but it’s not bad for 20 minutes of work and a lazy 73 lines of code. The GMPLIB is quite a nice piece of work.

Gauss-Legendre algorithm – Wikipedia, the free encyclopedia

And here’s my source code:

Barely commented, but consider the above link the comments.

[tags]Mathematics,Pi,Big Numbers[/tags]

Okay, most of the ones of these I’ve written have been at least a little obscure or obfuscated. This one is entirely straightforward, and is further unremarkable because you need to have the gnu bignum library to make it work. Shrug. It’s still kind of fun, and it demonstrates just how straightforward the library is to use.

Here’s the source code.

I found that I had to set my process limits up (apparently the print routines allocate memory on the stack, and that causes problems when your number has almost 10 million digits). In csh, that means “limit stacksize unlimited”. I had space to include that in the comment.

If you want something REALLY clever, try Fabrice Bellard’s obfuscated program to do the same thing. It’s a little too large for a .signature file, but it requires no additional libraries, and is positively diabolical. I must admit that I had a tiny bit of problems compiling this one, because the LL qualifier appearing on the second line seems to parse badly when it is broken up onto two lines. Still, very neat (and very fast too).

[tags]Prime Numbers,Mathematics[/tags]

This is the kind of paper I read when I have a lunch hour by myself at the local sushi bar.

[math/0110143] The Kruskal Count

The Kruskal Count is a card trick invented by Martin J. Kruskal in which a magician “guesses” a card selected by a subject according to a certain counting procedure. With high probability the magician can correctly “guess” the card. The success of the trick is based on a mathematical principle related to coupling methods for Markov chains. This paper analyzes in detail two simplified variants of the trick and estimates the probability of success. The model predictions are compared with simulation data for several variants of the actual trick.

[tags]Mathematics[/tags]

Through some convergence of factors (namely, my experience trying to read the accelerometers in the Wii remote, last night’s earthquake, and this mornings blear eyed checking of the Make Blog feed) I found the application below, which turns your Macbook or Macbook Pro into a seismograph by reading the accelerometer inside. I tried it on my budget Macbook. It works. Silly me, I though only the pros had accelerometers. Great, now I have something else to play with.

Suitable Systems / SeisMac

sieve completed... 0.000000 user, 0.000000 system
2**3 - 1 is prime 0.000000 user, 0.000000 system
2**5 - 1 is prime 0.000000 user, 0.000000 system
2**7 - 1 is prime 0.000000 user, 0.000000 system
2**13 - 1 is prime 0.000000 user, 0.000000 system 
2**17 - 1 is prime 0.000000 user, 0.000000 system 
2**19 - 1 is prime 0.000000 user, 0.000000 system
2**31 - 1 is prime 0.000000 user, 0.000000 system
2**61 - 1 is prime 0.000000 user, 0.000000 system
2**89 - 1 is prime 0.000000 user, 0.000000 system
2**107 - 1 is prime 0.000000 user, 0.000000 system 
2**127 - 1 is prime 0.000000 user, 0.000000 system 
2**521 - 1 is prime 0.012000 user, 0.000000 system
2**607 - 1 is prime 0.020001 user, 0.000000 system
2**1279 - 1 is prime 0.192012 user, 0.000000 system 
2**2203 - 1 is prime 1.092068 user, 0.000000 system 
2**2281 - 1 is prime 1.240077 user, 0.000000 system
2**3217 - 1 is prime 3.972248 user, 0.000000 system
2**4253 - 1 is prime 10.808675 user, 0.044002 system 
2**4423 - 1 is prime 12.324770 user, 0.044002 system 
2**9689 - 1 is prime 204.052752 user, 0.712044 system 
2**9941 - 1 is prime 225.430088 user, 0.744046 system 
2**11213 - 1 is prime 339.633225 user, 1.136071 system
2**19937 - 1 is prime 2574.836917 user, 8.732545 system
2611.599u 8.828s 1:27:19.73 50.0%       0+0k 0+0io 0pf+0w

Addendum: Sorry, I was too tired by the time I finished this last night to really comment. I wrote up a simple C program using the GNU multiprecision arithmetic library that used a simple sieve to find all numbers < 20000 which are prime, and then tested each one using the Lucas-Lehmer test to determine if it was, in fact, actually prime. It spews out the total runtime needed each time it finds one. You can see that it took 2611 seconds (43 minutes or so? my math might be off) of runtime (I was running another number cruncher at the same time, so the elapsed time was twice as long). The highest Mersenne prime I checked (2^19937 -1) was discovered in 1971. That means that in 43 minutes I was able to duplicate every computation on Mersenne primes done throughout human history up to 1971 myself. If I get some time later, I'm going to try to use the numbers generated to provide projections on how long it would take me to compute larger, more absurd primes using only my machine.

A pretty good jolt this time, really close (epicenter only a couple of miles away) with some pretty good rolling.

Recent Earthquakes – Info for event nc40194055

It seems that our neighbors to the north have sobered up on the whole anti-terrorism laws issue. Their Supreme Court recently revoked a law that allowed indefinite detention of foreign terrorism suspects violated their Charter of Rights and Freedoms.

BBC NEWS | Americas | Canada rejects anti-terror laws

Take the numbers from 82 stepping down to one, and write them all out together. You get a very big number.

828180797877767574737271706968676665646362616059585756555453525150494847464544\\
43424140393837363534333231302928272625242322212019181716151413121110987654321

It’s prime. What’s slightly odd is that there are no other numbers constructed the same way for n < 500 which are also prime.

[tags]Mathematics,Prime Numbers[/tags]

Addendum: I got the number from Prime Numbers by David Wells. I wrote a little Python script that did the Rabin-Miller primality test, and first verified that the number specified above was (at least very likely to be) prime, and then generated all similar numbers up to 500, and found they were all composite.

Addendum²: In 1964, Gillies discovered the 21st, 22nd and 23rd Mersenne primes using the Illiac II supercomputer. It took them 2 hours and 15 minutes to verify that 2**11213-1 is prime. I implemented the Lucas Lehmer test in Python, and verified the result in 40 seconds on my AMD64. I love Moore’s Law.

Perhaps this could be adapted to help me get exercise….

YouTube – Cat wheel II by Richard Norton

My daily sail upon the seas of the Internet blew me ashore on An Introduction to Scheme and its Implementation, a book on implementing Scheme by Paul Wilson. Scheme has a bit of nostalgia for me, having taken a class in Scheme from Will Clinger back in my graduate days at the University of Oregon.

Addendum: I previously blogged about Mark Feely’s presentation on writing a scheme to C compiler in 90 minutes. I’ll have to watch the video again.

A great collection of links to Classic Texts in Computer Science.

C’mon folks. It’s 2007 for pete’s sake. Why are there still programs which can’t be built on 64 bit machines? Yes, it can be a pain. When we ported the RenderMan with all its mumble mumble lines of code to the DEC Alpha, it was tiring, gritty, annoying, but absolutely necessary. And in the grand scheme of things, it’s not really that hard.