File Coverage

blib/lib/Math/PlanePath/HIndexing.pm

Criterion	Covered	Total	%
statement	141	168	83.9
branch	38	54	70.3
condition	6	14	42.8
subroutine	15	20	75.0
pod	5	5	100.0
total	205	261	78.5

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


# http://theinf1.informatik.uni-jena.de/~niedermr/publications.html
#
# Rolf Niedermeier
# http://fpt.akt.tu-berlin.de/niedermr/publications.html
#
#
# H second part down per paper
# |
# | *--*  *  *-
# | |  |  |  |
# | *  *--*  *
# | |        |
# | *  *--*  *
# | |  |  |  |
# | O  *  *--*
# |
# +------------
#
# eight similar to AlternatePaper
#
#                |
#    *--*  *--*  *  *-
#    |  |  |  |  |  |
#  --*  *  *  *--*  *--*
#       |  |           |
#       *  *  *--*--*--*
#    |  |  |
# *--*  *  O  *--*--*--*
# |                    |
# *--*--*--*  *  *  *--*
#             |  |  |
# *--*--*--*  *  *  *-
# |           |  |
# *--*  *--*  *  *  *-
#    |  |  |  |  |  |
#          *--*  *--*
#

package Math::PlanePath::HIndexing;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 128;

use Math::PlanePath;
use Math::PlanePath::Base::NSEW;
@ISA = ('Math::PlanePath::Base::NSEW',
        'Math::PlanePath');

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'round_down_pow',
  'round_up_pow',
  'digit_split_lowtohigh';
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;


use constant n_start => 0;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant diffxy_maximum => 0; # upper octant X<=Y so X-Y<=0
use constant 1.02 _UNDOCUMENTED__dxdy_list_at_n => 9;


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

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

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

  {
    # ENHANCE-ME: get direction without full N+1 calculation
    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 $low = _divrem_mutate ($n, 2);
  ### $low
  ### $n

  my @digits = digit_split_lowtohigh($n,4);
  my $len = ($n*0 + 2) ** scalar(@digits);  # inherit bignum 2

  my $x = 0;
  my $y = 0;
  my $rev = 0;
  my $xinvert = 0;
  my $yinvert = 0;
  while (@digits) {
    my $digit = pop @digits;

    ### $len
    ### $rev
    ### $digit

    my $new_xinvert = $xinvert;
    my $new_yinvert = $yinvert;
    my $xo = 0;
    my $yo = 0;
    if ($rev) {
      if ($digit == 1) {
        $xo = $len-1;
        $yo = $len-1;
        $rev ^= 1;
        $new_yinvert = $yinvert ^ 1;
      } elsif ($digit == 2) {
        $xo = 2*$len-2;
        $yo = 0;
        $rev ^= 1;
        $new_xinvert = $xinvert ^ 1;
      } elsif ($digit == 3) {
        $xo = $len;
        $yo = $len;
      }

    } else {
      if ($digit == 1) {
        $xo = $len-2;
        $yo = $len;
        $rev ^= 1;
        $new_xinvert = $xinvert ^ 1;
      } elsif ($digit == 2) {
        $xo = 1;
        $yo = 2*$len-1;
        $rev ^= 1;
        $new_yinvert = $yinvert ^ 1;
      } elsif ($digit == 3) {
        $xo = $len;
        $yo = $len;
      }
    }

    ### $xo
    ### $yo

    if ($xinvert) {
      $x -= $xo;
    } else {
      $x += $xo;
    }
    if ($yinvert) {
      $y -= $yo;
    } else {
      $y += $yo;
    }

    $xinvert = $new_xinvert;
    $yinvert = $new_yinvert;
    $len /= 2;
  }

  ### final: "$x,$y"

  if ($yinvert) {
    $y -= $low;
  } else {
    $y += $low;
  }

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

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

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

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

  if ($x < 0 || $y < 0 || $x > $y - ($y&1)) {
    return undef;
  }
  if (is_infinite($x)) {
    return $x;
  }
  my ($len, $level) = round_down_pow (int($y/1), 2);
  ### $len
  ### $level
  if (is_infinite($level)) {
    return $level;
  }

  my $n = 0;
  my $npower = $len*$len/2;
  my $rev = 0;
  while (--$level >= 0) {
    ### at: "$x,$y rev=$rev  len=$len n=$n"
    my $digit;
    my $new_rev = $rev;
    if ($y >= $len) {
      $y -= $len;
      if ($x >= $len) {
        ### digit 3 ...
        $digit = 3;
        $x -= $len;
      } else {
        my $yinv = $len-1-$y;
        ### digit 1 or 2: "y reduce to $y,  x cmp ".($yinv-($yinv&1))
        if ($x > $yinv-($yinv&1)) {
          ### digit 2, x invert to: $len-1-$x
          $digit = 2;
          $x = $len-1-$x;
        } else {
          ### digit 1, y invert to: $yinv
          $digit = 1;
          $y = $yinv;
        }
        $new_rev ^= 1;
      }
    } else {
      ### digit 0 ...
      $digit = 0;
    }

    if ($rev) {
      $digit = 3 - $digit;
      ### reversed digit: $digit
    }
    $rev = $new_rev;

    ### add n: $npower*$digit
    $n += $npower*$digit;
    $len /= 2;
    $npower /= 4;
  }

  ### end at: "$x,$y  n=$n rev=$rev"
  ### assert: $x == 0
  ### assert: $y == 0 || $y == 1

  return $n + $y^$rev;
}

# 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);
  ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
  ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
  ### HIndexing rect_to_n_range(): "$x1,$y1 to $x2,$y2"

  # y2 & 1 excluding the X=Y diagonal on odd Y rows
  if ($x2 < 0 || $y2 < 0 || $x1 > $y2 - ($y2&1)) {
    return (1, 0);
  }

  my ($len, $level) = round_down_pow (($y2||1), 2);
  return (0, 2*$len*$len-1);
}


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

sub level_to_n_range {
  my ($self, $level) = @_;
  return (0, 2*4**$level - 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, 2);
  my ($pow,$exp) = round_up_pow ($n+1, 4);
  return $exp;
}

sub _UNDOCUMENTED__level_to_area {
  my ($self, $level) = @_;
  return (2**$level - 1)**2;
}
sub _UNDOCUMENTED__level_to_area_Y {
  my ($self, $level) = @_;
  if ($level == 0) { return 0; }
  return 2**(2*$level-1)  - 2**$level;
}
sub _UNDOCUMENTED__level_to_area_up {
  my ($self, $level) = @_;
  if ($level == 0) { return 0; }
  return 2**(2*$level-1)  - 2**$level + 1;
}

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