File Coverage

blib/lib/Math/PlanePath/Hypot.pm

Criterion	Covered	Total	%
statement	121	186	65.0
branch	19	50	38.0
condition	5	17	29.4
subroutine	11	18	61.1
pod	9	9	100.0
total	165	280	58.9

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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 .



# A000328 Number of points of norm <= n^2 in square lattice.
#   1, 5, 13, 29, 49, 81, 113, 149, 197, 253, 317, 377, 441, 529, 613, 709, 797
#   a(n) = 1 + 4 * sum(j=0, n^2 / 4,    n^2 / (4*j+1) - n^2 / (4*j+3) )
# A014200 num points norm <= n^2, excluding 0, divided by 4
#
# A046109 num points norm == n^2
#
# A057655 num points x^2+y^2 <= n
# A014198 = A057655 - 1
#
# A004018 num points x^2+y^2 == n
#
# A057962 hypot count x-1/2,y-1/2 <= n
# is last point of each hypot in points=odd
#
# A057961 hypot count as radius increases
#

# points="square_horiz"
# points="square_vert"
# points="square_centre"
# A199015 square_centred partial sums
#


package Math::PlanePath::Hypot;
use 5.004;
use strict;
use Carp 'croak';

use vars '$VERSION', '@ISA';
$VERSION = 128;
use Math::PlanePath;
@ISA = ('Math::PlanePath');

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

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


use constant parameter_info_array =>
  [ { name            => 'points',
      share_key       => 'points_aeo',
      display         => 'Points',
      type            => 'enum',
      default         => 'all',
      choices         => ['all','even','odd'],
      choices_display => ['All','Even','Odd'],
      description     => 'Which X,Y points visit, either all of them or just X+Y=even or odd.',
    },
    Math::PlanePath::Base::Generic::parameter_info_nstart1(),
  ];

{
  my %x_negative_at_n = (all  => 3,
                         even => 2,
                         odd  => 2);
  sub x_negative_at_n {
    my ($self) = @_;
    return $self->n_start + $x_negative_at_n{$self->{'points'}};
  }
}
{
  my %y_negative_at_n = (all  => 4,
                                        even => 3,
                                        odd  => 3);
  sub y_negative_at_n {
    my ($self) = @_;
    return $self->n_start + $y_negative_at_n{$self->{'points'}};
  }
}
sub rsquared_minimum {
  my ($self) = @_;
  return ($self->{'points'} eq 'odd'
          ? 1     # odd at X=1,Y=0
          : 0);   # even,all at X=0,Y=0
}
# points=even includes X=Y so abs(X-Y)>=0
# points=odd doesn't include X=Y so abs(X-Y)>=1
*absdiffxy_minimum = \&rsquared_minimum;
*sumabsxy_minimum  = \&rsquared_minimum;

use constant turn_any_right => 0; # always left or straight
sub turn_any_straight {
  my ($self) = @_;
  return ($self->{'points'} ne 'all');  # points=all is left always
}


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

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

  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }

  my $points = ($self->{'points'} ||= 'all');
  if ($points eq 'all') {
    $self->{'n_to_x'} = [0];
    $self->{'n_to_y'} = [0];
    $self->{'hypot_to_n'} = [0];
    $self->{'y_next_x'} = [1, 1];
    $self->{'y_next_hypot'} = [1, 2];
    $self->{'x_inc'} = 1;
    $self->{'x_inc_factor'} = 2;
    $self->{'x_inc_squared'} = 1;
    $self->{'y_factor'} = 2;
    $self->{'opposite_parity'} = -1;

  } elsif ($points eq 'even') {
    $self->{'n_to_x'} = [0];
    $self->{'n_to_y'} = [0];
    $self->{'hypot_to_n'} = [0];
    $self->{'y_next_x'} = [2, 1];
    $self->{'y_next_hypot'} = [4, 2];
    $self->{'x_inc'} = 2;
    $self->{'x_inc_factor'} = 4;
    $self->{'x_inc_squared'} = 4;
    $self->{'y_factor'} = 2;
    $self->{'opposite_parity'} = 1;

  } elsif ($points eq 'odd') {
    $self->{'n_to_x'} = [];
    $self->{'n_to_y'} = [];
    $self->{'hypot_to_n'} = [];
    $self->{'y_next_x'} = [1];
    $self->{'y_next_hypot'} = [1];
    $self->{'x_inc'} = 2;
    $self->{'x_inc_factor'} = 4;
    $self->{'x_inc_squared'} = 4;
    $self->{'y_factor'} = 2;
    $self->{'opposite_parity'} = 0;

  } elsif ($points eq 'square_centred') {
    $self->{'n_to_x'} = [];
    $self->{'n_to_y'} = [];
    $self->{'hypot_to_n'} = [];
    $self->{'y_next_x'} = [undef,1];
    $self->{'y_next_hypot'} = [undef,2];
    $self->{'x_inc'} = 2;
    $self->{'x_inc_factor'} = 4;  # ((x+2)^2 - x^2) = 4*x+4
    $self->{'x_inc_squared'} = 4;
    $self->{'y_start'} = 1;
    $self->{'y_inc'} = 2;
    $self->{'opposite_parity'} = -1;

  } else {
    croak "Unrecognised points option: ", $points;
  }
  return $self;
}

sub _extend {
  my ($self) = @_;
  ### _extend() n: scalar(@{$self->{'n_to_x'}})
  ### y_next_x: $self->{'y_next_x'}

  my $n_to_x       = $self->{'n_to_x'};
  my $n_to_y       = $self->{'n_to_y'};
  my $hypot_to_n   = $self->{'hypot_to_n'};
  my $y_next_x     = $self->{'y_next_x'};
  my $y_next_hypot = $self->{'y_next_hypot'};
  my $y_start      = $self->{'y_start'} || 0;
  my $y_inc        = $self->{'y_inc'} || 1;

  # set @y to the Y with the smallest $y_next_hypot[$y], and if there's some
  # Y's with equal smallest hypot then all those Y's
  my @y = ($y_start);
  my $hypot = $y_next_hypot->[$y_start] || 99;
  for (my $y = $y_start+$y_inc; $y < @$y_next_x; $y += $y_inc) {
    if ($hypot == $y_next_hypot->[$y]) {
      push @y, $y;
    } elsif ($hypot > $y_next_hypot->[$y]) {
      @y = ($y);
      $hypot = $y_next_hypot->[$y];
    }
  }

  ### chosen y list: @y

  # if the endmost of the @$y_next_x, @$y_next_hypot arrays are used then
  # extend them by one
  if ($y[-1] == $#$y_next_x) {
    ### grow y_next_x ...
    my $y = $#$y_next_x + $y_inc;
    my $x = $y + ($self->{'points'} eq 'odd');
    $y_next_x->[$y] = $x;
    $y_next_hypot->[$y] = $x*$x+$y*$y;
    ### $y_next_x
    ### $y_next_hypot
    ### assert: $y_next_hypot->[$y] == $y**2 + $x*$x
  }

  # @x is the $y_next_x[$y] for each of the @y smallests, and step those
  # selected elements next X and hypot for that new X,Y
  my @x = map {
    my $y = $_;
    my $x = $y_next_x->[$y];
    $y_next_x->[$y] += $self->{'x_inc'};
    $y_next_hypot->[$y]
      += $self->{'x_inc_factor'} * $x + $self->{'x_inc_squared'};
    ### assert: $y_next_hypot->[$y] == ($x+$self->{'x_inc'})**2 + $y**2
    $x
  } @y;
  ### $hypot
  ### base octant: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)

  # transpose X,Y to Y,X
  {
    my @base_x = @x;
    my @base_y = @y;
    unless ($y[0]) { # no transpose of x,0
      shift @base_x;
      shift @base_y;
    }
    if ($x[-1] == $y[-1]) { # no transpose of x,x
      pop @base_x;
      pop @base_y;
    }
    push @x, reverse @base_y;
    push @y, reverse @base_x;
  }
  ### with transpose q1: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)

  # rotate +90 quadrant 1 into quadrant 2
  {
    my @base_y = @y;
    push @y, @x;
    push @x, map {-$_} @base_y;
  }
  ### with rotate q2: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)

  # rotate +180 quadrants 1+2 into quadrants 2+3
  push @x, map {-$_} @x;
  push @y, map {-$_} @y;

  ### store: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)
  ### at n: scalar(@$n_to_x)
  ### hypot_to_n: "h=$hypot n=".scalar(@$n_to_x)
  $hypot_to_n->[$hypot] = scalar(@$n_to_x);
  push @$n_to_x, @x;
  push @$n_to_y, @y;

  # ### hypot_to_n now: join(' ',map {defined($hypot_to_n->[$_]) && "h=$_,n=$hypot_to_n->[$_]"} 0 .. $#$hypot_to_n)


  # my $x = $y_next_x->[0];
  #
  # $x = $y_next_x->[$y];
  # $n_to_x->[$next_n] = $x;
  # $n_to_y->[$next_n] = $y;
  # $xy_to_n{"$x,$y"} = $next_n++;
  #
  # $y_next_x->[$y]++;
  # $y_next_hypot->[$y] = $y*$y + $y_next_x->[$y]**2;
}

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

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

  my $int = int($n);
  $n -= $int;  # fraction part

  my $n_to_x = $self->{'n_to_x'};
  my $n_to_y = $self->{'n_to_y'};

  while ($int >= $#$n_to_x) {
    _extend($self);
  }

  my $x = $n_to_x->[$int];
  my $y = $n_to_y->[$int];
  return ($x + $n * ($n_to_x->[$int+1] - $x),
          $y + $n * ($n_to_y->[$int+1] - $y));
}

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

  if ($self->{'opposite_parity'} >= 0) {
    $x = round_nearest ($x);
    $y = round_nearest ($y);
    if ((($x%2) ^ ($y%2)) == $self->{'opposite_parity'}) {
      return 0;
    }
  }
  if ($self->{'points'} eq 'square_centred') {
    unless (($y%2) && ($x%2)) {
      return 0;
    }
  }
  return 1;
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### Hypot xy_to_n(): "$x, $y"
  ### hypot_to_n last: $#{$self->{'hypot_to_n'}}

  $x = round_nearest ($x);
  $y = round_nearest ($y);

  if ((($x%2) ^ ($y%2)) == $self->{'opposite_parity'}) {
    return undef;
  }
  if ($self->{'points'} eq 'square_centred') {
    unless (($y%2) && ($x%2)) {
      return undef;
    }
  }

  my $hypot = $x*$x + $y*$y;
  if (is_infinite($hypot)) {
    ### infinity
    return undef;
  }

  my $n_to_x = $self->{'n_to_x'};
  my $n_to_y = $self->{'n_to_y'};

  my $hypot_to_n = $self->{'hypot_to_n'};
  while ($hypot > $#$hypot_to_n) {
    _extend($self);
  }

  my $n = $hypot_to_n->[$hypot];
  for (;;) {
    if ($x == $n_to_x->[$n] && $y == $n_to_y->[$n]) {
      return $n + $self->{'n_start'};
    }
    $n += 1;

    if ($n_to_x->[$n]**2 + $n_to_y->[$n]**2 != $hypot) {
      ### oops, hypot_to_n no good ...
      return undef;
    }
  }

  # if ($x < 0 || $y < 0) {
  #   return undef;
  # }
  # my $h = $x*$x + $y*$y;
  #
  # while ($y_next_x[$y] <= $x) {
  #   _extend($self);
  # }
  # return $xy_to_n{"$x,$y"};
}

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

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

  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }

  # circle area pi*r^2, with r^2 = $x2**2 + $y2**2
  return ($self->{'n_start'},
          $self->{'n_start'} + int (3.2 * (($x2+1)**2 + ($y2+1)**2)));
}

1;
__END__