File Coverage

blib/lib/Math/PlanePath/SquareSpiral.pm

Criterion	Covered	Total	%
statement	88	114	77.1
branch	27	36	75.0
condition	2	2	100.0
subroutine	13	19	68.4
pod	7	7	100.0
total	137	178	76.9

# Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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://d4maths.lowtech.org/mirage/ulam.htm
# http://d4maths.lowtech.org/mirage/img/ulam.gif
#     sample gif of primes made by APL or something
#
# http://www.sciencenews.org/view/generic/id/2696/title/Prime_Spirals
#     Ulam's spiral of primes
#
# http://web.archive.org/web/20160424015805id_/http://yoyo.cc.monash.edu.au/%7Ebunyip/primes/primeSpiral.htm
# http://yoyo.cc.monash.edu.au/%7Ebunyip/primes/triangleUlam.htm
#     Pulchritudinous Primes of Ulam spiral.

# http://mathworld.wolfram.com/PrimeSpiral.html
#
# Mark C. Chu-Carroll "The Surprises Never Eend: The Ulam Spiral of Primes"
# http://scienceblogs.com/goodmath/2010/06/the_surprises_never_eend_the_u.php
#
# http://yoyo.cc.monash.edu.au/%7Ebunyip/primes/index.html
# including image highlighting the lines

# S. M. Ellerstein, The Square Spiral, Journal of Recreational
# Mathematics 29 (#3, 1998) 188; 30 (#4, 1999-2000), 246-250.
#
# M. Stein and S. M. Ulam. "An Observation on the
# Distribution of Primes." Amer. Math. Monthly 74, 43-44,
# 1967.
#
# M.L. Stein; S. M. Ulam; and M. B. Wells. "A Visual
# Display of Some Properties of the Distribution of Primes."
# Amer. Math. Monthly 71, 516-520, 1964.

# cf sides alternately prime and Fibonacci
# A160790 corner N
# A160791 side lengths, alternately integer and triangular adding that integer
# A160792 corner N
# A160793 side lengths, alternately integer and sum primes
# A160794 corner N
# A160795 side lengths, alternately primes and fibonaccis


package Math::PlanePath::SquareSpiral;
use 5.004;
use strict;
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

use vars '$VERSION', '@ISA';
$VERSION = 128;
use Math::PlanePath;
*_sqrtint = \&Math::PlanePath::_sqrtint;
use Math::PlanePath::Base::NSEW;
@ISA = ('Math::PlanePath::Base::NSEW',
        'Math::PlanePath');

use Math::PlanePath::Base::Generic
  'round_nearest';

# uncomment this to run the ### lines
#use Smart::Comments '###';


# Note: this shared by other paths
use constant parameter_info_array =>
  [
   { name        => 'wider',
     display     => 'Wider',
     type        => 'integer',
     minimum     => 0,
     default     => 0,
     width       => 3,
     description => 'Wider path.',
   },
   Math::PlanePath::Base::Generic::parameter_info_nstart1(),
  ];

use constant xy_is_visited => 1;

#    2w+4 -- 2w+3 ----- w+2
#      |                 |
#    2w+5      0------- w+1
#      |
#    2w+6 ---
#                  ^
#                 X=0
#
sub x_negative_at_n {
  my ($self) = @_;
  return $self->n_start + ($self->{'wider'} ? 0 : 4);
}
sub y_negative_at_n {
  my ($self) = @_;
  return $self->n_start + 2*$self->{'wider'} + 6;
}
sub _UNDOCUMENTED__dxdy_list_at_n {
  my ($self) = @_;
  return $self->n_start + 2*$self->{'wider'} + 4;
}

use constant turn_any_right => 0; # only left or straight

sub _UNDOCUMENTED__turn_any_left_at_n {
  my ($self) = @_;
  return $self->n_start + $self->{'wider'} + 1;
}


#------------------------------------------------------------------------------

sub new {
  my $self = shift->SUPER::new (@_);

  # parameters
  $self->{'wider'} ||= 0;  # default
  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }

  return $self;
}

# wider==0
# base from bottom-right corner
#   d = [ 1,  2,  3,  4 ]
#   N = [ 2, 10, 26, 50 ]
#   N = (4 d^2 - 4 d + 2)
#   d = 1/2 + sqrt(1/4 * $n + -4/16)
#
# wider==1
# base from bottom-right corner
#   d = [ 1,  2,  3,  4 ]
#   N = [ 3, 13, 31, 57 ]
#   N = (4 d^2 - 2 d + 1)
#   d = 1/4 + sqrt(1/4 * $n + -3/16)
#
# wider==2
# base from bottom-right corner
#   d = [ 1,  2,  3, 4 ]
#   N = [ 4, 16, 36, 64 ]
#   N = (4 d^2)
#   d = 0 + sqrt(1/4 * $n + 0)
#
# wider==3
# base from bottom-right corner
#   d = [ 1,  2,  3 ]
#   N = [ 5, 19, 41 ]
#   N = (4 d^2 + 2 d - 1)
#   d = -1/4 + sqrt(1/4 * $n + 5/16)
#
# N = 4*d^2 + (-4+2*w)*d + (2-w)
#   = 4*$d*$d + (-4+2*$w)*$d + (2-$w)
# d = 1/2-w/4 + sqrt(1/4*$n + b^2-4ac)
# (b^2-4ac)/(2a)^2 = [ (2w-4)^2 - 4*4*(2-w) ] / 64
#                  = [ 4w^2 - 16w + 16 - 32 + 16w ] / 64
#                  = [ 4w^2 - 16 ] / 64
#                  = [ w^2 - 4 ] / 16
# d = 1/2-w/4 + sqrt(1/4*$n + (w^2 - 4) / 16)
#   = 1/4 * (2-w + sqrt(4*$n + w^2 - 4))
#   = 0.25 * (2-$w + sqrt(4*$n + $w*$w - 4))
#
# then offset the base by +4*$d+$w-1 for top left corner for +/- remainder
# rem = $n - (4*$d*$d + (-4+2*$w)*$d + (2-$w) + 4*$d + $w - 1)
#     = $n - (4*$d*$d + (-4+2*$w)*$d + 2 - $w + 4*$d + $w - 1)
#     = $n - (4*$d*$d + (-4+2*$w)*$d + 1 - $w + 4*$d + $w)
#     = $n - (4*$d*$d + (-4+2*$w)*$d + 1 + 4*$d)
#     = $n - (4*$d*$d + (2*$w)*$d + 1)
#     = $n - ((4*$d + 2*$w)*$d + 1)
#

sub n_to_xy {
  my ($self, $n) = @_;
  #### SquareSpiral n_to_xy: $n

  $n = $n - $self->{'n_start'};  # starting $n==0, warn if $n==undef
  if ($n < 0) {
    #### before n_start ...
    return;
  }

  my $w = $self->{'wider'};
  my $w_right = int($w/2);
  my $w_left = $w - $w_right;
  if ($n <= $w+1) {
    #### centre horizontal
    # n=0 at w_left
    # x = $n - int(($w+1)/2)
    #   = $n - int(($w+1)/2)
    return ($n - $w_left,  # n=0 at w_left
            0);
  }

  my $d = int ((2-$w + _sqrtint(4*$n + $w*$w)) / 4);
  #### d frac: ((2-$w + sqrt(int(4*$n) + $w*$w)) / 4)
  #### $d

  #### base: 4*$d*$d + (-4+2*$w)*$d + (2-$w)
  $n -= ((4*$d + 2*$w)*$d);
  #### remainder: $n

  if ($n >= 0) {
    if ($n <= 2*$d) {
      ### left vertical
      return (-$d - $w_left,
              -$n + $d);
    } else {
      ### bottom horizontal
      return ($n - $w_left - 3*$d,
              -$d);
    }
  } else {
    if ($n >= -2*$d-$w) {
      ### top horizontal
      return (-$n - $d - $w_left,
              $d);
    } else {
      ### right vertical
      return ($d + $w_right,
              $n + 3*$d + $w);
    }
  }
}

sub xy_to_n {
  my ($self, $x, $y) = @_;

  my $w = $self->{'wider'};
  my $w_right = int($w/2);
  my $w_left = $w - $w_right;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### xy_to_n: "x=$x, y=$y"
  ### $w_left
  ### $w_right

  my $d;
  if (($d = $x - $w_right) > abs($y)) {
    ### right vertical
    ### $d
    #
    # base bottom right per above
    ### BR: 4*$d*$d + (-4+2*$w)*$d + (2-$w)
    # then +$d-1 for the y=0 point
    # N_Y0  = 4*$d*$d + (-4+2*$w)*$d + (2-$w) + $d-1
    #       = 4*$d*$d + (-3+2*$w)*$d + (2-$w) + -1
    #       = 4*$d*$d + (-3+2*$w)*$d +  1-$w
    ### N_Y0: (4*$d + -3 + 2*$w)*$d + 1-$w
    #
    return (4*$d + -3 + 2*$w)*$d - $w + $y + $self->{'n_start'};
  }

  if (($d = -$x - $w_left) > abs($y)) {
    ### left vertical
    ### $d
    #
    # top left per above
    ### TL: 4*$d*$d + (2*$w)*$d + 1
    # then +$d for the y=0 point
    # N_Y0  = 4*$d*$d + (2*$w)*$d + 1 + $d
    #       = 4*$d*$d + (1 + 2*$w)*$d + 1
    ### N_Y0: (4*$d + 1 + 2*$w)*$d + 1
    #
    return (4*$d + 1 + 2*$w)*$d - $y + $self->{'n_start'};
  }

  $d = abs($y);
  if ($y > 0) {
    ### top horizontal
    ### $d
    #
    # top left per above
    ### TL: 4*$d*$d + (2*$w)*$d + 1
    # then -($d+$w_left) for the x=0 point
    # N_X0  = 4*$d*$d + (2*$w)*$d + 1 + -($d+$w_left)
    #       = 4*$d*$d + (-1 + 2*$w)*$d + 1 - $w_left
    ### N_Y0: (4*$d - 1 + 2*$w)*$d + 1 - $w_left
    #
    return (4*$d - 1 + 2*$w)*$d - $w_left - $x + $self->{'n_start'};
  }

  ### bottom horizontal, and centre y=0
  ### $d
  #
  # top left per above
  ### TL: 4*$d*$d + (2*$w)*$d + 1
  # then +2*$d to bottom left, +$d+$w_left for the x=0 point
  # N_X0  = 4*$d*$d + (2*$w)*$d + 1 + 2*$d + $d+$w_left)
  #       = 4*$d*$d + (3 + 2*$w)*$d + 1 + $w_left
  ### N_Y0: (4*$d + 3 + 2*$w)*$d + 1 + $w_left
  #
  return (4*$d + 3 + 2*$w)*$d + $w_left + $x + $self->{'n_start'};
}

# hi is exact but lo is not
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;

  $x1 = round_nearest ($x1);
  $y1 = round_nearest ($y1);
  $x2 = round_nearest ($x2);
  $y2 = round_nearest ($y2);

  # ENHANCE-ME: find actual minimum if rect doesn't cover 0,0
  return ($self->{'n_start'},
          max ($self->xy_to_n($x1,$y1),
               $self->xy_to_n($x2,$y1),
               $self->xy_to_n($x1,$y2),
               $self->xy_to_n($x2,$y2)));

  # my $w = $self->{'wider'};
  # my $w_right = int($w/2);
  # my $w_left = $w - $w_right;
  #
  # my $d = 1 + max (abs($y1),
  #                  abs($y2),
  #                  $x1 - $w_right, -$x1 - $w_left,
  #                  $x2 - $w_right, -$x2 - $w_left,
  #                  1);
  # ### $d
  # ### is: $d*$d
  #
  # # ENHANCE-ME: find actual minimum if rect doesn't cover 0,0
  # return (1,
  #         (4*$d - 4 + 2*$w)*$d + 2);  # bottom-right
}


# [ 1, 2, 3,  4,  5 ],
# [ 1, 3, 7, 13, 21 ]
# N = (d^2 - d + 1)
#   = ($d**2 - $d + 1)
#   = (($d - 1)*$d + 1)
# d = 1/2 + sqrt(1 * $n + -3/4)
#   = (1 + sqrt(4*$n - 3)) / 2
#
# wider=3
# [ 2, 3,  4,  5 ],
# [ 6, 13, 22, 33 ]
# N = (d^2 + 2 d - 2)
#   = ($d**2 + 2*$d - 2)
#   = (($d + 2)*$d - 2)
# d = -1 + sqrt(1 * $n + 3)
#
# wider=5
# [ 2, 3,  4,  5 ],
# [ 8, 17, 28, 41 ]
# N = (d^2 + 4 d - 4)
#   = ($d**2 + 4*$d - 4)
#   = (($d + 4)*$d - 4)
# d = -2 + sqrt(1 * $n + 8)
#
# wider=7
# [ 2, 3,  4,  5 ],
# [ 10, 21, 34, 49 ]
# N = (d^2 + 6 d - 6)
#   = ($d**2 + 6*$d - 6)
#   = (($d + 6)*$d - 6)
# d = -3 + sqrt(1 * $n + 15)
#
#
# N = (d^2 + (w-1)*d + 1-w)
# d = (1-w)/2 + sqrt($n + (w^2 + 2w - 3)/4)
#   = (1-w + sqrt(4*$n + (w-3)(w+1))) / 2
#
# extra subtract d+w-1
# Nbase = (d^2 + (w-1)*d + 1-w) + d+w-1
#       = d^2 + w*d

sub n_to_dxdy {
  my ($self, $n) = @_;
  ### n_to_dxdy(): $n

  $n = $n - $self->{'n_start'};  # starting $n==0, warn if $n==undef
  if ($n < 0) {
    #### before n_start ...
    return;
  }

  my $w = $self->{'wider'};
  my $d = int((1-$w + _sqrtint(4*$n + ($w+2)*$w+1)) / 2);

  my $int = int($n);
  $n -= $int;  # fraction 0 <= $n < 1
  $int -= ($d+$w)*$d-1;

  ### $d
  ### $w
  ### $n
  ### $int

  my ($dx, $dy);
  if ($int <= 0) {
    if ($int < 0) {
      ### horizontal ...
      $dx = 1;
      $dy = 0;
    } else {
      ### corner horiz to vert ...
      $dx = 1-$n;
      $dy = $n;
    }
  } else {
    if ($int < $d) {
      ### vertical ...
      $dx = 0;
      $dy = 1;
    } else {
      ### corner vert to horiz ...
      $dx = -$n;
      $dy = 1-$n;
    }
  }

  unless ($d % 2) {
    ### rotate +180 for even d ...
    $dx = -$dx;
    $dy = -$dy;
  }

  ### result: "$dx, $dy"
  return ($dx,$dy);
}



# old bit:
#
# wider==0
# base from two-way diagonal top-right and bottom-left
# s even for top-right diagonal doing top leftwards then left downwards
# s odd for bottom-left diagonal doing bottom rightwards then right pupwards
#   s = [ 0,  1,   2,   3,   4,   5,   6 ]
#   N = [ 1,  1,   3,   7,  13,  21,  31 ]
#         +0  +2  +4  +6  +8  +10
#            2   2   2   2   2
#
#   n = (($d - 1)*$d + 1)
#   s = 1/2 + sqrt(1 * $n + -3/4)
#     = .5 + sqrt ($n - .75)
#
#

#------------------------------------------------------------------------------

sub _NOTDOCUMENTED_n_to_figure_boundary {
  my ($self, $n) = @_;
  ### _NOTDOCUMENTED_n_to_figure_boundary(): $n

  # adjust to N=1 at origin X=0,Y=0
  $n = $n - $self->{'n_start'} + 1;

  if ($n < 1) {
    return undef;
  }

  my $wider = $self->{'wider'};
  if ($n <= $wider) {
    # single block row
    # +---+-----+----+
    # | 1 | ... | $n |  boundary = 2*N + 2
    # +---+-----+----+
    return 2*$n + 2;
  }

  my $d = int((_sqrtint(4*$n + $wider*$wider - 2) - $wider) / 2);
  ### $d
  ### $wider
  ### cmp: $d*($d+1+$wider) + $wider + 1

  if ($n > $d*($d+1+$wider)) {
    $wider++;
    ### increment for +2 after turn ...
  }
  return 4*$d + 2*$wider + 2;
}

#------------------------------------------------------------------------------
1;
__END__