#!/usr/bin/perl -w

# Copyright 2010, 2011, 2012 Kevin Ryde

# This file is part of Math-PlanePath.
#
# Math-PlanePath is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-PlanePath is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.


# Usage: perl ulam-spiral-xpm.pl >/tmp/foo.xpm     # write image file
#        xzgv /tmp/foo.xpm                         # view file
#
# This is a bit of fun drawing Ulam's spiral of primes in the SquareSpiral
# path.  The output is XPM format (which is plain text) and any good image
# viewer program should display it.
#
# Optional args
#
#     perl ulam-spiral-xpm.pl SIZE
# or
#     perl ulam-spiral-xpm.pl SIZE SCALE
#
# make the image SIZExSIZE pixels, and SCALE to expand each point to a
# SCALExSCALE square instead of a single pixel.
#

use 5.004;
use strict;
use Math::PlanePath::SquareSpiral;

my $size = 200;
my $scale = 1;

if (@ARGV >= 2) {
  $scale = $ARGV[1];
}
if (@ARGV >= 1) {
  $size = $ARGV[0];
}

my $path = Math::PlanePath::SquareSpiral->new;
my $x_origin = int($size / 2);
my $y_origin = int($size / 2);

my ($n_lo, $n_hi)
  = $path->rect_to_n_range (-$x_origin, -$y_origin,
                            -$x_origin+$size, -$y_origin+$size);

# Find the prime numbers 2 to $n_hi by sieve of Eratosthenes.
# Could also use Math::Prime::TiedArray or Math::Prime::XS.
#
my @primes = (0,    # 0
              0,    # 1
              1,    # 2  prime
              1,    # 3  prime
              (0,1) x ($n_hi/2));  # rest alternately even/odd
my $i = 3;
foreach my $i (3 .. int(sqrt($n_hi)) + 1) {
  next unless $primes[$i];
  foreach (my $j = 2*$i; $j <= $n_hi; $j += $i) {
    $primes[$j] = 0;
  }
}

# Draw the primes into an array of rows strings.
#
my @rows = (' ' x $size) x $size;

foreach my $n ($n_lo .. $n_hi) {
  next unless $primes[$n];

  my ($x, $y) = $path->n_to_xy ($n);

  $x = $x + $x_origin;
  $y = $y_origin - $y;  # inverted

  # $n_hi is an over-estimate in general, check x,y actually in desired size
  if ($x >= 0 && $x < $size && $y >= 0 && $y < $size) {
    substr ($rows[$y], $x,1) = '*';
  }
}

# Expand @rows points by $scale, horizontally and vertically.
#
if ($scale > 1) {
  foreach (@rows) {
    s{(.)}{$1 x $scale}eg;             # expand horizontally
  }
  @rows = map { ($_) x $scale} @rows;  # expand vertically

  $size *= $scale;
}

# XPM format is easy to print.
# Output is about 1 byte per pixel.
#
print <<"HERE";
/* XPM */
static char *ulam_spiral_xpm_pl[] = {
"$size $size 2 1",
" 	c black",
"*	c white",
HERE
foreach my $row (@rows) {
  print "\"$row\",\n";
}
print "};\n";

exit 0;