File Coverage

blib/lib/Math/PlanePath/PeanoCurve.pm

Criterion	Covered	Total	%
statement	174	198	87.8
branch	38	64	59.3
condition	22	28	78.5
subroutine	19	26	73.0
pod	7	7	100.0
total	260	323	80.5

# Copyright 2010, 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 .


# cf
#
# http://www.cut-the-knot.org/Curriculum/Geometry/PeanoComplete.shtml
#     applet, directions in 9 sub-parts
#
# math-image --path=PeanoCurve,radix=5 --all --output=numbers
# math-image --path=PeanoCurve,radix=5 --lines
#
# -----------
# Peano:
# T = 0.a1 a2 a3 a4 ...
#        x  y  x  y
#
# X = 0.b1        b2 ...
#       a1        a3.k(a2)
#
# Y = 0.c1        c2 ...
#       a2.k(a1)  a4.k(a1,a3)
#
# b1=a1
# c1 = a2 comp(a1)
# b2 = a3 comp(a2)
# c2 = a4 comp(a1+a3)
#
# bn = a[2n-1] comp a2+a4+...+a[2n-2]
# cn = a[2n] comp a1+a3+...+a[2n-1]
#
# Brouwer(?) no continuous one-to-one between R and RxR, so line and plane
# are distinguished.
#


package Math::PlanePath::PeanoCurve;
use 5.004;
use strict;
#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',
  'digit_join_lowtohigh';
use Math::PlanePath::Base::NSEW;

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


use constant n_start => 0;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
*xy_is_visited = \&Math::PlanePath::Base::Generic::xy_is_visited_quad1;

use constant parameter_info_array =>
  [ { name      => 'radix',
      display   => 'Radix',
      share_key => 'radix_3',
      type      => 'integer',
      minimum   => 2,
      default   => 3,
      width     => 3,
    } ];

# shared by WunderlichSerpentine
sub dx_minimum {
  my ($self) = @_;
  return ($self->{'radix'} % 2
          ? -1      # odd
          : undef); # even, unlimited
}
sub dx_maximum {
  my ($self) = @_;
  return ($self->{'radix'} % 2
          ? 1         # odd
          : undef);   # even, unlimited
}

# shared by WunderlichSerpentine
sub _UNDOCUMENTED__dxdy_list {
  my ($self) = @_;
  return ($self->{'radix'} % 2
          ? Math::PlanePath::Base::NSEW->_UNDOCUMENTED__dxdy_list
          : ());   # even, unlimited
}
#  *---  b^2-1 -- b^2 ---- b^2+b-1 = (b+1)b-1
#  |                          |
#  *-------
#          |
#  0 ----- b
#
sub _UNDOCUMENTED__dxdy_list_at_n {
  my ($self) = @_;
  return ($self->{'radix'} + 1) * $self->{'radix'} - 1;
}

# shared by WunderlichSerpentine
*dy_minimum = \&dx_minimum;
*dy_maximum = \&dx_maximum;

*dsumxy_minimum = \&dx_minimum;
*dsumxy_maximum = \&dx_maximum;

*ddiffxy_minimum = \&dx_minimum;
*ddiffxy_maximum = \&dx_maximum;

sub dir_maximum_dxdy {
  my ($self) = @_;
  return ($self->{'radix'} % 2
          ? (0,-1)   # odd, South
          : (0,0));  # even, supremum
}

sub _UNDOCUMENTED__turn_any_left_at_n {
  my ($self) = @_;
  return $self->{'radix'} - 1;
}
sub _UNDOCUMENTED__turn_any_right_at_n {
  my ($self) = @_;
  return ($self->{'radix'} == 2 ? 5
          : 2*$self->{'radix'} - 1);
}


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

sub new {
  my $self = shift->SUPER::new(@_);

  if (! $self->{'radix'} || $self->{'radix'} < 2) {
    $self->{'radix'} = 3;
  }
  return $self;
}

sub _n_to_xykk {
  my ($self, $n) = @_;
  my $radix = $self->{'radix'};
  my $radix_minus_1 = $radix - 1;

  my @ndigits = digit_split_lowtohigh($n,$radix);
  if (scalar(@ndigits) & 1) {
    push @ndigits, 0;            # so even number of entries
  }
  ### @ndigits

  my $xk = 0;
  my $yk = 0;
  my @ydigits;
  my @xdigits;

  for (my $i = $#ndigits >> 1; @ndigits; $i--) {    # high to low
    ### $i
    {
      my $ndigit = pop @ndigits;  # high to low
      $xk ^= $ndigit;
      $ydigits[$i] = ($yk & 1 ? $radix_minus_1-$ndigit : $ndigit);
    }
    {
      my $ndigit = pop @ndigits;
      $yk ^= $ndigit;
      $xdigits[$i] = ($xk & 1 ? $radix_minus_1-$ndigit : $ndigit);
    }
  }
  my $zero = $n*0;
  return ((map {digit_join_lowtohigh($_, $radix, $zero)} \@xdigits, \@ydigits),
          $xk,$yk);
}


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

  {
    # ENHANCE-ME: for odd radix the ends join and the direction can be had
    # without a 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 ($x,$y) = _n_to_xykk($self,$n);
  return ($x,$y);
}

sub _xykk_to_n {
  my ($self, $x,$y, $offset_xk,$offset_yk) = @_;
  ### PeanoCurve _xykk_to_n(): "$x, $y offset $offset_xk,$offset_yk"

  if (($offset_xk    && ($x-=$offset_xk) < 0)
      || ($offset_yk && ($y-=$offset_yk) < 0)) {
    return;    # offset goes negative
  }

  my $radix = $self->{'radix'};
  my @x = digit_split_lowtohigh ($x, $radix);
  my @y = digit_split_lowtohigh ($y, $radix);

  my $radix_minus_1 = $radix - 1;
  my $xk = 0;
  my $yk = 0;

  my @n;  # stored low to high, generated from high to low
  my $i_high = max($#x,$#y);
  my $npos = 2*$i_high+1;

  foreach my $i (reverse 0 .. $i_high) {  # high to low
    {
      my $digit = $y[$i] || 0;
      if ($yk & 1) {
        $digit = $radix_minus_1 - $digit;  # reverse digit
      }
      $n[$npos--] = $digit;
      $xk ^= $digit;
    }
    {
      my $digit = $x[$i] || 0;
      if ($xk & 1) {
        $digit = $radix_minus_1 - $digit;  # reverse digit
      }
      $n[$npos--] = $digit;
      $yk ^= $digit;
    }
  }
  ### final n: @n
  ### final xkyk: ($xk&1).' '.($yk&1)
  return ((! defined $offset_xk    || ($xk&1) == $offset_xk)
          && (! defined $offset_yk || ($yk&1) == $offset_yk)
          ? (digit_join_lowtohigh (\@n, $radix,
                                   $x*0*$y))  # inherit bignum 0
          : ());
}

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

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

  if ($x < 0 || $y < 0) { return undef; }
  if (is_infinite($x)) { return $x; }
  if (is_infinite($y)) { return $y; }

  return _xykk_to_n($self, $x,$y);
}

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

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

  my $radix = $self->{'radix'};

  my ($power, $level) = round_down_pow (max($x2,$y2), $radix);
  if (is_infinite($level)) {
    return (0, $level);
  }

  my $n_power = $power * $power * $radix;
  my $max_x = 0;
  my $max_y = 0;
  my $max_n = 0;
  my $max_xk = 0;
  my $max_yk = 0;

  my $min_x = 0;
  my $min_y = 0;
  my $min_n = 0;
  my $min_xk = 0;
  my $min_yk = 0;

  # l<=c
  #     l>c2 or h-1
  #     l>c2 or h<=c1
  # so does overlap if
  #     l<=c2 and h>c1
  #
  my $radix_minus_1 = $radix - 1;
  my $overlap = sub {
    my ($c,$ck,$digit, $c1,$c2) = @_;
    if ($ck & 1) {
      $digit = $radix_minus_1 - $digit;
    }
    ### overlap consider: "inv".($ck&1)."digit=$digit ".($c+$digit*$power)."<=c<".($c+($digit+1)*$power)." cf $c1 to $c2 incl"
    return ($c + $digit*$power <= $c2
            && $c + ($digit+1)*$power > $c1);
  };

  while ($level-- >= 0) {
    ### $power
    ### $n_power
    ### $max_n
    ### $min_n
    {
      my $digit;
      for ($digit = $radix_minus_1; $digit > 0; $digit--) {
        last if &$overlap ($max_y,$max_yk,$digit, $y1,$y2);
      }
      $max_n += $n_power * $digit;
      $max_xk ^= $digit;
      if ($max_yk&1) { $digit = $radix_minus_1 - $digit; }
      $max_y += $power * $digit;
      ### max y digit (complemented): $digit
      ### $max_y
      ### $max_n
    }
    {
      my $digit;
      for ($digit = 0; $digit < $radix_minus_1; $digit++) {
        last if &$overlap ($min_y,$min_yk,$digit, $y1,$y2);
      }
      $min_n += $n_power * $digit;
      $min_xk ^= $digit;
      if ($min_yk&1) { $digit = $radix_minus_1 - $digit; }
      $min_y += $power * $digit;
      ### min y digit (complemented): $digit
      ### $min_y
      ### $min_n
    }

    $n_power = int($n_power/$radix);
    {
      my $digit;
      for ($digit = $radix_minus_1; $digit > 0; $digit--) {
        last if &$overlap ($max_x,$max_xk,$digit, $x1,$x2);
      }
      $max_n += $n_power * $digit;
      $max_yk ^= $digit;
      if ($max_xk&1) { $digit = $radix_minus_1 - $digit; }
      $max_x += $power * $digit;
      ### max x digit (complemented): $digit
      ### $max_x
      ### $max_n
    }
    {
      my $digit;
      for ($digit = 0; $digit < $radix_minus_1; $digit++) {
        last if &$overlap ($min_x,$min_xk,$digit, $x1,$x2);
      }
      $min_n += $n_power * $digit;
      $min_yk ^= $digit;
      if ($min_xk&1) { $digit = $radix_minus_1 - $digit; }
      $min_x += $power * $digit;
      ### min x digit (complemented): $digit
      ### $min_x
      ### $min_n
    }

    $power = int($power/$radix);
    $n_power = int($n_power/$radix);
  }
  ### is: "$min_n at $min_x,$min_y  to  $max_n at $max_x,$max_y"
  return ($min_n, $max_n);
}

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

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

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