File Coverage

blib/lib/Math/PlanePath/ToothpickReplicate.pm

Criterion	Covered	Total	%
statement	110	225	48.8
branch	31	110	28.1
condition	1	15	6.6
subroutine	15	20	75.0
pod	6	6	100.0
total	163	376	43.3

# Copyright 2012, 2013, 2014, 2015 Kevin Ryde

# This file is part of Math-PlanePath-Toothpick.
#
# Math-PlanePath-Toothpick 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-Toothpick 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-Toothpick.  If not, see .

# block_order => 'AB123'
# block_order => 'A1B32' is depth first and finite parts first,
# in parts=1 where single infinite spine
#
# maybe tree methods same structure as ToothpickTree
#

# cf A175262 odd binary length and middle digit 1
#    A175263 odd binary length and middle digit 0
#

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

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


# return ($quotient, $remainder)
sub _divrem {
  my ($n, $d) = @_;
  if (ref $n && $n->isa('Math::BigInt')) {
    my ($quot,$rem) = $n->copy->bdiv($d);
    if (! ref $d || $d < 1_000_000) {
      $rem = $rem->numify;  # plain remainder if fits
    }
    return ($quot, $rem);
  }
  my $rem = $n % $d;
  return (int(($n-$rem)/$d), # exact division stays in UV
          $rem);
}


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

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

use Math::PlanePath::ToothpickTree;
*new = \&Math::PlanePath::ToothpickTree::new;
*x_negative = \&Math::PlanePath::ToothpickTree::x_negative;
*y_negative = \&Math::PlanePath::ToothpickTree::y_negative;
*rect_to_n_range = \&Math::PlanePath::ToothpickTree::rect_to_n_range;
*x_minimum = \&Math::PlanePath::ToothpickTree::x_minimum;
*y_minimum = \&Math::PlanePath::ToothpickTree::y_minimum;
*sumxy_minimum    = \&Math::PlanePath::ToothpickTree::sumxy_minimum;
*sumabsxy_minimum = \&Math::PlanePath::ToothpickTree::sumabsxy_minimum;
*rsquared_minimum = \&Math::PlanePath::ToothpickTree::rsquared_minimum;

use constant parameter_info_array =>
  [ { name      => 'parts',
      share_key => 'parts_toothpickreplicate',
      display   => 'Parts',
      type      => 'enum',
      default   => '4',
      choices   => ['4','3','2','1'],
      choices_display => ['4','3','2','1'],
      description => 'Which parts of the pattern to generate.',
    },
  ];

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

{
  my @x_negative_at_n = (undef,
                         undef,  # 1
                         3,      # 2
                         6,      # 3
                         5,      # 4
                        );
  sub x_negative_at_n {
    my ($self) = @_;
    return $x_negative_at_n[$self->{'parts'}];
  }
}
{
  my @y_negative_at_n = (undef,
                         undef,  # 1
                         undef,  # 2
                         2,      # 3
                         2,      # 4
                        );
  sub y_negative_at_n {
    my ($self) = @_;
    return $y_negative_at_n[$self->{'parts'}];
  }
}

# parts=1 same as parts=4
# parts=2 same as parts=4
# parts=3 same as parts=4
# parts=4    33,-12
#           133,-30
#           333,-112
#          1333,-230
#          3332,-1112    -> 3,-1
use constant dir_maximum_dxdy => (3,-1);

#------------------------------------------------------------------------------
# Fraction covered
# Xlevel = 2^(level+1) - 1
# Ylevel = 2^(level+1)
# Nend = (2*4^(level+1) + 1)/3 - 1
#
# Nend / (Xlevel*Ylevel)
#  -> ((2*4^(level+1) + 1)/3 - 1) / 4^(level+1)
#  -> (2*4^(level+1) + 1)/3 / 4^(level+1)
#  -> 2*4^(level+1)/3 / 4^(level+1)
#  -> 2/3

# Leading diagonal 1,3, 7,11,
#                  23,25,29,43,  +22,22,22,32
#                  87,89,93,97,  +86,86,86,86
#                  109,111,115,171,  +86,128
#                  343
# part2start = (4^level + 5)/3    = 3,7,23,87,343
# sums of part2start(level), but +2 in second half of each
# (3)/3=1
# (3+ 1+5)/3=3
# (3+ 1+5 + 4+5)/3=9


# v              v
# |      ->      |     part 3
# +---h      h---+
#
# +---v      h
# |      ->  |         part 1 rot then part 3
# h          +---v
#
#     v      v
#     |  ->  |         part 3 then part 3 again
# h---+      +---h
#

# v          +---v
# |      ->  |         part 1
# +---h      h
#
#     v      v---+
#     |  ->      |     part 3 then part 1 rot is +90
# h---+          h

# N = (2*4^level + 1)/3 + 1   is first of "level"
# 3N-3 = 2*4^level + 1
# 2*4^level = 3N-4
# 4^(level+1) = 6N-8
#
# part = (2*4^level - 2)/3  many points in "level"
# above = (2*4^(level+1) - 2)/3
#       = (4*2*4^level - 2)/3
#       = 4*(2*4^level - 2/4)/3
#       = 4*(2*4^level - 2)/3 + 4*(+ 2 - 2/4)/3
#       = 4*(2*4^level - 2)/3 + 2
#       = 4*part + 2
# part = (above-2)/4

my @quadrant_to_hdx = (1,-1, -1,1);
my @quadrant_to_vdy = (1, 1, -1,-1);

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

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

  {
    my $int = int($n);
    ### $int
    ### $n
    if ($n != $int) {
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+1);
      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 $parts = $self->{'parts'};
  my $x = 0;
  my $y = 0;
  my $hdx = 1;
  my $hdy = 0;
  my $vdx = 0;
  my $vdy = 1;

  if ($parts eq '2') {
    if ($n == 0) {
      return (0,1);
    }

    # first of a replication level
    # Nlevel = 2*(2*4^level - 2)/3 + 1
    #        = (4*4^level - 4)/3 + 1
    #        = (4*4^level - 4 + 3)/3
    #        = (4*4^level - 1)/3     = 5,21
    # 3N = 4*4^level - 1
    # 4^(level+1) = 3N+1

    my ($len,$level) = round_down_pow(3*$n+1, 4);
    my $three_parts = $len/2;

    ### $len
    ### $level
    ### $three_parts
    ### start this level: ($len-1)/3
    ### n reduced: $n-($len-1)/3

    (my $quadrant, $n) = _divrem ($n-($len-1)/3, $three_parts);
    ### $quadrant
    ### n remainder: $n
    ### assert: $quadrant >= 0
    ### assert: $quadrant <= 1

    $n += ($len/2-2)/3;
    if ($quadrant) { $hdx = -1; }
    ### n in quarter: $n

  } elsif ($parts == 3) {
    if ($n <= 1) {
      return (0,$n);
    }
    # Nend = 3*(2*4^level - 2)/3 + 2
    #      = (2*4^level - 2) + 2
    #      = 2*4^level     = 2,8,32
    # N-1 = 2*4^level
    # 4^(level+1) = 2N-2

    my ($len,$level) = round_down_pow(2*$n, 4);
    my $three_parts = $len/2;

    ### $len
    ### $level
    ### $three_parts
    ### start this level: ($len/2+1)
    ### n reduced: $n-($len/2+1)

    (my $quadrant, $n) = _divrem ($n-$len/2, $three_parts);
    ### $quadrant
    ### n remainder: $n
    ### assert: $quadrant >= 0
    ### assert: $quadrant <= 2

    $n += ($len/2-2)/3;
    ### n in quarter: $n

    if ($quadrant == 0) {
      $hdx = 0;  # rotate -90
      $hdy = -1;
      $vdx = 1;
      $vdy = 0;
      $x = -1; # offset
    } elsif ($quadrant == 2) {
      $hdx = -1;  # mirror
    }

  } elsif ($parts == 4) {
    if ($n <= 2) {
      if ($n == 0) { return (0,0); }
      if ($n == 1) { return (0,1); }
      return (0,-1);  # N==2
    }
    # first of a replication level
    # Nlevel = 4*(2*4^level - 2)/3 + 3
    #        = (8*4^level - 8)/3 + 3
    #        = (8*4^level - 8 + 9)/3
    #        = (8*4^level+1)/3           11,43,171
    # 3N = 8*4^level+1
    # 8*4^level = 3N-1
    # 4^(level+2) = 6N-2
    #
    # first of this level, using level+2
    # Nlevel = (4^(level+2)/2+1)/3
    #        = (4^(level+2)+2)/6
    #
    # three count = 3*(2*4^level - 2)/3 + 2
    #             = 2*4^level
    # 43-11 = 32
    # 172-44 = 128

    # getting level+2 and len = 4^(level+2)
    my ($len,$level) = round_down_pow(6*$n-2, 4);
    my $three_parts = $len/8;

    ### all breakdown ...
    ### $level
    ### $len
    ### $three_parts
    ### Nlevel base: ($len+2)/6

    (my $quadrant, $n) = _divrem ($n-($len+2)/6, $three_parts);
    ### $quadrant
    ### n remainder: $n
    ### assert: $quadrant >= 0
    ### assert: $quadrant <= 3

    # quarter middle
    # Nquarter = (2*4^level - 2)/3  = 2,10,42
    $n += ($len/8-2)/3;
    $hdx = $quadrant_to_hdx[$quadrant];
    $vdy = $quadrant_to_vdy[$quadrant];
    ### n in quarter: $n
  }

  # quarter first of a replication level
  # Nlevel = 4*(2*4^level - 2)/3 + 2
  #        = (8*4^level - 8)/3 + 2
  #        = (8*4^level - 8 + 6)/3
  #        = (8*4^level - 2)/3           2,10,42
  # 3N = 8*4^level-2
  # 8*4^level = 3N+2
  # 4^(level+2) = 6N+4
  #
  # using level+1
  # Nlevel = (8*4^level - 2)/3
  #        = (2*4^(level+1) - 2)/3


  # getting level+2 and 16*len
  my ($len,$level) = round_down_pow(6*$n+4, 4);
  my $part_n = (2*$len-2)/3;
  ### $level
  ### $part_n

  $len = 2**$level;
  for ( ;
        $level-- >= 0;
        $len /= 2,  $part_n = ($part_n-2)/4) {

    ### at: "x=$x,y=$y level=$level hxy=$hdx,$hdy vxy=$vdx,$vdy   n=$n"
    ### $len
    ### $part_n
    ### assert: $len == 2 ** ($level+1)
    ### assert: $part_n == (2 * 4 ** ($level+1) - 2)/3

    if ($n < $part_n) {
      ### part 0, no change ...
      next;
    }

    $n -= $part_n;
    $x += $len * ($hdx + $vdx);  # diagonal
    $y += $len * ($hdy + $vdy);

    if ($n == 0) {
      ### toothpick A ...
      last;
    }
    if ($n == 1) {
      ### toothpick B ...
      $x += $vdx;
      $y += $vdy;
      last;
    }
    $n -= 2;

    if ($n < $part_n) {
      ### part 1, rotate ...
      $x -= $hdx; # offset
      $y -= $hdy;
      ($hdx,$hdy, $vdx,$vdy)    # rotate 90 in direction v toward h
        = (-$vdx,-$vdy, $hdx,$hdy);
      next;
    }
    $n -= $part_n;

    if ($n < $part_n) {
      ### part 2 ...
      next;
    }
    $n -= $part_n;

    ### part 3, mirror ...
    $hdx = -$hdx;
    $hdy = -$hdy;
  }

  ### assert: $n == 0 || $n == 1

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

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

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

  my $parts = $self->{'parts'};
  my $rotated = ($parts == 3 && $x >= 0 && $y < 0);
  if ($rotated) {
    ($x,$y) = (-$y,$x+1);  # rotate +90 and shift up
    ### rotated: "x=$x y=$y"
  }

  my ($len,$level) = round_down_pow (max(abs($x), abs($y)-1),
                                     2);
  if (is_infinite($level)) {
    return $level;
  }
  ### $level
  ### $len

  my $zero = $x * 0 * $y;
  my $n = $zero;

  if ($parts == 2) {
    if ($x == 0) {
      if ($y == 1) { return 0; }
    }
    $n += (2*$len*$len+1)/3;   # +1,+3,+11,+43
    if ($x < 0) {
      $x = -$x;
      $n += 2*$len*$len;  # second quad, +2,+8,+32
    }

  } elsif ($parts == 3) {
    ### 3/4 ...
    if ($x == 0) {
      if ($y == 0) { return 0; }
      if ($y == 1) { return 1; }
    }
    $n += (10*$len*$len+2)/3;   # +4,+14,+54,+214,+854,+3414
    if ($rotated) {
      $n -= 2*$len*$len;  # fourth quad, -2, -8, -32
    } elsif ($x < 0) {
      $x = -$x;
      if ($y > 0) {
        $n += 2*$len*$len;  # second quad, +2, +8, +32
      } else {
        return undef;  # third quad, empty
      }
    }
  } elsif ($parts == 4) {
    if ($x == 0) {
      if ($y == 0)  { return 0; }
      if ($y == 1)  { return 1; }
      if ($y == -1) { return 2; }
    }
    $n += (2*$len*$len+1);
    if ($x < 0) {
      $x = -$x;
      if ($y > 0) {
        $n += 2*$len*$len;  # second quad, +2, +8, +32
      } else {
        $n += 4*$len*$len;  # third quad, +4,+16
        $y = -$y;
      }
    } else {
      if ($y < 0) {
        $n += 6*$len*$len;  # fourth quad
        $y = -$y;
      }
    }
  }

  #                              2^(level+1)-1
  #                              v
  #          +-----------+---------+
  #          |           |         | <- 2^(level+1)
  #          |   3             2   |
  #          | mirror        same  |
  #          |         --B--       | <- 2^level + 1
  #          |           |         |
  #          +--         A       --+ <- 2^level
  #                      |         |
  #                          1     |
  #                         rot    |
  #             0           +90    |
  #                    |           |
  #                    +-----------+
  #                      ^
  #                     2^level

  my $part_n = (2*$len*$len - 2) / 3;
  ### $part_n

  while ($level-- > 0) {
    ### at: "x=$x,y=$y  len=$len part_n=$part_n   n=$n"
    ### assert: $len == 2 ** ($level+1)
    ### assert: $part_n == (2 * 4 ** ($level+1) - 2)/3

    if ($x == $len) {
      if ($y == $len) {
        ### toothpick A ...
        return $n + $part_n;
      }
      if ($y == $len+1) {
        ### toothpick B ...
        return $n + $part_n + 1;
      }
    }

    if ($y <= $len) {
      if ($x < $len) {
        ### part 0 ...
      } else {
        ### part 1, rotate ...
        $n += $part_n + 2;
        ($x,$y) = ($len-$y,$x-$len+1); # shift, rotate +90
      }
    } else {
      $y -= $len;
      if ($x > $len) {
        ### part 2 ...
        $n += 2*$part_n + 2;
        $x -= $len;
      } else {
        ### part 3 ...
        $n += 3*$part_n + 2;
        $x = $len-$x; # mirror
      }
    }

    $len /= 2;
    $part_n = ($part_n-2)/4;
  }

  ### end loop: "x=$x y=$y   n=$n"

  if ($x == 1) {
    if ($y == 1) {
      return $n;
    } elsif ($y == 2) {
      return $n + 1;
    }
  }

  return undef;
}

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

# parts=1
# LevelPoints[k] = 4*LevelPoints[k] + 2  starting LevelPoints[0] = 2
# LevelPoints[k] = 2 + 2*4 + 2*4^2 + ... + 2*4^(k-1) + 4^k*LevelPoints[0]
# LevelPoints[k] = 2 + 2*4 + 2*4^2 + ... + 2*4^(k-1) + 2*4^k
# LevelPoints[k] = 2*(4^(k+1) - 1)/3

my %level_to_n_range = (4 => -2,
                        3 => -3,
                        2 => -4,
                        1 => -5,
                       );
sub level_to_n_range {
  my ($self, $level) = @_;
  return (0,
          (4**($level+1) * (2*$self->{'parts'})
           + $level_to_n_range{$self->{'parts'}}) / 3);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n) + 1;
  _divrem_mutate ($n, 2*$self->{'parts'});
  my ($pow, $exp) = round_down_pow (4*$n - $level_to_n_range{$self->{'parts'}}, 2);
  return $exp + 1;
}

# return $remainder, modify $n
# the scalar $_[0] is modified, but if it's a BigInt then a new BigInt is made
# and stored there, the bigint value is not changed
sub _divrem_mutate {
  my $d = $_[1];
  my $rem;
  if (ref $_[0] && $_[0]->isa('Math::BigInt')) {
    ($_[0], $rem) = $_[0]->copy->bdiv($d);  # quot,rem in array context
    if (! ref $d || $d < 1_000_000) {
      return $rem->numify;  # plain remainder if fits
    }
  } else {
    $rem = $_[0] % $d;
    $_[0] = int(($_[0]-$rem)/$d); # exact division stays in UV
  }
  return $rem;
}

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