File Coverage

blib/lib/Math/PlanePath/GosperIslands.pm

Criterion	Covered	Total	%
statement	186	217	85.7
branch	30	46	65.2
condition	4	5	80.0
subroutine	27	27	100.0
pod	4	4	100.0
total	251	299	83.9

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



# math-image --path=GosperIslands --lines --scale=10
# math-image --path=GosperIslands --all --output=numbers_dash
#
# fractal dimension  2*log(3) / log(7) = 1.12915, A113211 decimal

# Check:
# A017926   boundary length, unit equilateral triangles
# A229977   boundary length, initial hexagon=1


package Math::PlanePath::GosperIslands;
use 5.004;
use strict;
use POSIX 'ceil';
use Math::Libm 'hypot';
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

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

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

use Math::PlanePath::SacksSpiral;

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


use constant n_frac_discontinuity => 0;
use constant x_negative_at_n => 3;
use constant y_negative_at_n => 5;
use constant sumabsxy_minimum => 2; # minimum X=2,Y=0 or X=1,Y=1
use constant rsquared_minimum => 2; # minimum X=1,Y=1

# jump across rings is upwards, so dY maximum unbounded but minimum=-1
use constant dy_minimum => -1;

# dX and dY unbounded jumping between rings, with the jump position rotating
# around slowly with the twistiness of the ring.
use constant absdx_minimum => 1;

# jump across rings is ENE, so dSum maximum unbounded but minimum=-2
use constant dsumxy_minimum => -2;

# jump across rings is ENE, so dDiffXY minimum unbounded but maximum=+2
use constant ddiffxy_maximum => 2;

use constant dir_maximum_dxdy => (1,-1); # South-East
use constant turn_any_straight => 0; # never straight


#------------------------------------------------------------------------------
sub new {
  my $self = shift->SUPER::new (@_);
  $self->{'sides'} ||= 6; # default
  return $self;
}

# innermost hexagon level 0
#     level 0    len=6
#     level 1    len=18
# each sidelen = 3^level
# each ringlen = 6*3^level
#
# Nstart(level) = 1 + ringlen(0) + ... + ringlen(level-1)
#               = 1 + 6* [ 3^0 + ... + 3^(level-1) ]
#               = 1 + 6* [ (3^level - 1) / 2 ]
#               = 1 + 3*(3^level - 1)
#               = 1 + 3*3^level - 3
#               = 3^(level+1) - 2
#
# 3^(level+1) = N+2
# level+1 = log3(N+2)
# level = log3(N+2) - 1
#
# sides=3
# ringlen = 3*3^level
# Nstart(level) = 1 + 3* [ 3^0 + ... + 3^(level-1) ]
#               = 1 + 3* [ (3^level - 1) / 2 ]
#               = 1 + 3*(3^level - 1)/2
#               = 1 + (3*3^level - 3)/2
#               = (3*3^level - 3 + 2)/2
#               = (3^(level+1) - 1)/2
# 3^(level+1) - 1 = 2*N
# 3^(level+1) = 2*N+1

my @_xend = (2, 5);
my @_yend = (0, 1);
my @_hypotend = (4, 28);
sub _ends_for_level {
  my ($level) = @_;
  if ($#_xend < $level) {
    my $x = $_xend[-1];
    my $y = $_yend[-1];
    do {
      ($x,$y) = ((5*$x - 3*$y)/2,   # 2*$x + rotate +60
                 ($x + 5*$y)/2);    # 2*$y + rotate +60
      ### _ends_for_level() push: scalar(@_xend)."  $x,$y"
      # ### assert: "$x,$y" eq join(','__PACKAGE__->n_to_xy(scalar(@xend) ** 3))
      push @_xend, $x;
      push @_yend, $y;
      push @_hypotend, $_hypotend[-1] * 7;
    } while ($#_xend < $level);
  }
}

my @level_x = (0);
my @level_y = (0);

sub n_to_xy {
  my ($self, $n) = @_;
  ### GosperIslands n_to_xy(): $n
  if ($n < 1) {
    return;
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }

  my $sides = $self->{'sides'};
  my ($pow, $level) = round_down_pow (($sides == 6 ? $n+2 : 2*$n+1),
                                      3);
  my $base = ($sides == 6 ? $pow-2 : ($pow-1)/2);
  ### $level
  ### $base

  $n -= $base;  # remainder
  my $sidelen = $pow / 3;
  $level--;
  my $side = int ($n / $sidelen);
  my ($x, $y) = _side_n_to_xy ($n - $side*$sidelen);

  _ends_for_level($level);
  ### raw xy: "$x,$y"
  ### pos: "$_xend[$level],$_yend[$level]"

  if ($sides == 6) {
    ($x,$y) = (($x+3*$y)/-2,             # rotate +120
               ($x-$y)/2);
    ### rotate xy: "$x,$y"

    $x += $_xend[$level];
    $y += $_yend[$level];
    if ($side >= 3) {
      $x = -$x;   # rotate 180
      $y = -$y;
      $side -= 3;
    }
    if ($side == 1) {
      ($x,$y) = (($x-3*$y)/2,   # rotate +60
                 ($x+$y)/2);
    } elsif ($side == 2) {
      ($x,$y) = (($x+3*$y)/-2,  # rotate +120
                 ($x-$y)/2);
    }
  } else {
    my $xend = $_xend[$level];
    my $yend = $_yend[$level];
    my $xp = 0;
    my $yp = 0;
    foreach (0 .. $level-1) {
      $xp +=$_xend[$_];
      $yp +=$_yend[$_];
    }
    ### ends: "$xend,$yend prev $xp,$yp"
    ($xp,$yp) = (($xp-3*$yp)/-2,   # rotate -120
                 ($xp+$yp)/-2);
    # ($xp,$yp) = (($xp-3*$yp)/2,   # rotate +60
    #              ($xp+$yp)/2);
    if ($side == 0) {
      $x += $xp;
      $y += $yp;
    } elsif ($side == 1) {
      ($x,$y) = (($x+3*$y)/-2,  # rotate +120
                 ($x-$y)/2);
      $x += $xp;
      $y += $yp;
      $x += $xend;
      $y += $yend;
    } elsif ($side == 2) {
      ($x,$y) = (($x-3*$y)/-2,   # rotate +240==-120
                 ($x+$y)/-2);
      $x += $xp;
      $y += $yp;
      $x += $xend;
      $y += $yend;
      ($xend,$yend) = (($xend+3*$yend)/-2,  # rotate +120
                       ($xend-$yend)/2);
      $x += $xend;
      $y += $yend;
    }
  }
  return ($x,$y);
}

# ENHANCE-ME: share with GosperSide ?
sub _side_n_to_xy {
  my ($n) = @_;
  ### _side_n_to_xy(): $n
  if ($n < 0) {
    return;
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }

  my $x;
  my $y = 0;
  {
    my $int = int($n);
    $x = 2*($n - $int);
    $n = $int;  # BigFloat int() gives BigInt, use that
  }
  my $xend = 2;
  my $yend = 0;

  foreach my $digit (digit_split_lowtohigh($n,3)) {
    my $xend_offset = 3*($xend-$yend)/2;   # end and end +60
    my $yend_offset = ($xend+3*$yend)/2;

    ### at: "$x,$y"
    ### $digit
    ### $xend
    ### $yend
    ### $xend_offset
    ### $yend_offset

    if ($digit == 1) {
      ($x,$y) = (($x-3*$y)/2  + $xend,   # rotate +60
                 ($x+$y)/2    + $yend);
    } elsif ($digit == 2) {
      $x += $xend_offset;   # offset and offset +60
      $y += $yend_offset;
    }
    $xend += $xend_offset;   # offset and offset +60
    $yend += $yend_offset;
  }

  ### final: "$x,$y"
  return ($x, $y);
}


#------------------------------------------------------------------------------
# xy_to_n()

# innermost origin 0,0 N=0 level 0,
# first hexagon level 1
#
# Nstart(level) = 3^level - 2
# sidelength(level) = 3^(level-1)
# so Nstart(level) + side * sidelength(level)
#  = 3^level - 2 + side * 3^(level-1)
#  = (3 + side) * 3^(level-1) - 2
#
sub xy_to_n {
  my ($self, $x_full, $y_full) = @_;
  $x_full = round_nearest($x_full);
  $y_full = round_nearest($y_full);
  ### GosperIslands xy_to_n(): "$x_full,$y_full"

  foreach my $overflow (3*$x_full+3*$y_full, 3*$x_full-3*$y_full) {
    if (is_infinite($overflow)) { return $overflow; }
  }
  if (($x_full + $y_full) % 2) {
    ### odd point, not reached ...
    return undef;
  }

  my $r = hypot($x_full,$y_full);
  my $level_limit = ceil(log($r+1)/log(sqrt(7)));
  ### $r
  ### $level_limit
  if (is_infinite($level_limit)) {
    return $level_limit;
  }

  for my $level (0 .. $level_limit + 1) {
    _ends_for_level($level);
    my $xend = $_xend[$level];
    my $yend = $_yend[$level];

    my $x = $x_full - $xend;
    my $y = $y_full - $yend;
    ### level end: "$xend,$yend"
    ### level subtract to: "$x,$y"
    ($x,$y) = ((3*$y-$x)/2,  # rotate -120;
               ($x+$y)/-2);
    ### level rotate to: "$x,$y"

    foreach my $side (0 .. 5) {
      ### try: "level=$level side=$side  $x,$y"
      if (defined (my $n = _xy_to_n_in_level($x,$y,$level))) {
        return ($side+3)*3**$level + $n - 2;
      }
      $x -= $xend;
      $y -= $yend;
      ### subtract to: "$x,$y"
      ($x,$y) = (($x+3*$y)/2,   # rotate -60
                 ($y-$x)/2);

    }
  }
  return undef;
}

sub _xy_to_n_in_level {
  my ($x, $y, $level) = @_;

  _ends_for_level($level);
  my @pending_n = (0);
  my @pending_x = ($x);
  my @pending_y = ($y);
  my @pending_level = ($level);

  while (@pending_n) {
    my $n = pop @pending_n;
    $x = pop @pending_x;
    $y = pop @pending_y;
    $level = pop @pending_level;
    ### consider: "$x,$y  n=$n level=$level"

    if ($level == 0) {
      if ($x == 0 && $y == 0) {
        return $n;
      }
      ### level=0 and not 0,0
      next;
    }
    if (($x*$x+3*$y*$y) > $_hypotend[$level] * 1.1) {
      ### radius out of range: ($x*$x+3*$y*$y)." cf end ".$_hypotend[$level]
      next;
    }

    $level--;
    $n *= 3;
    my $xend = $_xend[$level];
    my $yend = $_yend[$level];
    ### descend: "end=$xend,$yend on level=$level"

    # digit 0
    push @pending_n, $n;
    push @pending_x, $x;
    push @pending_y, $y;
    push @pending_level, $level;
    ### push: "$x,$y  digit=0"

    # digit 1
    $x -= $xend;
    $y -= $yend;
    ($x,$y) = (($x+3*$y)/2,   # rotate -60
               ($y-$x)/2);
    push @pending_n, $n + 1;
    push @pending_x, $x;
    push @pending_y, $y;
    push @pending_level, $level;
    ### push: "$x,$y  digit=1"

    # digit 2
    $x -= $xend;
    $y -= $yend;
    ($x,$y) = (($x-3*$y)/2,   # rotate +60
               ($x+$y)/2);
    push @pending_n, $n + 2;
    push @pending_x, $x;
    push @pending_y, $y;
    push @pending_level, $level;
    ### push: "$x,$y  digit=2"
  }

  return undef;
}

# Each
#           *---
#          /
#      ---*
# is width=5 heightflat=1 is
#     hypot^2 = 5*5 + 3 * 1*1
#             = 25+3
#             = 28
#     hypot = 2*sqrt(7)
#
# comes in closer to
#     level=1   n=9,x=2,y=2  is hypot=sqrt(2*2+3*2*2) = sqrt(16) = 4
#     level=2   n=31,x=2,y=6 is hypot=sqrt(2*2+3*6*6) = sqrt(112) = sqrt(7)*4
# so
#     radius = 4 * sqrt(7)^(level-1)
#     radius/4 = sqrt(7)^(level-1)
#     level-1 = log(radius/4) / log(sqrt(7))
#     level = log(radius/4) / log(sqrt(7)) + 1
#
# Or
#     level=1   n=9,x=2,y=2  is h=2*2+3*2*2 = 16
#     level=2   n=31,x=2,y=6 is h=2*2+3*6*6 = 112 = 7*16
#     h = 16 * 7^(level-1)
#     h/16 = 7^(level-1)
#     level-1 = log(h/16) / log(7)
#     level = log(h/16) / log(7) + 1
#
# is the next outer level, so high covering is the end of the previous
# level,
#
# Nstart(level) - 1 = 3^(level+2) - 2 - 1
#                   = 3^(level+2) - 3
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  $y1 *= sqrt(3);
  $y2 *= sqrt(3);
  my ($r_lo, $r_hi) = Math::PlanePath::SacksSpiral::_rect_to_radius_range
    ($x1,$y1, $x2,$y2);
  $r_hi *= 2;
  my $level_plus_1 = ceil( log(max(1,$r_hi/4)) / log(sqrt(7)) ) + 2;
  return (1, 3**$level_plus_1 - 3);
}

1;
__END__