## 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.

A color pdf version of the source is available here, and a downloadable version should be here

```
```#!/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[0])){
die "Please provide a png file to save to.";
}
my $pngfile = 'ulam_spiral.png';
if($ARGV[0] =~ m/\.png$/i){
$pngfile = $ARGV[0];
} 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[0];
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{$_[0]}) and return @{$memo{$_[0]}};
$memo{$_[0]}= [int(cos($_[0]*PI/2)),int(sin($_[0]*PI/2))];
return @{$memo{$_[0]}};
}
}

Advertisements

Sorry, the comment form is closed at this time.