Allan’s Musings

The Ulam Spiral

Posted by Allan Peda on May 11, 2009

A few years ago I wrote a perl program to generate a map of prime numbers based on on a what is known as an Ulam Spiral. Sure enough the patterns are readily recognizible. The origin is barely discernible in the center of the image below as a lone red pixel, though you’ll need to zoom in to see it.

``` #!/usr/bin/env perl # -*- cperl -*- BEGIN { \$^W = 1 } use strict; use diagnostics; use warnings; use Carp; use English; use Math::Trig; use GD; use constant PI => 4*atan(1); # Allan Peda # June 30, 2007 # Ulam prime number spiral # http://mathworld.wolfram.com/PrimeSpiral.html # http://en.wikipedia.org/wiki/Ulam_spiral # # we iterate over a grid, but we want to follow a spiral. # Direction is cartesian based x being horizontal motion, # and y being vertical. Up and right are postive motion. unless(defined(\$ARGV)){ die "Please provide a png file to save to."; } my \$pngfile = 'ulam_spiral.png'; if(\$ARGV =~ m/\.png\$/i){ \$pngfile = \$ARGV; } else { print "Supplied file not suffixed with .png defaulting to \"\$pngfile\".\n"; } # how many turns to make in the spiral my \$num_corners = 1600; my @turning_points = mkturns(\$num_corners); # which directions do we follow my \$max = \$turning_points[\$#turning_points]; print "Generating map for first \$max primes.\n"; my (\$x_width, \$y_width) = (0,0); \$x_width = \$y_width = int(sqrt(\$max+0.5)); print "Dimensions are: \$x_width x \$y_width\n"; my %primes = %{genprimes(\$max)}; my \$im = new GD::Image(\$x_width, \$y_width); my \$white = \$im->colorAllocate(255,255,255); my \$black = \$im->colorAllocate(0,0,0); my \$origin_color = \$im->colorAllocate(255,0,0); my \$position=\$turning_points; my (\$x, \$y) = (0,0); my (\$delx,\$dely) = (1,0); # start out moving to the right sub gdx { return( shift()+(\$x_width/2) ) } sub gdy { return( (\$y_width/2)-shift()) } \$im->setPixel(gdx(\$x),gdy(\$y),\$origin_color); foreach my \$next_turn (@turning_points) { for (my \$i=\$position;\$isetPixel(gdx(\$x),gdy(\$y),\$black); } } \$position = \$next_turn+1; } open(PNGFH, ">".\$pngfile) or die "Could not open image file \"\$pngfile\"."; print PNGFH \$im->png; close PNGFH; # an array of run lengths up to each turn # equation was guessed at after inspection of coordinates # from a hand drawn cartesian spiral sub mkturns{ my \$ulim = shift; my \$len = 0; my @rv; foreach my \$i (1..\$ulim) { my \$delta = (\$i%2?\$i+1:\$i)/2; \$len = \$len+\$delta; push( @rv, \$len); } return @rv; } # generates prime numbers as hash keys sub genprimes { my \$max = shift; my \$sieve = ''; my %primes = ( 2 => 1 ); GUESS: for (my \$guess = 2; \$guess <= \$max; \$guess++) { next GUESS if vec(\$sieve,\$guess,1); \$primes{\$guess} = 1; for (my \$mults = \$guess * \$guess; \$mults <= \$max; \$mults += \$guess) { vec(\$sieve,\$mults,1) = 1; } } return \%primes; } # uses memoized trig function # function modded to go past 180 degrees { my %memo = (); sub nextDirection { my (\$x, \$y) = (shift, shift); my \$nextoffset = undef; if( defined(\$memo{\$x}{\$y})){ \$nextoffset = \$memo{\$x}{\$y}; } else { my \$acos = acos(\$x); my \$asin = asin(\$y); if( \$asin \$asin?\$acos:\$asin)/PI); # hard coded count of four directions in spiral (right, up, left down) \$nextoffset = \$memo{\$x}{\$y} = (\$offset+1)%4; } return to_xy(\$nextoffset); } } # uses memoized trig function { my %memo = (); sub to_xy { defined(\$memo{\$_}) and return @{\$memo{\$_}}; \$memo{\$_}= [int(cos(\$_*PI/2)),int(sin(\$_*PI/2))]; return @{\$memo{\$_}}; } } ```