File Coverage

blib/lib/Math/PlanePath/HypotOctant.pm

Criterion	Covered	Total	%
statement	98	145	67.5
branch	16	40	40.0
condition	2	14	14.2
subroutine	14	18	77.7
pod	7	7	100.0
total	137	224	61.1

# Copyright 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 .


# Circle drop splash rings from
# math-image --path=HypotOctant --values=DigitProductSteps,values_type=count
# math-image --path=Hypot --values=DigitProduct
# math-image --path=Hypot --values=DigitCount
# math-image --path=Hypot --values=Modulo,modulus=1000
# http://stefan.guninski.com/oeisposter/
#
# pi*r^2 - pi*(r-1)^2 = pi*(2r-1)
# octant is 1/8 of that pi*(2x-1)/8
# pi*(2x-1)/8=100k
# 2x-1 = 100k*8/pi
# x = 100*4/pi*k
#
# 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) )
#
# A057655 num points norm <= n in square lattice.
#
# A036702 num points |z=a+bi| <= n with 0<=a, 0<=b<=a, so octant
# A036703 num points n-1 < z <= n, first diffs?



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

use vars '$VERSION', '@ISA';
$VERSION = 129;
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 X+Y odd.',
    },
  ];

use constant class_x_negative => 0;
use constant class_y_negative => 0;

sub x_minimum {
  my ($self) = @_;
  return ($self->{'points'} eq 'odd'
          ? 1    # odd, line X=Y not included
          : 0);  # octant Y<=X so X-Y>=0
}
# points=odd X=1,Y=0
# otherwise  X=0,Y=0
*sumabsxy_minimum  = \&x_minimum;
*diffxy_minimum    = \&x_minimum;  # X>=Y so X-Y>=0
*absdiffxy_minimum = \&x_minimum;
*rsquared_minimum  = \&x_minimum;

sub absdy_minimum {
  my ($self) = @_;
  return ($self->{'points'} eq 'all'
          ? 0
          : 1);  # never same Y
}

sub dir_minimum_dxdy {
  my ($self) = @_;
  return ($self->{'points'} eq 'all'
          ? (1,0)    # all i=1 to X=1,Y=0
          : (1,1));  # odd,even always at least NE
}
# max direction SE diagonal as anything else is at most tangent to the
# eighth of a circle
use constant dir_maximum_dxdy => (1,-1); # South-East


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

# my @n_to_x = (undef, 0);
# my @n_to_y = (undef, 0);
# my @hypot_to_n = (1);
# my @y_next_x = (1, 1);
# my @y_next_hypot = (1, 2);

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

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

  } elsif ($points eq 'even') {
    $self->{'n_to_x'} = [undef, 0];
    $self->{'n_to_y'} = [undef, 0];
    $self->{'hypot_to_n'} = [1];
    $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->{'opposite_parity'} = 1;

  } elsif ($points eq 'odd') {
    $self->{'n_to_x'} = [undef];
    $self->{'n_to_y'} = [undef];
    $self->{'hypot_to_n'} = [undef];
    $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->{'opposite_parity'} = 0;

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


# at h=x^2+y^2
# step to (x+k)^2+y^2
# is add 2*x*k+k*k

sub _extend {
  my ($self) = @_;
  ### _extend() n: scalar(@{$self->{'n_to_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 = (0);
  my $hypot = $y_next_hypot->[0];
  for (my $i = 1; $i < @$y_next_x; $i++) {
    if ($hypot == $y_next_hypot->[$i]) {
      push @y, $i;
    } elsif ($hypot > $y_next_hypot->[$i]) {
      @y = ($i);
      $hypot = $y_next_hypot->[$i];
    }
  }

  if ($y[-1] == $#$y_next_x) {
    my $y = scalar(@$y_next_x);
    my $x = $y + ($self->{'points'} eq 'odd');
    $y_next_x->[$y] = $x;
    $y_next_hypot->[$y] = $x*$x+$y*$y;
    ### assert: $y_next_hypot->[$y] == $y**2 + $y_next_x->[$y]**2
  }

  ### store: join(' ',map{"$n_to_x->[$_],$n_to_y->[$_]"} 0 .. $#$n_to_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_y, @y;
  push @$n_to_x,
    map {
      my $x = $y_next_x->[$_];
      $y_next_x->[$_] += $self->{'x_inc'};
      $y_next_hypot->[$_]
        += $self->{'x_inc_factor'} * $x + $self->{'x_inc_squared'};
      ### assert: $y_next_hypot->[$_] == $_**2 + $y_next_x->[$_]**2
      $x
    } @y;

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

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

  if ($n < 1) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  {
    my $int = int($n);
    if ($n != $int) {
      my $frac = $n - $int;  # inherit possible BigFloat/BigRat
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+1);
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
  }

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

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

  return ($n_to_x->[$n], $n_to_y->[$n]);
}

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;
  }

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

  if ($x < 0 || $y < 0 || $y > $x) {
    ### outside first octant ...
    return undef;
  }

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

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

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

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

# 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);
  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }

  if ($x2 < 0 || $y2 < 0) {
    return (1, 0);
  }

  # circle area pi*r^2, with r^2 = $x2**2 + $y2**2
  return (1, 1 + int (3.2/8 * (($x2+1)**2 + ($y2+1)**2)));
}

1;
__END__