File Coverage

blib/lib/Math/PlanePath/QuadricIslands.pm

Criterion	Covered	Total	%
statement	75	121	61.9
branch	6	26	23.0
condition	0	9	0.0
subroutine	22	26	84.6
pod	6	6	100.0
total	109	188	57.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 .


# math-image --path=QuadricIslands --lines --scale=10
# math-image --path=QuadricIslands --all --output=numbers_dash --size=132x50


package Math::PlanePath::QuadricIslands;
use 5.004;
use strict;
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

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

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

use Math::PlanePath::QuadricCurve;

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


use constant n_frac_discontinuity => 0;
use constant x_negative_at_n => 1;
use constant y_negative_at_n => 1;
use constant sumabsxy_minimum => 1;   # minimum X=1/2,Y=1/2
use constant rsquared_minimum => 0.5; # minimum X=1/2,Y=1/2

use constant dx_maximum => 1;
use constant dy_maximum => 1;
use constant dsumxy_maximum => 1;
use constant ddiffxy_minimum => -1; # dDiffXY=+1 or -1
use constant ddiffxy_maximum => 1;
use constant dir_maximum_dxdy => (0,-1); # South

# N=1,2,3,4  gcd(1/2,1/2) = 1/2
use constant gcdxy_minimum => 1/2;

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

# N=1 to 4      4 of, level=0
# N=5 to 36    12 of, level=1
# N=37 to ..   48 of, level=3
#
# each loop = 4*8^level
#
#     n_base = 1 + 4*8^0 + 4*8^1 + ... + 4*8^(level-1)
#            = 1 + 4*[ 8^0 + 8^1 + ... + 8^(level-1) ]
#            = 1 + 4*[ (8^level - 1)/7 ]
#            = 1 + 4*(8^level - 1)/7
#            = (4*8^level - 4 + 7)/7
#            = (4*8^level + 3)/7
#
#     n >= (4*8^level + 3)/7
#     7*n = 4*8^level + 3
#     (7*n - 3)/4 = 8^level
#
#    nbase(k+1)-nbase(k)
#       = (4*8^(k+1)+3  - (4*8^k+3)) / 7
#       = (4*8*8^k - 4*8^k) / 7
#       = (4*8-4) * 8^k / 7
#       = 28 * 8^k / 7
#       = 4 * 8^k
#
#    nbase(0) = (4*8^0 + 3)/7 = (4+3)/7 = 1
#    nbase(1) = (4*8^1 + 3)/7 = (4*8+3)/7 = (32+3)/7 = 35/7 = 5
#    nbase(2) = (4*8^2 + 3)/7 = (4*64+3)/7 = (256+3)/7 = 259/7 = 37
#
### loop 1: 4* 8**1
### loop 2: 4* 8**2
### loop 3: 4* 8**3

# sub _level_to_base {
#   my ($level) = @_;
#   return (4*8**$level + 3) / 7;
# }
# ### level_to_base(1): _level_to_base(1)
# ### level_to_base(2): _level_to_base(2)
# ### level_to_base(3): _level_to_base(3)

# base = (4 * 8**$level + 3)/7
#      = (4 * 8**($level+1) / 8 + 3)/7
#      = (8**($level+1) / 2 + 3)/7
sub _n_to_base_and_level {
  my ($n) = @_;
  my ($base,$level) = round_down_pow ((7*$n - 3)*2, 8);
  return (($base/2 + 3)/7,
          $level - 1);
}

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

  my ($base, $level) = _n_to_base_and_level($n);
  ### $level
  ### $base
  ### level: "$level"
  ### next base would be: (4 * 8**($level+1) + 3)/7

  my $rem = $n - $base;
  ### $rem

  ### assert: $n >= $base
  ### assert: $n < 8**($level+1)
  ### assert: $rem>=0
  ### assert: $rem < 4 * 8 ** $level

  my $sidelen = 8**$level;
  my $side = int($rem / $sidelen);
  ### $sidelen
  ### $side
  ### $rem
  $rem -= $side*$sidelen;
  ### assert: $side >= 0 && $side < 4
  my ($x, $y) = Math::PlanePath::QuadricCurve::n_to_xy ($self, $rem);

  my $pos = 4**$level / 2;
  ### side calc: "$x,$y   for pos $pos"
  ### $x
  ### $y

  if ($side < 1) {
    ### horizontal rightwards
    return ($x - $pos,
            $y - $pos);
  } elsif ($side < 2) {
    ### right vertical upwards
    return (-$y + $pos,     # rotate +90, offset
            $x - $pos);
  } elsif ($side < 3) {
    ### horizontal leftwards
    return (-$x + $pos,     # rotate 180, offset
            -$y + $pos);
  } else {
    ### left vertical downwards
    return ($y - $pos,     # rotate -90, offset
            -$x + $pos);
  }
}

#   +-------+-------+-------+
#   |31     | 24 0,1|     23|
#   |       |       |       |
#   |   +-------+-------+   |
#   |   |4  |   |3  |   |   |
#   |   |   |   |   |   |   |
#   +---|--- ---|--- ---|---+   Y=0.5
#   |32 |   |   |   |   | 16|
#   |   |   |   |   |   |   |
#   |   +=======+=======+   |   Y=0
#   |   |1  |   |2  |   |   |
#   |   |   |   |   |   |   |
#   +---|--- ---|--- ---|---+   Y=-0.5
#   |   |   |   |   |   |   |
#   |   |   |   |   |   |   |
#   |   +-------+-------+   |   Y=-1
#   |       |       |       |
#   |7      |8      |     15|
#   +-------+-------+-------+
#
# -2 <= 2*x < 2, round to -2,-1,0,1
# then 4*yround -8,-4,0,4
# total -10 to 5 inclusive

my @inner_n_list = ([1,7], [1,8], [2,8], [2,15],  # Y=-1
                    [1,32], [1],   [2],  [2,16],  # Y=-0.5
                    [4,32], [4],   [3],  [3,16],  # Y=0
                    [4,31],[4,24],[3,24],[3,23]); # Y=0.5

sub xy_to_n {
  return scalar((shift->xy_to_n_list(@_))[0]);
}
sub xy_to_n_list {
  my ($self, $x, $y) = @_;
  ### QuadricIslands xy_to_n(): "$x, $y"

  if ($x >= -1 && $x < 1
      && $y >= -1 && $y < 1) {

    ### round 2x: floor(2*$x)
    ### round 2y: floor(2*$y)
    ### index: floor(2*$x) + 4*floor(2*$y) + 10
    ### assert: floor(2*$x) + 4*floor(2*$y) + 10 >= 0
    ### assert: floor(2*$x) + 4*floor(2*$y) + 10 <= 15

    return @{$inner_n_list[ floor(2*$x)
                            + 4*floor(2*$y)
                            + 10 ]};
  }

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

  my $high;
  if ($x >= $y + ($y>0)) {
    # +($y>0) to exclude the downward bump of the top side
    ### below leading diagonal ...
    if ($x < -$y) {
      ### bottom quarter ...
      $high = 0;
    } else {
      ### right quarter ...
      $high = 1;
      ($x,$y) = ($y, -$x);   # rotate -90
    }
  } else {
    ### above leading diagonal
    if ($y > -$x) {
      ### top quarter ...
      $high = 2;
      $x = -$x;   # rotate 180
      $y = -$y;
    } else {
      ### right quarter ...
      $high = 3;
      ($x,$y) = (-$y, $x);   # rotate +90
    }
  }
  ### rotate to: "$x,$y   high=$high"

  # ymax = (10*4^(l-1)-1)/3
  # ymax < (10*4^(l-1)-1)/3+1
  # (10*4^(l-1)-1)/3+1 > ymax
  # (10*4^(l-1)-1)/3 > ymax-1
  # 10*4^(l-1)-1 > 3*(ymax-1)
  # 10*4^(l-1) > 3*(ymax-1)+1
  # 10*4^(l-1) > 3*(ymax-1)+1
  # 10*4^(l-1) > 3*ymax-3+1
  # 10*4^(l-1) > 3*ymax-2
  # 4^(l-1) > (3*ymax-2)/10
  #
  # (2*4^(l-1) + 1)/3 = ymin
  # 2*4^(l-1) + 1 = 3*y
  # 2*4^(l-1) = 3*y-1
  # 4^(l-1) = (3*y-1)/2
  #
  # ypos = 4^l/2 = 2*4^(l-1)


  # z = -2*y+x
  # (2*4**($level-1) + 1)/3 = z
  # 2*4**($level-1) + 1 = 3*z
  # 2*4**($level-1) = 3*z - 1
  # 4**($level-1) = (3*z - 1)/2
  #               = (3*(-2y+x)-1)/2
  #               = (-6y+3x - 1)/2
  #               = -3*y + (3x-1)/2

  # 2*4**($level-1) = -2*y-x
  # 4**($level-1) = -y-x/2
  # 4**$level = -4y-2x
  #
  # line slope y/x = 1/2 as an index
  my $z = -$y-$x/2;
  my ($len,$level) = round_down_pow ($z, 4);
  ### $z
  ### amin: 2*4**($level-1)
  ### $level
  ### $len
  if (is_infinite($level)) {
    return $level;
  }

  $len *= 2;
  $x += $len;
  $y += $len;
  ### shift to: "$x,$y"
  my $n = Math::PlanePath::QuadricCurve::xy_to_n($self, $x, $y);

  # Nmin = (4*8^l+3)/7
  # Nmin+high = (4*8^l+3)/7 + h*8^l
  #           = (4*8^l + 3 + 7h*8^l)/7 +
  #           = ((4+7h)*8^l + 3)/7
  #
  ### plain curve on: ($x+2*$len).",".($y+2*$len)."  give n=".(defined $n && $n)
  ### $high
  ### high: (8**$level)*$high
  ### base: (4 * 8**($level+1) + 3)/7
  ### base with high: ((4+7*$high) * 8**($level+1) + 3)/7

  if (defined $n) {
    return ((4+7*$high) * 8**($level+1) + 3)/7 + $n;
  } else {
    return;
  }
}

# level width extends
#    side = 4^level
#    ypos = 4^l / 2
#    width = 1 + 4 + ... + 4^(l-1)
#          = (4^l - 1)/3
#    ymin = ypos(l) - 4^(l-1) - width(l-1)
#         = 4^l / 2  - 4^(l-1) - (4^(l-1) - 1)/3
#         = 4^(l-1) * (2 - 1 - 1/3) + 1/3
#         = (2*4^(l-1) + 1) / 3
#
#    (2*4^(l-1) + 1) / 3 = y
#    2*4^(l-1) + 1 = 3*y
#    2*4^(l-1) = 3*y-1
#    4^(l-1) = (3*y-1)/2
#
# ENHANCE-ME: slope Y=X/2+1 or thereabouts for sides
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### QuadricIslands rect_to_n_range(): "$x1,$y1  $x2,$y2"

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

  my $m = max(abs($x1), abs($x2),
              abs($y1), abs($y2));

  my ($len,$level) = round_down_pow (6*$m-2, 4);
  ### $len
  ### $level
  return (1,
          (32*8**$level - 4)/7);
}

#    ymax = ypos(l) + 4^(l-1) + width(l-1)
#         = 4^l / 2  + 4^(l-1) + (4^(l-1) - 1)/3
#         = 4^(l-1) * (4/2 + 1 + 1/3) - 1/3
#         = 4^(l-1) * (2 + 1 + 1/3) - 1/3
#         = 4^(l-1) * 10/3 - 1/3
#         = (10*4^(l-1) - 1) / 3
#
#    (10*4^(l-1) - 1) / 3 = y
#    10*4^(l-1) - 1 = 3*y
#    10*4^(l-1) = 3*y+1
#    4^(l-1) = (3*y+1)/10
#
# based on max ??? ...
#
# my ($power,$level) = round_down_pow ((3*$m+1-3)/10, 4);
# ### $power
# ### $level
# return (1,
#         (4*8**($level+3) + 3)/7 - 1);

#------------------------------------------------------------------------------
# Nstart(k) = (4*8^k + 3)/7
# Nend(k) = Nstart(k+1) - 1
#         = (4*8*8^k + 3)/7 - 1
#         = (4*8*8^k + 8*3 - 8*3 + 3)/7 - 1
#         = (4*8*8^k + 8*3)/7 + (-8*3 + 3)/7 - 1
#         = 8*Nstart(k) + (-8*3 + 3)/7 - 1
#         = 8*Nstart(k) - 4

sub level_to_n_range {
  my ($self, $level) = @_;
  my $n_lo = (4 * 8**$level + 3)/7;
  return ($n_lo, 8*$n_lo - 4);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 1) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  my ($base,$level) = _n_to_base_and_level($n);
  return $level;
}

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