While exploring some related links from the previous video on analog computers, I ran across this very interesting link on creating fractals using video feedback.
YouTube - How to make fractals without a computer.
A couple of years ago, I blogged about H. E. Dudeney’s Amusements in Mathematics. Today, I noticed that Project Gutenberg had released a copy of The Canterbury Puzzles by Henry Ernest Dudeney - Project Gutenberg. This book has quite a few more nominally mathematical puzzles than its sibling. In particular, it introduces the game Kayles, which makes appearances in most of the books I have on combinatorial game theory such as Conway’s On Numbers and Games.
Addendum: Some information on strategy for Kayles can be found here.
A while ago, I got interested in reimplementing various versions of the Sieve of Eratosthenes. I eventually tacked together a threaded version that could calculate all the primes up to a trillion in about 20 minutes, and then kind of got bored. Today, however, there’s an interesting link to an article on Lambda the Ultimate about sieves. In particular, it makes the claim that the normal “sieve” code that people use to illustrate the effectiveness of functional programming languages isn’t actually a sieve at all. Interesting!
As some of you may have noticed, I occasionally like to write small programs to compute odd little mathematical curiousities. Something I hadn’t done in a long while was to use the Sieve of Eratosthenes to compute a bunch of prime numbers. I suspect that I wrote such a program very early in my exploration of computers, maybe thirty years ago. The basic algorithm is pretty straightforward, but takes O(N) space to find the primes less than N. It’s not hard to reduce the storage to N bits, and with a trivial bit of work you can reduce it to N/2 bits, and with a little more, you can reduce it to 4/15 N bits. That was fun to work out.
But last night I did something a bit different: I implemented a simple “segmented sieve”. Basically, the idea is that to generate all the primes up to N, you find all the primes up to sqrt(N), and save those. Then, you process the rest of the numbers in buffers of sqrt(N) at a time, by sieving each buffer with your saved primes. It’s a really simple idea, but makes sieving larger numbers more practical. I implemented the simplest possible version of this last night, and set it going to compute all the prime numbers up to 1E12. Here’s the log from the run:
[chessboard] % time ./bigsieve 1000000000000 > output 20580.222u 18.741s 5:43:39.10 99.9% 0+0k 0+0io 0pf+0w
And here’s the output file:
78498 primes less than 1000000 148933 primes less than 2000000 216816 primes less than 3000000 283146 primes less than 4000000 348513 primes less than 5000000 412849 primes less than 6000000 476648 primes less than 7000000 539777 primes less than 8000000 602489 primes less than 9000000 664579 primes less than 10000000 ...999989 lines deleted... 37607912018 primes less than 1000000000000
The program doesn’t include any of the basic optimizations, I suspect it woudn’t be difficult to make this thing a factor of 2-4 faster without adding a ton of code. I’ll probably see if I can do that over the next few days. It’s a useless but fun project.
I do like programming “for fun”, and this includes writing programs that, well, have been written hundreds of times before, but which I never have bothered doing. Often, this is just an exercise: freely available libraries often exist for many tasks which are considerably better than anything that I write. Nevertheless, I don’t consider it pointless. There are situations where having simple, self-contained code that is free from any particular licensing issues.
For instance, I’ve never really bothered coding up an FFT implementation. (Actually, this isn’t entirely true: I helped work on an i860 based assembly language version twenty years ago when I was at Princeton, but I must admit that I didn’t do much of the heavy lifting.) When forced to use one, I typically just fire up FFTW and use it. It’s a fast, relatively simple and flexible FFT library, and well worth using. But there are some circumstances where it seems like overkill. So, I sat down and implemented a basic radix-2 in-place FFT in C using the algorithm from Crandall and Pomerance’s Prime Numbers. Running it on my MacBook, for an 8192 length complex forward transform, I get a time of about 3.95 ms per transform. The same size transform running with FFTW is almost a factor of ten faster: about .413ms per transform. Of course, my code is a big stomping 43 comfortable, loafing lines long.
But this got me thinking about how it might be possible to speed it up. I have a couple of books on my shelf at work that will probably useful (I looted them from people who were more knowledgeable than me at DSP who, for some reason, decided to give them away, but I haven’t read them). Even doubling the speed of this “naive” implementation would make it significantly more useful. And, it seems like it could be a reasonably fun thing to tinker with.
It’s too late today to make any significant progress, but I did dig up a useful link.
Efficient FFT Algorithm and Programming Tricks
Addendum: Updating the so-called “twiddle factors” in the inner loop was being done in a stupid way in my code before. It really only requires a complex multiply in the inner loop, which speeds it up considerably. Now, my 8192 element FFT takes a mere 1.16ms on my MacBook, a speedup of nearly 4x, and now it’s within a factor of 3 of the speed of FFTW. I don’t like the way it really looks yet, but it’s better than most. I’ll tidy it up soon and archive it here for future consumption.
Addendum: Now that it’s basically working, I’m left wondering about numerical accuracy. Accuracy of the discrete Fourier transform and the fast Fourier transform would appear to be relevent. It appears that my choice of recurrence relationship for the twiddle factors probably results in rather significant errors (although I also suspect that it’s more than entirely adequate for the normal audio tasks which I use these things for typically).
Here is a good explanation of the FFT. I have used FFT libraries before, but I never really bothered to code one myself.
In fact, at the moment, it’s the largest known prime, with over 9.8 million digits. As part of my pi day celebration yesterday, I was trying to review how I might speed up my C code which calculates pi to large numbers of digits. Most of the fast ways rely on fast multiplication, utilizing the FFT algorithm. I wasn’t sure how that really worked.
So… I decided to write a program to learn and test my knowledge.
Rather than compute pi though, I thought I might try a somewhat different but similar task that relied on big number computation. I remembered that Fabrice Bellard had an obfuscated C code entry that used an integer Fast Fourier Transform to print out the largest (then) known prime. It’s really quite nifty. I decided to try to implement a similar thing, but rather than starting with his obfuscated code, I decided that I’d try to use the FFTW library to do the same.
It’s a tiny bit tricky: after all, we are using a floating point FFT to multiply large integers, and getting the precision issues nailed isn’t entirely obvious. My original idea was to represent the large numbers in base 100000 (giving five decimal digits per place) but for reasons which aren’t entirely clear to me, when computing the current largest prime, I seemed to run out of precision. Limiting myself to base 10000 seems to have cured the problem.
My program seems to work (at least minimally) now. I can compute Bellard’s number (2**6972593-1) in about 9.7 seconds, and the largest known prime (2**32582657-1) in about 53 seconds. Ironically, such numbers are really easy to output in base two: they are just lots of ones. But to convert them to base ten for formatting in the way us humans like is much more complicated.
The final output looks like:
12,457,502,601,536,945,540,085,550,157,479,950,312,279,598,515,115,184 284,367,047,566,259,111,523,599,739,738,055,975,960,661,684,593,910,041 988,688,211,130,870,620,428,490,430,485,634,271,939,241,796,764,631,759 ... ... 181631 lines deleted for brevity ... 826,726,604,958,937,322,582,512,072,612,621,443,114,535,641,869,584,273 577,446,330,457,465,821,333,212,445,737,104,635,692,000,092,659,011,752 880,154,053,967,871
I have double checked this against the official value on mersenne.org, and it matches precisely. I know that the program hides at least one gross inefficiency still: I suspect that I can actually compute the large prime in just about twelve seconds. But my head hurts for the moment, and I think I’ll pass for now.
Tomorrow is 3/14, known as pi day, or Albert Einstein’s birthday. How better to celebrate than with a script in pi-thon?
#!/usr/bin/env python
from itertools import islice
def g(q, r, t, i):
while True:
u, y =3*(3*i+1)*(3*i+2), (q*(27*i-12)+5*r) // (5 * t)
yield y
q, r, t, i = 10*q*i*(2*i-1), 10*u*(q*(5*i-2)+r-y*t), t*u, i+1
def pig():
return g(1,180,60,2)
p = pig()
pistr = ''.join([ str(x) for x in islice(p, 10001)])
print "%s." % pistr[0]
def byn(s, n):
return [s[x:x+n] for x in range(0, len(s), n)]
for line in byn(pistr[1:], 50):
for group in byn(line, 5):
print group,
print
I can never remember these formulas, so I wrote this program. I’m putting it here so I won’t lose it, and so others may benefit.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* $Id$
*
* Written by Mark VandeWettering.
*
* Any copyright would be pretty stupid. Use this code as you see
* fit, for any purpose, commercial or non-commercial, with or
* without attribution.
*
* I wrote this little test program mostly as documentation for some
* interesting recurrence relationships which I keep forgetting and
* have to dig up from time to time.
*
* To generate a continuous tone, you basically are repeatedly
* evaluating sin(n * ang), for n = 0... however many you need (or,
* alternatively, cos). Writing the loop in this naive way is pretty
* horrible though, it takes many cycles to evaluate a sine or cosine.
* Another way to do so is to compute the complex number m = cos(i *
* ang) (where i is the sqrt(-1)) and then maintain an accumulator which
* is repeatedly multiplied by m. Unfortunately, eventually numerical
* roundoff makes the resulting vector diverge.
*
* Luckily, there is a third option: to use a recurrence relationship
* for cos(n * ang) and sin(n * ang).
*
* From section 5.4 of Numerical Recipes: The Art of Scientific Computing:
*
* cos(n*ang) = 2 cos(ang) cos((n-1) * ang) - cos((n-2)*ang)
* sin(n*ang) = 2 cos(ang) sin((n-1) * ang) - sin((n-2)*ang)
*
* Somewhat amazingly, these recurrence relationships are _incredibly_
* stable. The program below implements both of them (generating a
* sin and a cosine in each iteration) using only two multiplies and
* two adds per iteration. On the 3.4Ghz Intel machine I have, this
* entire program runs in about 11 seconds, and the vectors traced
* out by the program are still unit vectors to the precision of
* a normal printf. You can try replacing the "doubles", with "floats",
* and you will see the precision isn't as good, but might be adequate
* for your needs. It should be noted, however, that the program runs
* in exactly the same time in either case (doubles and floats are both
* single cycle operations) with my machine/compiler.
*/
/* This test program runs for a billion cycles, which at the sample rate
* specified is just about a full day. The idea here is to see if the
* generated sines and * cosines continue to trace out unit vectors. */
#define NCYCLES 1000000000
#define SAMPLE_RATE 11025
#define FREQ (800.0)
main()
{
double ang = 2.0 * M_PI * FREQ / SAMPLE_RATE ;
double c0, c1, c2 ;
double s0, s1, s2 ;
double cm = 2.0 * cos(ang) ;
int i ;
c0 = cos(-2.0 * ang) ;
c1 = cos(-1.0 * ang) ;
s0 = sin(-2.0 * ang) ;
s1 = sin(-1.0 * ang) ;
for (i=0; i<NCYCLES; i++) {
c2 = cm * c1 - c0 ; /* c2 is cos(i * ang) ; */
s2 = cm * s1 - s0 ; /* s2 is sin(i * ang) ; */
c0 = c1 ;
c1 = c2 ;
s0 = s1 ;
s1 = s2 ;
}
for (i=0; i<100; i++) {
c2 = cm * c1 - c0 ; /* c2 is cos(i * ang) ; */
s2 = cm * s1 - s0 ; /* s2 is sin(i * ang) ; */
printf("%f %f %f\n", s2, c2, s2*s2+c2*c2) ;
c0 = c1 ;
c1 = c2 ;
s0 = s1 ;
s1 = s2 ;
}
}
/* $Log$ */
Addendum: Jim Blinn has a chapter on circle drawing in his book, JIm Blinn’s Corner which is probably relevent. It’s been some time since I read it (and I can’t find it on my shelf at the moment), but I believe it also points you to Minksy’s circle drawing algorithm that appeared in HAKMEM. It’s clever, but doesn’t really draw circles, and is no faster than the method above. What’s really odd is that it uses less storage for temporaries, which makes it look as if it were buggy.
Addendum2: Well, I coded that program up, just for fun. It turns out that you can use it as well. Gosper explains how to compute the desired value of epsilon in the Minsky’s algorithm, and unlike the code above, it is stable when evaluated using floating point arithmetic. Unfortunately, it generates a fairly notable ellipse, which doesn’t change over repeated evaluation, but isn’t very round either. As the speed the points go around the curve decrease, the curve becomes more circular.
Still, pretty interesting.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/*
* According to Gosper, Item 151 in HAKMEM,
* 1-epsilon^2/2 = cos(theta), where theta
* is the desired angular increment.
* This means that...
* 1 - cos(theta) = epsilon^2 / 2
* epsilon = sqrt(2.0 * (1-cos(theta))) ;
*/
#define SAMPLE_RATE 11025
#define FREQ 800.0
main()
{
float x = 1.0, y = 0.0 ;
float ang = 2.0 * M_PI * FREQ / SAMPLE_RATE ;
int i ;
float epsilon = sqrt(2.0 * (1-cos(ang))) ;
for (i=0; i<1000000000; i++) {
x -= epsilon * y ;
y += epsilon * x ;
}
for (i=0; i<1000; i++) {
x -= epsilon * y ;
y += epsilon * x ;
printf("%f %f %f\n", x, y, sqrt(x*x+y*y)) ;
}
}
Addendum3: I’m not sure what I really said about the “complex multiply” method really is true. Using the builtin gcc complex support, it runs in about the same time (even though I think it should have twice as many multiplies, perhaps it pipelines better?) and works just as well with respect to roundoff.
For completeness:
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
#define SAMPLE_RATE 11025
#define FREQ 800.0
main()
{
int i ;
complex m = cexp(2.0i * M_PI * FREQ / SAMPLE_RATE) ;
complex acc = 1 ;
for (i=0; i<1000000000; i++)
acc *= m ;
for (i=0; i<1000; i++) {
printf("%f %f %f\n", creal(acc), cimag(acc), cabs(acc)) ;
acc *= m ;
}
}
It’s the third game of the World Series, and the score is tied 0-0. If you can’t find any enjoyment in the reality of the game, perhaps the theoretical niceties of the game might provide some distraction:
Me? I enjoy both.
[tags]Baseball, Mathematics[/tags]
I pride myself on being able to carry out calculations in my head that cause most people to go for the calculator. Need to calculate a 15% tip? Divide the total bill by 10, and then divide by 2, and then add the results. Easy. Sometimes while exercising on a bike I try to factor four or five digit numbers in my head. Or try to remember how Conway’s Doomsday Algorithm allows you to compute the day of the week for arbitrary dates. It’s just a way to exercise the mind.
Another calculation that frequently (well, not that frequently, but still) comes up is to convert temperatures between Celsius and Fahrenheit. The way that I remember this is that zero C is 32 F, so to convert Celsius to Fahrenheit, you multiply C by 9/5 and then add 32. Of course, to do the reverse, you have to do the subtraction first (F-32) and then multiply by 5/9. They aren’t hard to remember, but I was shocked to find there is a nicer way that avoids the assymetry.
The Citizen Scientist has a cute article that shows the easier way. The fact you have to remember is that -40C is the same temperature as -40F. How does that help? Well, to convert from C to F, you
To convert from F to C, you:
So, a balmy 86F is 30C. 15C is 59F. Neat!
Having difficulty multiplying by 9/5 or 5/9? You could just use 2 and 1/2 if you’d rather, it introduces about a 10% error, which might be close enough. (You’d get 23C and 50F above). You can get closer if you subtract 10% when you multiply by 2, and add 10% when you multiply by 1/2, but frankly, I think that reintroduces the problems of assymetry that this revamp gives you. Probably easier (and more accurate) just to work on multiplying and dividing by 9 and 5.
Anyway, cool article.
[tags]Mental Math,Mathematics[/tags]
Today is 3/14, also known as Pi Day. Wish a “Happy Pi Day!” to your coworkers, and be forever branded as a math geek.
[tags]Mathematics,Transcendental Numbers,Pi[/tags]
For the Lego fans out there, check out K’s Legobrot: 3D Lego Fractals. I ran across it while looking for something completely different.  The program itself is rather unremarkable, but I still thought it was kind of cute.
[tags]Lego,Fractal,Mandelbrot,PostScript[/tags]
I’ve been reviewing some papers and math regarding the use of the Kelly Criterion in gambling, and found this link to Ed Thorp’s paper on the subject. It’s in PDF format, but looks very interesting.
[tags]Gambling,Mathematics,Information Theory[/tags]
