File Coverage

blib/lib/Math/PlanePath/RationalsTree.pm

Criterion	Covered	Total	%
statement	176	222	79.2
branch	62	90	68.8
condition	11	28	39.2
subroutine	26	36	72.2
pod	18	18	100.0
total	293	394	74.3

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


# SB,CW N with same X,Y is those N which are palindromes below high 1-bit
# as noted Claudio Bonanno and Stefano Isola, ``Orderings of the Rationals
# and Dynamical Systems'', May 16, 2008.
# cf A006995 binary palindromes, so always odd
#    A178225 characteristic of binary palindromes
#    A048700 binary palindromes odd length
#    A048701 binary palindromes even length
# A044051 binary palindromes (B+1)/2, B odd so B+1 even
# A044051-1 = (B-1)/2 strips low 1-bit to be palindromes below high 1-bit

# Boyko B. Bantchev, "Fraction Space Revisited"
# http://www.math.bas.bg/bantchev/articles/fractions.pdf

# cf Ronald L. Graham, Donald E. Knuth, and Oren Patashnik, Concrete
# Mathematics: A Foundation for Computer Science, Second
# Edition. Addison-Wesley. 1994.
# On Stern-Brocot tree.

# cf A054429 permutation reverse within binary row
#    A057114 - permutation SB X -> X+1
#    A057115 - permutation SB X -> X-1
#

#                    high-to-low   low-to-high
# (X+Y)/Y  Y/(X+Y)     HCS            AYT
# X/(X+Y)  (X+Y)/Y      CW            SB    \ alt bit flips
# Y/(X+Y)  (X+Y)/X     Drib          Bird   /
#
#     9  10                    12  10
# 8      11                 8      14
#        12  13                     9  13
#            14                        11
#            15                        15
#
# Stern-Brocot              Calkin-Wilf


#------------------------------------------------------------------------------
# HCS turn left when even number of 1-bits in N+1
#     turn right when odd number of 1-bits in N+1
#
# A010059 start=0: 1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0
#   match 1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0
#   PlanePathTurn planepath=RationalsTree,tree_type=HCS,  turn_type=Left
#
# A010060 start=0: 0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,0,1,1,0,1
#   match 0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,0,1,1,0,1
#   PlanePathTurn planepath=RationalsTree,tree_type=HCS,  turn_type=Right
#
# 10  |     768        50                  58       896
#  9  |     384   49        52   60        57  448       640
#  8  |     192        27        31       224       320
#  7  |      96   25   26   30   29  112       160   41   42
#  6  |      48                  56        80
#  5  |      24   13   15   28        40   21   23   44
#  4  |      12        14        20        22        36
#  3  |       6    7        10   11        18   19        34
#  2  |       3         5         9        17        33
#  1  |       1    2    4    8   16   32   64  128  256  512
# Y=0 |
#     +-----------------------------------------------------
#       X=0   1    2    3    4    5    6    7    8    9   10
#
#                               1/1
#                  /------------- -------------\
#               2/1                             1/2               2,3 L,R
#          /----   ----\                   /----   ----\
#       3/1             3/2             1/3             2/3    4,5,6,7 L,L,R,R
#      /   \           /   \           /   \           /   \      8        12
#   4/1     5/2     4/3     5/3     1/4     2/5     3/4     3/5   L,L,R,L, R,R,L,R
#  /   \   /   \   /   \   /   \   /   \   /   \   /   \   /   \
# 5/1 7/2 7/3 8/3 5/4 7/5 7/4 8/5 1/5 2/7 3/7 3/8 4/5 5/7 4/7 5/8
#
#         *
#        / \                        U=0 = X+Y, Y           shear
#       /   *                       D=1 = Y,   X+Y         shear+transpose
#      /     \a  = 0.1^k.1
#     N
#      \     /b  = 1.0^k.0
#       \   *
#        \ / \c  = 1.0^k.1         c=even bits, left
#         *
#
# F[-1]=1 F[0]=0 F[1]=1 F[2]=1 F[3]=2 F[4]=3 F[5]=5 ...
# 1^k is F[k-1]*X+F[k]*Y, F[k]*X+F[k+1]*Y
#  X   ,     Y    0
#     Y,  X+ Y    1
#  X+ Y,  X+2Y    2
#  X+2Y, 2X+3Y    3
# 2X+3Y, 3X+5Y    4
#
# then aX = F[k]*X+F[k+1]*Y + F[k+1]*X+F[k+2]*Y
#         = (F[k]+F[k+1])*X + (F[k+1]+F[k+2])*Y
#         = F[k+2]*X + F[k+3]*Y
#      aY = F[k+1]*X + F[k+2]*Y                 near X=phi*Y big
#
# 0^k is X+k*Y, Y
# so bX = Y
#    bY = X+k*Y + Y = X+(k+1)*Y                 near Y axis
#
# c1X = Y
# c1Y = X+Y
# c2X = Y + k*(X+Y) = k*X + (k+1)*Y
# c2Y = X+Y
# cX = X+Y
# cY = k*X + (k+1)*Y + X+Y = (k+1)X + (k+2)Y    near X=Y

#
#         *
#        / \ /a  = 0.1^k.0
#       /   *
#      /     \b  = 0.1^k.1
#     N
#      \     /c  = 1.0^k.1     c=even bits, left
#       \   *
#        \ /
#         *
#------------------------------------------------------------------------------

package Math::PlanePath::RationalsTree;
use 5.004;
use strict;
use Carp 'croak';
#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',
  'bit_split_lowtohigh',
  'digit_join_lowtohigh';
*_divrem = \&Math::PlanePath::_divrem;

use Math::PlanePath::CoprimeColumns;
*_coprime = \&Math::PlanePath::CoprimeColumns::_coprime;

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


use constant parameter_info_array =>
  [ { name            => 'tree_type',
      share_key       => 'tree_type_rationalstree',
      display         => 'Tree Type',
      type            => 'enum',
      default         => 'SB',
      choices         => ['SB','CW','AYT','HCS','Bird','Drib','L'],
      choices_display => ['SB','CW','AYT','HCS','Bird','Drib','L'],
    },
  ];

use constant class_x_negative => 0;
use constant class_y_negative => 0;
sub x_minimum {
  my ($self) = @_;
  return ($self->{'tree_type'} eq 'L' ? 0 : 1);
}
use constant y_minimum => 1;
use constant gcdxy_maximum => 1;  # no common factor
use constant tree_num_children_list => (2); # complete binary tree
use constant tree_n_to_subheight => undef; # complete tree, all infinity

{
  my %absdy_minimum = (# SB   => 0,
                       CW   => 1,
                       # AYT  => 0,
                       # Bird => 0,
                       # Drib => 0,
                       L    => 1);
  sub absdy_minimum {
    my ($self) = @_;
    return $absdy_minimum{$self->{'tree_type'}} || 0;
  }
}

{
  # Drib apparent minimum dX=k dY=2*k+1 approaches dX=1,dY=2
  my %dir_minimum_dxdy = (CW   => [0,1],
                          Drib => [1,2],
                          L    => [1,1], # at N=0 dX=1,dY=1
                         );
  sub dir_minimum_dxdy {
    my ($self) = @_;
    return @{$dir_minimum_dxdy{$self->{'tree_type'}} || [1,0]};
  }
}
{
  my %dir_maximum_dxdy
    = (SB   => [1,-1],
       # CW   => [0,0],

       Bird => [1,-1],
       # Drib => [0,0],

       HCS  => [2,-1],
       # AYT  => [0,0],
       # L    => [0,0], # at 2^k-1 dX=k+1,dY=-1 so approach Dir=4
      );
  sub dir_maximum_dxdy {
    my ($self) = @_;
    return @{$dir_maximum_dxdy{$self->{'tree_type'}} || [0,0]};
  }
}

{
  my %turn_any_straight = (# SB   => 0,
                           # CW   => 0,
                           Bird => 1,   # straight at N=7 and N=8
                           # Drib => 0,
                           AYT  => 1,   # straight at N=7
                           # HCS  => 0,
                           # L    => 0,
                          );
  sub turn_any_straight {
    my ($self) = @_;
    return $turn_any_straight{$self->{'tree_type'}};
  }
}


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

my %attributes = (CW   => [ n_start => 1, ],
                  SB   => [ n_start => 1, reverse_bits => 1 ],
                  Drib => [ n_start => 1, alternating => 1 ],
                  Bird => [ n_start => 1, alternating => 1, reverse_bits => 1 ],
                  AYT  => [ n_start => 1, sep1s => 1 ],
                  HCS  => [ n_start => 1, sep1s => 1, reverse_bits => 1 ],
                  L    => [ n_start => 0 ],
                 );

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

  my $tree_type = ($self->{'tree_type'} ||= 'SB');
  my $attributes = $attributes{$tree_type}
    || croak "Unrecognised tree type: ",$tree_type;
  %$self = (%$self, @$attributes);

  return $self;
}

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

  if ($n < $self->{'n_start'}) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  # what to do for fractional $n?
  {
    my $int = int($n);
    if ($n != $int) {
      ### frac ...
      my $frac = $n - $int;  # inherit possible BigFloat/BigRat
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+1);
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      ### x1,y1: "$x1, $y1"
      ### x2,y2: "$x2, $y2"
      ### dx,dy: "$dx, $dy"
      ### result: ($frac*$dx + $x1).', '.($frac*$dy + $y1)
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int;
  }

  my $zero = ($n * 0);  # inherit bignum 0
  my $one = $zero + 1;  # inherit bignum 1

  if ($self->{'n_start'} == 0) {
    # L tree adjust;
    $n += 2;
  }

  my @nbits = bit_split_lowtohigh($n);
  pop @nbits;
  ### lowtohigh sans high: @nbits

  if (! $self->{'reverse_bits'}) {
    @nbits = reverse @nbits;
    ### reverse to: @nbits
  }

  my $x = $one;
  my $y = $one;

  if ($self->{'sep1s'}) {
    foreach my $nbit (@nbits) {
      $x += $y;
      if ($nbit) {
        ($x,$y) = ($y,$x);
      }
    }

  } elsif ($self->{'alternating'}) {
    foreach my $nbit (@nbits) {
      ($x,$y) = ($y,$x);
      if ($nbit) {
        $x += $y;     # (x,y) -> (x+y,x), including swap
      } else {
        $y += $x;     # (x,y) -> (y,x+y), including swap
      }
    }

  } elsif ($self->{'tree_type'} eq 'L') {
    my $sub = 2;
    foreach my $nbit (@nbits) {
      if ($nbit) {
        $y += $x;     # (x,y) -> (x,x+y)
        $sub = 0;
      } else {
        $x += $y;     # (x,y) -> (x+y,y)
      }
    }
    $x -= $sub;   # -2 at N=00...000 all zero bits

  } else {
    ### nbits apply CW: @nbits
    foreach my $nbit (@nbits) {   # high to low
      if ($nbit) {
        $x += $y;     # (x,y) -> (x+y,y)
      } else {
        $y += $x;     # (x,y) -> (x,x+y)
      }
    }
  }
  ### result: "$x, $y"
  return ($x,$y);
}

sub xy_is_visited {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  if ($self->{'tree_type'} eq 'L' && $x == 0 && $y == 1) {
    return 1;
  }
  if ($x < 1
      || $y < 1
      || ! _coprime($x,$y)) {
    return 0;
  }
  return 1;
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### RationalsTree xy_to_n(): "$x,$y   $self->{'tree_type'}"

  if ($x < $self->{'n_start'} || $y < 1) {
    return undef;
  }
  if (is_infinite($x)) {
    return $x;
  }
  if (is_infinite($y)) {
    return $y;
  }

  my @quotients = _xy_to_quotients($x,$y)
    or return undef;  # $x,$y have a common factor
  ### @quotients

  my @nbits;
  if ($self->{'sep1s'}) {
    $quotients[0]++;  # the integer part, making it 1 or more
    foreach my $q (@quotients) {
      push @nbits, (0) x ($q-1), 1;   # runs of "000..0001"
    }
    pop @nbits;  # no high 1-bit separator

  } else {
    if ($quotients[0] < 0) {   # X=0,Y=1 in tree_type="L"
      return $self->{'n_start'};
    }

    my $bit = 1;
    foreach my $q (@quotients) {
      push @nbits, ($bit) x $q;
      $bit ^= 1;     # alternate runs of "00000" or "11111"
    }
    ### nbits in quotient order: @nbits

    if ($self->{'alternating'}) {
      # Flip every second bit, starting from the second lowest.
      for (my $i = 1; $i <= $#nbits; $i += 2) {
        $nbits[$i] ^= 1;
      }
    }

    if ($self->{'tree_type'} eq 'L') {
      # Flip all bits.
      my $anyones = 0;
      foreach my $nbit (@nbits) {
        $nbit ^= 1;   # mutate array
        $anyones ||= $nbit;
      }
      unless ($anyones) {
        push @nbits, 0,0;
      }
    }
  }

  if ($self->{'reverse_bits'}) {
    @nbits = reverse @nbits;
  }
  push @nbits, 1;   # high 1-bit

  ### @nbits
  my $n = digit_join_lowtohigh (\@nbits, 2,
                                $x*0*$y);   # inherit bignum 0
  if ($self->{'tree_type'} eq 'L') {
    return $n-2;
  } else {
    return $n;
  }
}

# Return a list of the quotients from Euclid's greatest common divisor
# algorithm on X,Y.  This is also the terms of the continued fraction
# expansion of rational X/Y.
#
# The last term, the last in the list, is decremented since this is what the
# code above requires.  This term is the top-most quotient in for example
# gcd(7,1) is 7=7*1+0 with q=7 returned as 6.
#
# If $x,$y have a common factor then the return is an empty list.
# If $x,$y have no common factor then the returned list is always one or
# more quotients.
#
sub _xy_to_quotients {
  my ($x,$y) = @_;
  my @ret;
  for (;;) {
    my ($q, $r) = _divrem($x,$y);
    push @ret, $q;
    last unless $r;
    $x = $y;
    $y = $r;
  }

  if ($y > 1) {
    ### found Y>1 common factor, no N at this X,Y ...
    return;
  }
  $ret[-1]--;
  return @ret;
}


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

  $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;
  ### $x2
  ### $y2

  if ($x2 < 1 || $y2 < 1) {
    ### no values, rect below first quadrant
    if ($self->{'n_start'}) {
      return (1,0);
    } else {
      return (0,0);
    }
  }

  my $zero = ($x1 * 0 * $y1 * $x2 * $y2);  # inherit bignum
  ### $zero

  if ($x1 < 1) { $x1 = 1; }
  if ($y1 < 1) { $y1 = 1; }

  # # big x2, small y1
  # # big y2, small x1
  # my $level = _bingcd_max ($y2,$x1);
  # ### $level
  # {
  #   my $l2 = _bingcd_max ($x2,$y1);
  #   ### $l2
  #   if ($l2 > $level) { $level = $l2; }
  # }

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

  return ($self->{'n_start'},
          $self->{'n_start'} + (2+$zero) ** ($level + 3));
}

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

  if ($x < $y) { ($x,$y) = ($y,$x) }

  ### div: int($x/$y)
  ### bingcd: int($x/$y) + $y

  return int($x/$y) + $y + 1;
}

#   ### fib: _fib_log($y)
# # ENHANCE-ME: log base PHI, or something close for BigInt
# # 2*log2() means log base sqrt(2)=1.4 instead of PHI=1.6
# #
# # use constant 1.02; # for leading underscore
# # use constant _PHI => (1 + sqrt(5)) / 2;
# #
# sub _fib_log {
#   my ($x) = @_;
#   ### _fib_log(): $x
#   my $f0 = ($x * 0);
#   my $f1 = $f0 + 1;
#   my $count = 0;
#   while ($x > $f0) {
#     $count++;
#     ($f0,$f1) = ($f1,$f0+$f1);
#   }
#   return $count;
# }

#------------------------------------------------------------------------------
use constant tree_num_roots => 1;

# N=1 basis children 2N,2N+1
# N=S basis 2(N-(S-1))+(S-1)
#           = 2N - 2(S-1) + (S-1)
#           = 2N - (S-1)
sub tree_n_children {
  my ($self, $n) = @_;
  my $n_start = $self->{'n_start'};
  if ($n >= $n_start) {
    $n = 2*$n - $n_start;
    return ($n+1, $n+2);
  } else {
    return;
  }
}
sub tree_n_num_children {
  my ($self, $n) = @_;
  return ($n >= $self->{'n_start'} ? 2 : undef);
}
sub tree_n_parent {
  my ($self, $n) = @_;
  $n = $n - $self->{'n_start'}; # N=0 basis, and warn if $n==undef
  if ($n > 0) {
    return int(($n-1)/2) + $self->{'n_start'};
  } else {
    return undef;
  }
}
sub tree_n_to_depth {
  my ($self, $n) = @_;
  ### RationalsTree tree_n_to_depth(): $n
  $n = $n - $self->{'n_start'}; # N=0 basis, and warn if $n==undef
  unless ($n >= 0) {
    return undef;
  }
  my ($pow, $exp) = round_down_pow ($n+1, 2);
  return $exp;
}

sub tree_depth_to_n {
  my ($self, $depth) = @_;
  return ($depth >= 0
          ? 2**$depth + $self->{'n_start'}-1
          : undef);
}
# (2^(d+1)+s-1)-1 = 2^(d+1)+s-2
sub tree_depth_to_n_end {
  my ($self, $depth) = @_;
  return ($depth >= 0
          ? 2**($depth+1) + $self->{'n_start'}-2
          : undef);
}
sub tree_depth_to_n_range {
  my ($self, $depth) = @_;
  if ($depth >= 0) {
    my $pow = 2**$depth;
    return ($pow + $self->{'n_start'}-1, 2*$pow + $self->{'n_start'}-2);
  }
  return; # no such $depth
}
sub tree_depth_to_width {
  my ($self, $depth) = @_;
  return ($depth >= 0
          ? 2**$depth
          : undef);
}

1;
__END__