File Coverage

blib/lib/Math/PlanePath/QuintetCentres.pm

Criterion	Covered	Total	%
statement	223	247	90.2
branch	63	74	85.1
condition	15	20	75.0
subroutine	25	31	80.6
pod	9	9	100.0
total	335	381	87.9

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


# Boundary of unit squares:
# 4*3^n
# QuintetCentres unit squares boundary  a(n) = 4*3^(n-3)
# 12,36,108,324,972
# match 12,36,108,324,972
# A003946 G.f.: (1+x)/(1-3*x).
# A025579 a(1)=1, a(2)=2, a(n) = 4*3^(n-3) for n >= 3.
# A027327 a(n) = Sum{(k+1)*T(n,m-k)}, 0<=k<=m, where m=0 for n=0,1; m=n for n >= 2; T given by A026120.


package Math::PlanePath::QuintetCentres;
use 5.004;
use strict;
use POSIX 'ceil';
#use List::Util 'min','max';
*min = \&Math::PlanePath::_min;
*max = \&Math::PlanePath::_max;

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

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

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'digit_split_lowtohigh',
  'round_up_pow';

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


use constant n_start => 0;
use constant parameter_info_array => [ { name      => 'arms',
                                         share_key => 'arms_4',
                                         display   => 'Arms',
                                         type      => 'integer',
                                         minimum   => 1,
                                         maximum   => 4,
                                         default   => 1,
                                         width     => 1,
                                         description => 'Arms',
                                       } ];

{
  my @x_negative_at_n = (undef, 112, 9, 2, 2);
  sub x_negative_at_n {
    my ($self) = @_;
    return $x_negative_at_n[$self->{'arms'}];
  }
}
{
  my @y_negative_at_n = (undef, 2, 4, 6, 7);
  sub y_negative_at_n {
    my ($self) = @_;
    return $y_negative_at_n[$self->{'arms'}];
  }
}

use constant dx_minimum => -1;
use constant dx_maximum => 1;
use constant dy_minimum => -1;
use constant dy_maximum => 1;

*_UNDOCUMENTED__dxdy_list = \&Math::PlanePath::_UNDOCUMENTED__dxdy_list_eight;
{
  my @_UNDOCUMENTED__dxdy_list_at_n = (undef, 18, 14, 11, 11);
  sub _UNDOCUMENTED__dxdy_list_at_n {
    my ($self) = @_;
    return $_UNDOCUMENTED__dxdy_list_at_n[$self->{'arms'}];
  }
}

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

# N=9 first straight, then for other arms 18,27,36
sub _UNDOCUMENTED__turn_any_straight_at_n {
  my ($self) = @_;
  return 9*$self->arms_count;
}


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

sub new {
  my $self = shift->SUPER::new(@_);
  $self->{'arms'} = max(1, min(4, $self->{'arms'} || 1));
  return $self;
}

my @rot_to_x = (0,0,-1,-1);
my @rot_to_y = (0,1,1,0);
my @dir4_to_dx = (1,0,-1,0);
my @dir4_to_dy = (0,1,0,-1);
my @digit_reverse = (0,1,0,0,1);

sub n_to_xy {
  my ($self, $n) = @_;
  ### QuintetCentres n_to_xy(): "arms=$self->{'arms'}   $n"

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

  {
    my $int = int($n);
    if ($n != $int) {
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+$self->{'arms'});
      my $frac = $n - $int;  # inherit possible BigFloat
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int; # BigFloat int() gives BigInt, use that
  }
  my $zero = ($n * 0);   # inherit BigInt 0

  # arm as initial rotation
  my $rot = _divrem_mutate ($n, $self->{'arms'});

  my @digits = digit_split_lowtohigh($n,5);
  my @sx;
  my @sy;
  {
    my $sx = $zero + $dir4_to_dx[$rot];
    my $sy = $zero + $dir4_to_dy[$rot];
    foreach (@digits) {
      push @sx, $sx;
      push @sy, $sy;

      # 2*(sx,sy) + rot+90(sx,sy)
      ($sx,$sy) = (2*$sx - $sy,
                   2*$sy + $sx);
    }
    ### @digits
    my $rev = 0;
    for (my $i = $#digits; $i >= 0; $i--) {  # high to low
      ### digit: $digits[$i]
      if ($rev) {
        ### reverse: "$digits[$i] to ".(5 - $digits[$i])
        $digits[$i] = 4 - $digits[$i];
      }
      $rev ^= $digit_reverse[$digits[$i]];
      ### now rev: $rev
    }
  }
  ### reversed n: @digits


  my $x =
    my $y =
      my $ox =
        my $oy = $zero;

  while (defined (my $digit = shift @digits)) {  # low to high
    my $sx = shift @sx;
    my $sy = shift @sy;
    ### at: "$x,$y  digit $digit   side $sx,$sy"

    # if ($rot & 2) {
    #   ($sx,$sy) = (-$sx,-$sy);
    # }
    # if ($rot & 1) {
    #   ($sx,$sy) = (-$sy,$sx);
    # }

    if ($digit == 0) {
      $x -= $sx;   # left at 180
      $y -= $sy;

    } elsif ($digit == 1) {
      # centre
      ($x,$y) = (-$y,$x);      # rotate -90
      ### rotate to: "$x,$y"
      # $rot--;

    } elsif ($digit == 2) {
      $x += $sy;   # down at -90
      $y -= $sx;
      ### offset to: "$x,$y"

    } elsif ($digit == 3) {
      ($x,$y) = (-$y,$x);      # rotate -90
      $x += $sx;   # right at 0
      $y += $sy;
      # $rot++;

    } else {  # $digit == 4
      ($x,$y) = ($y,-$x);      # rotate +90
      $x -= $sy;   # up at +90
      $y += $sx;
      # $rot++;
    }

    $ox += $sx;
    $oy += $sy;
  }

  ### final: "$x,$y  origin $ox,$oy"
  return ($x + $ox + $rot_to_x[$rot],
          $y + $oy + $rot_to_y[$rot]);
}


# modulus 2*X+Y
#              3
#          0   2   4
#              1
#
#   0 is X=0,Y=0
#
my @modulus_to_x = (0,1,1,1,2);
my @modulus_to_y = (0,-1,0,1,0);

my @modulus_to_digit
  = (0,2,1,4,3,    0,0,10,30,20,     #  0  base
     0,4,3,1,2,    0,10,50,40,10,    # 10
     4,0,1,3,2,    60,20,40,50,20,   # 20  rotated +90
     2,1,3,4,0,    30,60,0,30,50,    # 30
     1,0,3,2,4,    30,20,70,40,40,   # 40
     3,4,1,2,0,    70,10,30,50,50,   # 50  rotated +180
     4,2,3,0,1,    60,60,20,70,10,   # 60
     2,3,1,0,4,    70,0,60,70,40,    # 70  rotated +270
    );
sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### QuintetCentres xy_to_n(): "$x, $y"

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

  foreach my $overflow (2*$x + 2*$y, 2*$x - 2*$y) {
    if (is_infinite($overflow)) { return $overflow; }
  }

  # my $level_limit = log($x*$x + $y*$y + 1) * 1 * 2;
  # if (is_infinite($level_limit)) { return $level_limit; }

  my @digits;
  my $arm;
  my $state;
  for (;;) {
    # if ($level_limit-- < 0) {
    #   ### oops, level limit ...
    #   return undef;
    # }
    if ($x == 0) {
      if ($y == 0) {
        ### found first arm 0,0 ...
        $arm = 0;
        $state = 0;
        last;
      }
      if ($y == 1) {
        ### found second arm 0,1 ...
        $arm = 1;
        $state = 20;
        last;
      }
    } elsif ($x == -1) {
      if ($y == 1) {
        ### found third arm -1,1 ...
        $arm = 2;
        $state = 50;
        last;
      }
      if ($y == 0) {
        ### found fourth arm -1,0 ...
        $arm = 3;
        $state = 70;
        last;
      }
    }
    my $m = (2*$x + $y) % 5;
    ### at: "$x,$y   digits=".join(',',@digits)
    ### mod remainder: $m

    $x -= $modulus_to_x[$m];
    $y -= $modulus_to_y[$m];
    push @digits, $m;

    ### digit: "$m  to $x,$y"
    ### shrink to: ((2*$x + $y) / 5).','.((2*$y - $x) / 5)
    ### assert: (2*$x + $y) % 5 == 0
    ### assert: (2*$y - $x) % 5 == 0

    # shrink
    # (2 -1)  inverse (2  1)
    # (1 2)           (-1 2)
    #
    ($x,$y) = ((2*$x + $y) / 5,
               (2*$y - $x) / 5);
  }

  ### @digits
  my $arms = $self->{'arms'};
  if ($arm >= $arms) {
    return undef;
  }

  my $n = 0;
  foreach my $m (reverse @digits) {  # high to low
    ### $m
    ### digit: $modulus_to_digit[$state + $m]
    ### state: $state
    ### next state: $modulus_to_digit[$state+5 + $m]

    $n = 5*$n + $modulus_to_digit[$state + $m];
    $state = $modulus_to_digit[$state+5 + $m];
  }
  ### final n along arm: $n

  return $n*$arms + $arm;
}

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

# whole plane covered when arms==4
sub xy_is_visited {
  my ($self, $x, $y) = @_;
  return ($self->{'arms'} == 4
          || defined($self->xy_to_n($x,$y)));
}

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

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

  my ($r_lo, $r_hi) = _rect_to_radius_range($x1,$y1, $x2,$y2);
  $r_hi *= 2;
  my $level_plus_1 = ceil( log(max(1,$r_hi/4)) / log(sqrt(5)) ) + 2;

  # Simple over-estimate would be: return (0, 5**$level_plus_1);

  my $level_limit = $level_plus_1;
  ### $level_limit
  if (is_infinite($level_limit)) { return ($level_limit,$level_limit); }

  $x1 = round_nearest ($x1);
  $y1 = round_nearest ($y1);
  $x2 = round_nearest ($x2);
  $y2 = round_nearest ($y2);
  ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
  ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
  ### sorted range: "$x1,$y1  $x2,$y2"

  my $rect_dist = sub {
    my ($x,$y) = @_;
    my $xd = ($x < $x1 ? $x1 - $x
              : $x > $x2 ? $x - $x2
              : 0);
    my $yd = ($y < $y1 ? $y1 - $y
              : $y > $y2 ? $y - $y2
              : 0);
    return ($xd*$xd + $yd*$yd);
  };

  my $arms = $self->{'arms'};
  ### $arms
  my $n_lo;
  {
    my @hypot = (4);
    my $top = 0;
    for (;;) {
    ARM_LO: foreach my $arm (0 .. $arms-1) {
        my $i = 0;
        my @digits;
        if ($top > 0) {
          @digits = ((0)x($top-1), 1);
        } else {
          @digits = (0);
        }

        for (;;) {
          my $n = 0;
          foreach my $digit (reverse @digits) { # high to low
            $n = 5*$n + $digit;
          }
          $n = $n*$arms + $arm;
          ### lo consider: "i=$i  digits=".join(',',reverse @digits)."  is n=$n"

          my ($nx,$ny) = $self->n_to_xy($n);
          my $nh = &$rect_dist ($nx,$ny);
          if ($i == 0 && $nh == 0) {
            ### lo found inside: $n
            if (! defined $n_lo || $n < $n_lo) {
              $n_lo = $n;
            }
            next ARM_LO;
          }

          if ($i == 0 || $nh > $hypot[$i]) {
            ### too far away: "nxy=$nx,$ny   nh=$nh vs ".$hypot[$i]

            while (++$digits[$i] > 4) {
              $digits[$i] = 0;
              if (++$i <= $top) {
                ### backtrack up ...
              } else {
                ### not found within this top and arm, next arm ...
                next ARM_LO;
              }
            }
          } else {
            ### lo descend ...
            ### assert: $i > 0
            $i--;
            $digits[$i] = 0;
          }
        }
      }

      # if an $n_lo was found on any arm within this $top then done
      if (defined $n_lo) {
        last;
      }

      ### lo extend top ...
      if (++$top > $level_limit) {
        ### nothing below level limit ...
        return (1,0);
      }
      $hypot[$top] = 5 * $hypot[$top-1];
    }
  }

  my $n_hi = 0;
 ARM_HI: foreach my $arm (reverse 0 .. $arms-1) {
    my @digits = ((4) x $level_limit);
    my $i = $#digits;
    for (;;) {
      my $n = 0;
      foreach my $digit (reverse @digits) { # high to low
        $n = 5*$n + $digit;
      }
      $n = $n*$arms + $arm;
      ### hi consider: "arm=$arm  i=$i  digits=".join(',',reverse @digits)."  is n=$n"

      my ($nx,$ny) = $self->n_to_xy($n);
      my $nh = &$rect_dist ($nx,$ny);
      if ($i == 0 && $nh == 0) {
        ### hi found inside: $n
        if ($n > $n_hi) {
          $n_hi = $n;
          next ARM_HI;
        }
      }

      if ($i == 0 || $nh > (4 * 5**$i)) {
        ### too far away: "$nx,$ny   nh=$nh vs ".(4 * 5**$i)

        while (--$digits[$i] < 0) {
          $digits[$i] = 4;
          if (++$i < $level_limit) {
            ### hi backtrack up ...
          } else {
            ### hi nothing within level limit for this arm ...
            next ARM_HI;
          }
        }

      } else {
        ### hi descend
        ### assert: $i > 0
        $i--;
        $digits[$i] = 4;
      }
    }
  }

  if ($n_hi == 0) {
    ### oops, lo found but hi not found
    $n_hi = $n_lo;
  }

  return ($n_lo, $n_hi);
}

#------------------------------------------------------------------------------
# levels

# level=0
# level=1  0 to 4
# level=2  0 to 24    is 5^level-1
#
# multiple arms the same full points of arms=1
# so arms*5^level points numbered starting 0
#        = 5^level*arms - 1
sub level_to_n_range {
  my ($self, $level) = @_;
  return (0,  5**$level * $self->{'arms'} - 1);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  _divrem_mutate ($n, $self->{'arms'});
  my ($pow, $exp) = round_up_pow ($n+1, 5);
  return $exp;
}

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