Allan’s Musings

Just another WordPress.com weblog

Posts Tagged ‘Ulam spiral’

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.

Prime numbers illustrated as Ulam spiral

Prime numbers Plotted as an Ulam Spiral


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

Posted in Mathematical Recreations | Tagged: , , | Leave a Comment »