libmath-planepath-perl 113-1 / usr / share / perl5 / Math / PlanePath / FilledRings.pm

# Copyright 2012, 2013 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/>.


# math-image --path=FilledRings --all --output=numbers_dash  --size=70x30
#
# offset 0 to 1
# same PixelRings fractional offset for midpoint
#
# default 0.5     |-------int-------|             0.5
#         |------------------|                    0
#                         |------------------|   0.99 or 1.0
#
# default  0 |-------int-------|--------int-------|
#                   |------------------|             0.5
#                             |------------------|   0.99 or 1.0
#
# innermost
#   |-------0-------|-------1--------|
#
#
# offset -0.5 to +0.5
# h-0.5 < R < h+0.5
# h-0.5 <= R-0.5 < h+0.5
# h-0.5 < R+0.5 <= h+0.5

# A036702 count Gaussian |z| <= n
# A036706 count Gaussian n-1/2 < |z| < n+1/2 with a>0,b>=0, so 1/4
# A036707 count Gaussian |z| < n+1/2 with b>=0, so 1/2 plane
# A036708 count Gaussian n-1/2 < |z| < n+1/2 with b>=0, so 1/4
#
# A000328 num points <= circle radius n
#         1, 5, 13, 29, 49, 81, 113, 149, 197, 253, 317, 377, 441,
# A046109 num points == circle radius n
# A051132 num points <  circle radius n
#         0, 1, 9, 25, 45, 69, 109, 145, 193, 249, 305, 373, 437,
# A057655 num points x^2+y^2 <= n
#         1, 5, 9, 9, 13, 21, 21, 21, 25, 29, 37, 37, 37, 45, 45,
#

package Math::PlanePath::FilledRings;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 113;
use Math::PlanePath;
@ISA = ('Math::PlanePath');
*_divrem = \&Math::PlanePath::_divrem;

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

use Math::PlanePath::SacksSpiral;
*_rect_to_radius_corners = \&Math::PlanePath::SacksSpiral::_rect_to_radius_corners;

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


# N(r) = 1 + 4*sum  floor(r^2/(4i+1)) - floor(r^2/(4i+3))
#
# N(r+1) - N(r)
#   = 1 + 4*sum  floor((r+1)^2/(4i+1)) - floor((r+1)^2/(4i+3))
#     - 1 + 4*sum  floor(r^2/(4i+1)) - floor(r^2/(4i+3))
#   = 4*sum  floor(((r+1)^2-r^2)/(4i+1)) - floor(((r+1)^2-r^2)/(4i+3))
#   = 4*sum  floor((2r+1)/(4i+1)) - floor((2r+1)/(4i+3))
#
# _cumul[0] index=0 is r=1/2
#  r = index+1/2
#  2r+1 = 2(index+1/2)+1
#       = 2*index+1+1
#       = 2*index+2
#
#  2r+1 >= 4i+1
#  2r >= 4i
#  i <= (2*index+2)/2
#  i <= index+1
#
#  r=3.5
#  sqrt(3*3+3*3) = 4.24 out
#  sqrt(3*3+2*2) = 3.60 out
#  sqrt(3*3+1*1) = 3.16 in
#
#      * * *
#    * * * * *
#  * * * * * * *
#  * * * o * * *   3+5+7+7+7+5+3 = 37
#  * * * * * * *
#    * * * * *
#      * * *
#
# N(r) = 1 + 4*( floor(12.25/1)-floor(12.25/3)
#          + floor(12.25/5)-floor(12.25/7)
#          + floor(12.25/9)-floor(12.25/11) )
#      = 37
#
# (index+1/2)^2 = index^2 + index + 1/4
#               >= index*(index+1)
# (end+1 + 1/2)^2
#   = (end+3/2)^2
#   = end^2 + 3*end + 9/4
#   = end*(end+3) + 2 + 1/4
#
# (r+1/2)^2 = r^2+r+1/4  floor=r*(r+1)
# (r-1/2)^2 = r^2-r+1/4  ceil=r*(r-1)+1

use constant n_frac_discontinuity => 0;
use constant xy_is_visited => 1;

use constant dx_minimum => -1;
use constant dx_maximum => 1;
use constant dy_minimum => -1;
use constant dy_maximum => 1;
use constant dsumxy_minimum => -2; # diagonals
use constant dsumxy_maximum => 2;
use constant ddiffxy_minimum => -2;
use constant ddiffxy_maximum => 2;
use constant dir_maximum_dxdy => (1,-1); # South-East

use constant parameter_info_array =>
  [
   Math::PlanePath::Base::Generic::parameter_info_nstart1(),
  ];


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

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

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

  # internals
  $self->{'cumul'} = [ 1 ];  # N=0 basis
  $self->{'offset'} ||= 0;

  return $self;
}

sub _cumul_extend {
  my ($self) = @_;
  ### _cumul_extend() ...

  my $cumul = $self->{'cumul'};
  my $r2 = ($#$cumul + 3) * $#$cumul + 2;
  my $c = 0;
  for (my $d = 1; $d <= $r2; $d += 4) {
    $c += int($r2/$d) - int($r2/($d+2));
  }
  push @$cumul, 4*$c + 1;
  ### $cumul
}

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

  $n = $n - $self->{'n_start'};  # to N=0 basis, warning if $n==undef
  if ($n <= 1) {
    if ($n < 0) {
      return;
    } else {
      return ($n, 0);   # 0<=N<=1
    }
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }

  {
    # ENHANCE-ME: direction of N+1 from the cumulative lookup
    my $int = int($n);
    if ($n != $int) {
      my $frac = $n - $int;
      my ($x1,$y1) = $self->n_to_xy($int + $self->{'n_start'});
      my ($x2,$y2) = $self->n_to_xy($int+1 + $self->{'n_start'});
      if ($y2 == 0 && $x2 > 0) { $x2 -= 1; }
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int;
  }

  ### search cumul for: "n=$n"
  my $cumul = $self->{'cumul'};
  my $r = 1;
  for (;;) {
    if ($r > $#$cumul) {
      _cumul_extend ($self);
    }
    if ($cumul->[$r] > $n) {
      last;
    }
    $r++;
  }
  ### $r

  $n -= $cumul->[$r-1];
  my $len = $cumul->[$r] - $cumul->[$r-1];   # length of this ring

  ### cumul: "$cumul->[$r-1] to $cumul->[$r]"
  ### $len
  ### n rem: $n

  $len /= 4;     # length of a quadrant of this ring
  (my $quadrant, $n) = _divrem ($n, $len);

  ### len of quadrant: $len
  ### $quadrant
  ### n into quadrant: $n
  ### assert: $quadrant >= 0
  ### assert: $quadrant < 4

  my $rev;
  if ($rev = ($n > $len/2)) {
    $n = $len - $n;
  }
  ### $rev
  ### $n

  # my $rhi = ($r+1)*$r;
  # my $rlo = ($r-1)*$r+1;

  my $rlo = ($r-0.5+$self->{'offset'})**2;
  my $rhi = ($r+0.5+$self->{'offset'})**2;
  my $x = $r;
  my $y = 0;
  while ($n > 0) {
    ### at: "$x,$y n=$n"

    $y++;
    ### inc y to: $y

    if ($x*$x + $y*$y > $rhi) {
      $x--;
      ### dec x to: $x
      ### assert: $x*$x + $y*$y <= $rhi
      ### assert: $x*$x + $y*$y >= $rlo
    }
    $n--;
    last if $n <= 0;

    if (($x-1)*($x-1) + $y*$y >= $rlo) {
      ### another dec x to: $x
      $x--;
      $n--;
      last if $n <= 0;
    }
  }

  # if ($n) {
  #   ### n frac: $n
  # }

  if ($rev) {
    ($x,$y) = ($y,$x);
  }
  if ($quadrant & 2) {
    $x = -$x;
    $y = -$y;
  }
  if ($quadrant & 1) {
    ($x,$y) = (-$y, $x);
  }
  ### return: "$x, $y"
  return ($x, $y);
}


# h=x^2+y^2
# h >= (r-1/2)^2
# sqrt(h) >= r-1/2
# sqrt(h)+1/2 >= r
# r = int (sqrt(h)+1/2)
#   = int( (2*sqrt(h)+1)/2 }
#   = int( (sqrt(4*h) + 1)/2 }

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### FilledRings xy_to_n(): "$x, $y"
  $x = round_nearest ($x);
  $y = round_nearest ($y);

  if ($x == 0 && $y == 0) {
    return $self->{'n_start'};
  }

  my $r = int ((sqrt(4*($x*$x+$y*$y)) + 1) / 2);
  ### $r
  if (is_infinite($r)) {
    return undef;
  }

  my $cumul = $self->{'cumul'};
  while ($#$cumul < $r) {
    _cumul_extend ($self);
  }
  my $n = $cumul->[$r-1];
  ### n base: $n

  my $len = $cumul->[$r] - $n;
  ### $len
  $len /= 4;
  ### len/4: $len

  if ($y < 0) {
    ### y neg, rotate 180
    $y = -$y;
    $x = -$x;
    $n += 2*$len;
  }

  if ($x < 0) {
    $n += $len;
    ($x,$y) = ($y,-$x);
    ### neg x, rotate 90
    ### n base now: $n
  }

  ### assert: $x >= 0
  ### assert: $y >= 0

  my $rev;
  if ($rev = ($x < $y)) {
    ### top octant, reverse: "x=$x len/4=".($len/4)." gives ".($len/4 - $x)
    ($x,$y) = ($y,$x);
  }

  my $offset = 0;
  my $rhi = ($r+1)*$r;
  my $rlo = ($r-1)*$r+1;
  ### assert: $x*$x + $y*$y <= $rhi
  ### assert: $x*$x + $y*$y >= $rlo

  my $tx = $r;
  my $ty = 0;
  while ($ty < $y) {
    ### at: "$tx,$ty offset=$offset"

    $ty++;
    ### inc ty to: $ty
    if ($tx*$tx + $ty*$ty > $rhi) {
      $tx--;
      ### dec tx to: $tx
      ### assert: $tx*$tx + $ty*$ty <= $rhi
      ### assert: $tx*$tx + $ty*$ty >= $rlo
    }
    $offset++;
    last if $x == $tx && $y == $ty;

    if (($tx-1)*($tx-1) + $ty*$ty >= $rlo) {
      ### another dec tx to: "tx=$tx"
      $tx--;
      $offset++;
      last if $y == $ty;
    }
  }

  $n += $self->{'n_start'};
  if ($rev) {
    return $n + $len - $offset;
  } else {
    return $n + $offset;
  }
}

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

  ($x1,$y1, $x2,$y2) = _rect_to_radius_corners ($x1,$y1, $x2,$y2);
  ### radius range: "$x1,$y1 $x2,$y2"

  if ($x1 >= 1) { $x1 -= 1; }
  if ($y1 >= 1) { $y1 -= 1; }
  $x2 += 1;
  $y2 += 1;

  return (int((21*($x1*$x1 + $y1*$y1)) / 7) + $self->{'n_start'},
          int((22*($x2*$x2 + $y2*$y2)) / 7) + $self->{'n_start'} - 1);
}

1;
__END__

=for stopwords Ryde Math-PlanePath OEIS

=head1 NAME

Math::PlanePath::FilledRings -- concentric filled lattice rings

=head1 SYNOPSIS

 use Math::PlanePath::FilledRings;
 my $path = Math::PlanePath::FilledRings->new;
 my ($x, $y) = $path->n_to_xy (123);

=head1 DESCRIPTION

This path puts points on integer X,Y pixels of filled rings with radius 1
unit each ring.

                    110-109-108-107-106                        6
                   /                   \
            112-111  79--78--77--76--75 105-104                5
              |    /                   \      |
        114-113  80  48--47--46--45--44  74 103-102            4
          |    /      |               |    \      |
        115  81  50--49  27--26--25  43--42  73 101            3
       /   /      |    /           \      |    \   \
    116  82  52--51  28  14--13--12  24  41--40  72 100        2
      |   |   |    /   /           \   \      |   |   |
    117  83  53  29  15   5-- 4-- 3  11  23  39  71  99        1
      |   |   |   |   |   |       |   |   |   |   |   |
    118  84  54  30  16   6   1-- 2  10  22  38  70  98   <- Y=0
      |   |   |   |   |   |        /   /   /   /   /   /
    119  85  55  31  17   7-- 8-- 9  21  37  69  97 137       -1
      |   |   |    \   \           /   /      |   |   |
    120  86  56--57  32  18--19--20  36  67--68  96 136       -2
       \   \      |    \           /      |    /   /
        121  87  58--59  33--34--35  65--66  95 135           -3
          |    \      |               |    /      |
        122-123  88  60--61--62--63--64  94 133-134           -4
              |    \                   /      |
            124-125  89--90--91--92--93 131-132               -5
                   \                   /
                    126-127-128-129-130

                              ^
     -6  -5  -4  -3  -2  -1  X=0  1   2   3   4   5   6

For example the ring N=22 to N=37 is all the points

    2.5 < hypot(X,Y) < 3.5
    where hypot(X,Y) = sqrt(X^2+Y^2)

=head2 N Start

The default is to number points starting N=1 as shown above.  An optional
C<n_start> can give a different start with the same shape.  For example to
start at 0,

=cut

# math-image --path=FilledRings,n_start=0 --all --output=numbers_dash --size=37x16

=pod

          26-25-24          n_start => 0
         /        \
       27 13-12-11 23
      /  /        \  \
    28 14  4--3--2 10 22
     |  |  |     |  |  |
    29 15  5  0--1  9 21
     |  |  |      /  /  /
    30 16  6--7--8 20 36
      \  \        /  /
       31 17-18-19 35
         \        /
        8 32-33-34

The only effect is to push the N values by a constant amount but can help
match N on the axes to counts of X,Y points E<lt> R or similar.

=head1 FUNCTIONS

See L<Math::PlanePath/FUNCTIONS> for the behaviour common to all path
classes.

=over 4

=item C<$path = Math::PlanePath::FilledRings-E<gt>new ()>

=item C<$path = Math::PlanePath::FilledRings-E<gt>new (n_start =E<gt> $n)>

Create and return a new path object.

=back

=head1 OEIS

Entries in Sloane's Online Encyclopedia of Integer Sequences related to this
path include,

=over

L<http://oeis.org/A036704> (etc)

=back

    A036705  first diffs of N on X axis,
               being count of X,Y points n-1/2 < norm <= n+1/2
    A036706  1/4 of those diffs

    n_start=1 (the default)
      A036707  N/2+X-1 along X axis,
                 being count norm <= n+1/2 in half plane
      A036708  (N(X,0)-N(X-1,0))/2+1,
                 first diffs of the half plane count

    n_start=0
      A036704  N on X axis, from X=1 onwards
                 count of X,Y points norm <= n+1/2

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::PixelRings>,
L<Math::PlanePath::Hypot>,
L<Math::PlanePath::MultipleRings>

=head1 HOME PAGE

L<http://user42.tuxfamily.org/math-planepath/index.html>

=head1 LICENSE

Copyright 2012, 2013 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/>.

=cut