File Coverage

blib/lib/Math/PlanePath/BetaOmega.pm

Criterion	Covered	Total	%
statement	125	135	92.5
branch	50	62	80.6
condition	4	4	100.0
subroutine	16	16	100.0
pod	3	3	100.0
total	198	220	90.0

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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=BetaOmega --lines --scale=20
#
# math-image --path=BetaOmega --all --output=numbers_dash

# http://www.upb.de/pc2/papers/files/pdfps399main.toappear.ps   # gone
# http://www.uni-paderborn.de/pc2/papers/files/pdfps399main.toappear.ps
# http://wwwcs.upb.de/pc2/papers/files/399.ps   # gone
#
# copy ?
# http://www.cs.uleth.ca/~wismath/cccg/papers/27l.ps


package Math::PlanePath::BetaOmega;
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',
  'bit_split_lowtohigh',
  'digit_split_lowtohigh',
  'digit_join_lowtohigh';

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




use constant n_start => 0;
use constant class_x_negative => 0;
use constant y_negative_at_n => 4;
*xy_is_visited = \&Math::PlanePath::Base::Generic::_xy_is_visited_x_positive;
use constant 1.02 _UNDOCUMENTED__dxdy_list_at_n => 4;


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

# tables generated by tools/beta-omega-table.pl
#
my @next_state = (28, 8,36,88,  8,28,32,76,  4,16,44,64, 16, 4,40,84,
                  12,24,52,72, 24,12,48,92, 20, 0,60,80,  0,20,56,68,
                  68, 4,40,60, 64, 0,60,40, 76,12,48,36, 72, 8,36,48,
                  84,20,56,44, 80,16,44,56, 92,28,32,52, 88,24,52,32,
                  28, 8,36,48,  8,28,32,52,  4,16,44,56, 16, 4,40,60,
                  12,24,52,32, 24,12,48,36, 20, 0,60,40,  0,20,56,44);
my @digit_to_x = (0,0,1,1, 0,1,1,0, 1,0,0,1, 1,1,0,0,
                  1,1,0,0, 1,0,0,1, 0,1,1,0, 0,0,1,1,
                  1,1,0,0, 0,1,1,0, 1,0,0,1, 0,0,1,1,
                  0,0,1,1, 1,0,0,1, 0,1,1,0, 1,1,0,0,
                  0,0,1,1, 0,1,1,0, 1,0,0,1, 1,1,0,0,
                  1,1,0,0, 1,0,0,1, 0,1,1,0, 0,0,1,1);
my @digit_to_y = (0,1,1,0, 0,0,1,1, 0,0,1,1, 0,1,1,0,
                  1,0,0,1, 1,1,0,0, 1,1,0,0, 1,0,0,1,
                  0,1,1,0, 1,1,0,0, 1,1,0,0, 0,1,1,0,
                  1,0,0,1, 0,0,1,1, 0,0,1,1, 1,0,0,1,
                  0,1,1,0, 0,0,1,1, 0,0,1,1, 0,1,1,0,
                  1,0,0,1, 1,1,0,0, 1,1,0,0, 1,0,0,1);
my @xy_to_digit = (0,1,3,2, 0,3,1,2, 1,2,0,3, 3,2,0,1,
                   2,3,1,0, 2,1,3,0, 3,0,2,1, 1,0,2,3,
                   3,2,0,1, 3,0,2,1, 2,1,3,0, 0,1,3,2,
                   1,0,2,3, 1,2,0,3, 0,3,1,2, 2,3,1,0,
                   0,1,3,2, 0,3,1,2, 1,2,0,3, 3,2,0,1,
                   2,3,1,0, 2,1,3,0, 3,0,2,1, 1,0,2,3);
my @min_digit = (0,0,3,0, 0,2,1,1, 2,undef,undef,undef,
                 0,0,1,0, 0,1,3,2, 2,undef,undef,undef,
                 1,0,0,1, 0,0,2,2, 3,undef,undef,undef,
                 3,0,0,2, 0,0,2,1, 1,undef,undef,undef,
                 2,1,1,2, 0,0,3,0, 0,undef,undef,undef,
                 2,2,3,1, 0,0,1,0, 0,undef,undef,undef,
                 3,2,2,0, 0,1,0,0, 1,undef,undef,undef,
                 1,1,2,0, 0,2,0,0, 3,undef,undef,undef,
                 3,0,0,2, 0,0,2,1, 1,undef,undef,undef,
                 3,2,2,0, 0,1,0,0, 1,undef,undef,undef,
                 2,2,3,1, 0,0,1,0, 0,undef,undef,undef,
                 0,0,3,0, 0,2,1,1, 2,undef,undef,undef,
                 1,1,2,0, 0,2,0,0, 3,undef,undef,undef,
                 1,0,0,1, 0,0,2,2, 3,undef,undef,undef,
                 0,0,1,0, 0,1,3,2, 2,undef,undef,undef,
                 2,1,1,2, 0,0,3,0, 0,undef,undef,undef,
                 0,0,3,0, 0,2,1,1, 2,undef,undef,undef,
                 0,0,1,0, 0,1,3,2, 2,undef,undef,undef,
                 1,0,0,1, 0,0,2,2, 3,undef,undef,undef,
                 3,0,0,2, 0,0,2,1, 1,undef,undef,undef,
                 2,1,1,2, 0,0,3,0, 0,undef,undef,undef,
                 2,2,3,1, 0,0,1,0, 0,undef,undef,undef,
                 3,2,2,0, 0,1,0,0, 1,undef,undef,undef,
                 1,1,2,0, 0,2,0,0, 3);
my @max_digit = (0,3,3,1, 3,3,1,2, 2,undef,undef,undef,
                 0,1,1,3, 3,2,3,3, 2,undef,undef,undef,
                 1,1,0,2, 3,3,2,3, 3,undef,undef,undef,
                 3,3,0,3, 3,1,2,2, 1,undef,undef,undef,
                 2,2,1,3, 3,1,3,3, 0,undef,undef,undef,
                 2,3,3,2, 3,3,1,1, 0,undef,undef,undef,
                 3,3,2,3, 3,2,0,1, 1,undef,undef,undef,
                 1,2,2,1, 3,3,0,3, 3,undef,undef,undef,
                 3,3,0,3, 3,1,2,2, 1,undef,undef,undef,
                 3,3,2,3, 3,2,0,1, 1,undef,undef,undef,
                 2,3,3,2, 3,3,1,1, 0,undef,undef,undef,
                 0,3,3,1, 3,3,1,2, 2,undef,undef,undef,
                 1,2,2,1, 3,3,0,3, 3,undef,undef,undef,
                 1,1,0,2, 3,3,2,3, 3,undef,undef,undef,
                 0,1,1,3, 3,2,3,3, 2,undef,undef,undef,
                 2,2,1,3, 3,1,3,3, 0,undef,undef,undef,
                 0,3,3,1, 3,3,1,2, 2,undef,undef,undef,
                 0,1,1,3, 3,2,3,3, 2,undef,undef,undef,
                 1,1,0,2, 3,3,2,3, 3,undef,undef,undef,
                 3,3,0,3, 3,1,2,2, 1,undef,undef,undef,
                 2,2,1,3, 3,1,3,3, 0,undef,undef,undef,
                 2,3,3,2, 3,3,1,1, 0,undef,undef,undef,
                 3,3,2,3, 3,2,0,1, 1,undef,undef,undef,
                 1,2,2,1, 3,3,0,3, 3);

sub n_to_xy {
  my ($self, $n) = @_;
  ### BetaOmega n_to_xy(): $n
  ### hex: sprintf "%#X", $n

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

  my $int = int($n);
  $n -= $int;  # remaining fraction, preserve possible BigFloat/BigRat

  my $zero = $int * 0;  # inherit bignum
  my @ndigits = digit_split_lowtohigh($int,4);
  ### ndigits: join(', ',@ndigits)."   count ".scalar(@ndigits)

  my $state = ($#ndigits & 1 ? 28 : 0);
  my $dirstate   = ($#ndigits & 1 ? 0 : 28); # default if all $ndigit==3
  my @xbits;
  my @ybits;

  foreach my $i (reverse 0 .. $#ndigits) {
    my $ndigit = $ndigits[$i];    # high to low
    $state += $ndigit;
    if ($ndigit != 3) {
      $dirstate = $state;  # lowest non-3 digit
    }

    ### $ndigit
    ### $state
    ### $dirstate
    ### digit_to_x: $digit_to_x[$state]
    ### digit_to_y: $digit_to_y[$state]
    ### next_state: $next_state[$state]

    $xbits[$i] = $digit_to_x[$state];
    $ybits[$i] = $digit_to_y[$state];
    $state = $next_state[$state];
  }

  ### $dirstate
  ### frac: $n
  ### Ymin: - (((4+$zero)**int($#ndigits/2) - 1) * 2 / 3)

  # with $n fractional part
  return ($n * ($digit_to_x[$dirstate+1] - $digit_to_x[$dirstate])
          + digit_join_lowtohigh(\@xbits, 2, $zero),

          $n * ($digit_to_y[$dirstate+1] - $digit_to_y[$dirstate])
          + (digit_join_lowtohigh(\@ybits, 2, $zero)

             # Ymin = - (4^floor(level/2) - 1) * 2 / 3
             - (((4+$zero)**int(scalar(@ndigits)/2) - 1) * 2 / 3)));
}


# ($len,$level) rounded down for $y ...
sub _y_round_down_len_level {
  my ($y) = @_;
  my $pos;
  if ($pos = ($y >= 0)) {
    # eg. 1 becomes 3, or 5 becomes 15, 2^k-1
    $y = 3 * $y;
  } else {
    # eg. -2 becomes 7, or -10 becomes 31, 2^k-1
    $y = 1 - 3*$y;
  }
  my ($len, $level) = round_down_pow($y,2);

  # Make positive y give even level, and negative y give odd level.
  # If positive and odd then reduce, or if negative and even then reduce.
  if (($level & 1) == $pos) {
    $level--;
    $len /= 2;
  }

  return ($len, $level);
}

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

  $x = round_nearest ($x);
  if ($x < 0) {
    return undef;
  }
  if (is_infinite($x)) {
    return $x;
  }
  my @xbits = bit_split_lowtohigh($x);

  $y = round_nearest ($y);
  my $zero = ($x * 0 * $y);
  my ($len, $level) = _y_round_down_len_level ($y);
  ### y: "len=$len  level=$level"

  if ($#xbits > $level) {
    ### increase level to xbits ...
    $level = $#xbits;
    $len = (2+$zero) ** $level;
  }
  ### $len
  ### $level

  $y += (($level&1 ? 4 : 2) * $len - 2) / 3;
  ### offset y to: $y
  if (is_infinite($y)) {
    return $y;
  }
  my @ybits = bit_split_lowtohigh($y);
  my $state = ($level & 1 ? 28 : 0);

  my @ndigits;
  foreach my $i (reverse 0 .. $level) {   # high to low
    ### at: "i=$i state=$state  xbit=".($xbits[$i]||0)." ybit=".($ybits[$i]||0)

    my $ndigit = $xy_to_digit[$state + 2*($xbits[$i]||0) + ($ybits[$i]||0)];
    $ndigits[$i] = $ndigit;
    $state = $next_state[$state+$ndigit];
  }

  return digit_join_lowtohigh(\@ndigits, 4, $zero);
}

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

  $x1 = round_nearest ($x1);
  $x2 = round_nearest ($x2);
  ($x1,$x2) = ($x2,$x1) if $x1 > $x2;

  if ($x2 < 0) {
    return (1, 0);
  }

  $y1 = round_nearest ($y1);
  $y2 = round_nearest ($y2);
  ($y1,$y2) = ($y2,$y1) if $y1 > $y2;

  my ($len, $level) = round_down_pow ($x2, 2);
  ### x len/level: "$len  $level"

  # If y1/y2 both positive or both negative then only look at the bigger of
  # the two.  If y1 negative and y2 positive then consider both.
  foreach my $y (($y2 > 0 ? ($y2) : ()),
                 ($y1 < 0 ? ($y1) : ())) {
    my ($ylen, $ylevel) = _y_round_down_len_level ($y);
    ### y len/level: "$ylen  $ylevel"
    if ($ylevel > $level) {
      $level = $ylevel;
      $len = $ylen;
    }
  }
  if (is_infinite($len)) {
    return (0, $len);
  }

  my $n_min = my $n_max = 0;
  my $y_min = my $y_max = - (4**int(($level+1)/2) - 1) * 2 / 3;
  my $x_min = my $x_max = 0;
  my $min_state = my $max_state = ($level & 1 ? 28 : 0);
  ### $x_min
  ### $y_min

  while ($level >= 0) {
    ### $level
    ### $len
    {
      my $x_cmp = $x_min + $len;
      my $y_cmp = $y_min + $len;
      my $digit = $min_digit[3*$min_state
                             + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0)
                             + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0)];

      # my $xr = ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0);
      # my $yr = ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0);
      # my $key = 3*$min_state + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0) + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0);
      # ### min at: "min_state=$min_state  $x_min,$y_min   cmp $x_cmp,$y_cmp"
      # ### min_state: state_string($min_state)
      # ### $xr
      # ### $yr
      # ### $key
      # ### min digit: $digit
      # ### min key: $key
      # ### y offset: $digit_to_y[$max_state+$digit]

      $n_min = 4*$n_min + $digit;
      $min_state += $digit;
      if ($digit_to_x[$min_state]) { $x_min += $len; }
      $y_min += $len * $digit_to_y[$min_state];
      $min_state = $next_state[$min_state];
    }
    {
      my $x_cmp = $x_max + $len;
      my $y_cmp = $y_max + $len;
      my $digit = $max_digit[3*$max_state
                             + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0)
                             + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0)];

      # my $xr = ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0);
      # my $yr = ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0);
      # my $key = 3*$min_state + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0) + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0);
      # ### max at: "max_state=$max_state  $x_max,$y_max   cmp $x_cmp,$y_cmp"
      # ### $x_cmp
      # ### $y_cmp
      # ### $xr
      # ### $yr
      # ### $key
      # ### max digit: $digit
      # ### x offset: $digit_to_x[$max_state+$digit]
      # ### y offset: $digit_to_y[$max_state+$digit]
      # ### y digit offset: $digit_to_y[$max_state+$digit]
      # ### y min shift part: - ($level&1)

      $n_max = 4*$n_max + $digit;
      $max_state += $digit;
      if ($digit_to_x[$max_state]) { $x_max += $len; }
      $y_max += $len * $digit_to_y[$max_state];
      $max_state = $next_state[$max_state];
    }

    $len = int($len/2);
    $level--;
  }

  return ($n_min, $n_max);
}

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

use Math::PlanePath::HilbertCurve;
*level_to_n_range = \&Math::PlanePath::HilbertCurve::level_to_n_range;
*n_to_level       = \&Math::PlanePath::HilbertCurve::n_to_level;

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