File Coverage

blib/lib/Math/PlanePath/CfracDigits.pm

Criterion	Covered	Total	%
statement	124	150	82.6
branch	24	40	60.0
condition	5	18	27.7
subroutine	25	30	83.3
pod	7	7	100.0
total	185	245	75.5

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


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

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::RationalsTree;
*_xy_to_quotients = \&Math::PlanePath::RationalsTree::_xy_to_quotients;

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

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


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

use constant n_start => 0;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant x_minimum => 1;
use constant y_minimum => 2;
use constant diffxy_maximum => -1; # upper octant X<=Y-1 so X-Y<=-1
use constant gcdxy_maximum => 1;  # no common factor

# FIXME: believe this is right, but check N+1 always changes Y
sub absdy_minimum {
  my ($self) = @_;
  return ($self->{'radix'} < 3 ? 0 : 1);
}

# radix=1 N=1       has dir4=0
# radix=2 N=5628    has dir4=0 dx=9,dy=0
# radix=3 N=1189140 has dir4=0 dx=1,dy=0
# radix=4 N=169405  has dir4=0 dx=2,dy=0
# always eventually 0 ?
# use constant dir_minimum_dxdy => (1,0);  # the default

# radix=1 N=4    dX=1,dY=-1 for dir4=3.5
# radix=2 N=4413 dX=9,dY=-1
# radix=3 N=9492 dX=3,dY=-1
# ENHANCE-ME: suspect believe approaches 360 degrees, eventually, but proof?
# use constant dir_maximum_dxdy => (0,0);  # the default

sub turn_any_straight {
  my ($self) = @_;
  return ($self->{'radix'} != 1);  # radix=1 never straight
}

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'};
}


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

sub new {
  my $self = shift->SUPER::new (@_);
  unless ($self->{'radix'} && $self->{'radix'} >= 1) {
    $self->{'radix'} = 2;
  }
  return $self;
}

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

  if ($n < 0) { return; }
  if (is_infinite($n)) { return ($n,$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 $radix = $self->{'radix'};
  my $zero = ($n * 0);  # inherit bignum 0
  my $x = $zero;
  my $y = 1 + $zero;    # inherit bignum 1

  foreach my $q (_n_to_quotients_bottomtotop($n,$radix,$zero)) {  # bottom to top
    ### at: "$x,$y   q=$q"

    # 1/(q + X/Y) = 1/((qY+X)/Y)
    #             = Y/(qY+X)
    ($x,$y) = ($y, $q*$y + $x);
  }

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

# Return a list of quotients bottom to top.  The base3 digits of N are split
# by "3" delimiters and the parts adjusted so the first bottom-most q>=2 and
# the rest q>=1.  The values are ready to be used as continued fraction
# terms.
#
sub _n_to_quotients_bottomtotop {
  my ($n, $radix, $zero) = @_;
  ### _n_to_quotients_bottomtotop(): $n

  my $radix_plus_1 = $radix + 1;
  my @ret;
  my @group;
  foreach my $digit (_digit_split_1toR_lowtohigh($n,$radix_plus_1)) {
    if ($digit == $radix_plus_1) {
      ### @group
      push @ret, _digit_join_1toR_destructive(\@group, $radix, $zero) + 1;
      @group = ();
    } else {
      push @group, $digit;
    }
  }
  ### final group: @group
  push @ret, _digit_join_1toR_destructive(\@group, $radix, $zero) + 1;

  $ret[0] += 1;  # bottom-most is +2 rather than +1

  ### _n_to_quotients_bottomtotop result: @ret
  return @ret;
}

# Return a list of digits 1 <= d <= R which is $n written in $radix, low to
# high digits.
sub _digit_split_1toR_lowtohigh {
  my ($n, $radix) = @_;
  ### assert: $radix >= 1
  ### assert: $n >= 0

  if ($radix == 1) {
    return (1) x $n;
  }
  my @digits = digit_split_lowtohigh($n,$radix);

  # mangle 0 -> R
  my $borrow = 0;
  foreach my $digit (@digits) {   # low to high
    if ($borrow = (($digit -= $borrow) <= 0)) {  # modify array contents
      $digit += $radix;
    }
  }
  if ($borrow) {
    ### assert: $digits[-1] == $radix
    pop @digits;
  }

  return @digits;
}

sub _digit_join_1toR_destructive {
  my ($aref, $radix, $zero) = @_;
  ### assert: $radix >= 1

  if ($radix == 1) {
    return scalar(@$aref);
  }

  # mangle any digit==$radix down to digit=0
  my $carry = 0;
  foreach my $digit (@$aref) {   # low to high
    if ($carry = (($digit += $carry) >= $radix)) {  # modify array contents
      $digit -= $radix;
    }
  }
  if ($carry) {
    push @$aref, 1;
  }

  ### _digit_join_1toR_destructive() result: digit_join_lowtohigh($aref, $radix, $zero)
  return digit_join_lowtohigh($aref, $radix, $zero);
}

sub xy_is_visited {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  return (! ($x < 1 || $y < 2 || $x >= $y)
          && _coprime($x,$y));
}

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

  if (is_infinite($x)) { return $x; }
  if (is_infinite($y)) { return $y; }
  if ($x < 1 || $y < 2 || $x >= $y) {
    return undef;
  }

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

  # drop initial 0 integer part
  ### assert: $quotients[0] == 0
  shift @quotients;

  return _cfrac_join_toptobottom(\@quotients,
                                 $self->{'radix'},
                                 $x * 0 * $y);   # inherit bignum 0
}

# $aref is a list of continued fraction quotients from top-most to
# bottom-most.  There's no initial integer term in $aref.  Each quotient is
# q >= 1 except the bottom-most which q-1 and so also >=1.
#
sub _cfrac_join_toptobottom {
  my ($aref, $radix, $zero) = @_;
  ### _cfrac_join_toptobottom(): $aref

  my @digits;
  foreach my $q (reverse @$aref) {
    ### assert: $q >= 1
    push @digits, _digit_split_1toR_lowtohigh($q-1, $radix), $radix+1;
  }
  pop @digits;  # no high delimiter
  ### groups digits 1toR: @digits
  return _digit_join_1toR_destructive(\@digits, $radix+1, $zero);
}


# X/Y = F[k]/F[k+1] quotients all 1
# N = all delimiter digits R,R,...,R
#   = 1222...2221
#   = R^k + 2*(R^k+1)/(R-1) - 1
#   = (RR^k - R^k + 2R^k + 2 - R + 1) / (R-1)
#   = (RR^k + R^k - R + 3) / (R-1)
#   = ((R+1)R^k - R + 3) / (R-1)
# take high as "12" = R+2
# k = log(Y)/log(phi)
# N = (R+2) * R ** k
# N = Y ** (log(R)/log(phi))
#
# 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

  #   |    /
  #   |   / x1
  #   |  /  +-----y2
  #   | /   |
  #   |/    +-----
  #
  if ($x2 < 1 || $y2 < 2 || $x1 >= $y2) {
    ### no values, rect outside upper octant ...
    return (1,0);
  }

  my $zero = ($x1 * 0 * $y1 * $x2 * $y2);  # inherit bignum
  my $radix = $self->{'radix'};

  return (0,
          ($radix+3)
          * ($radix+1 + $zero) ** ($radix == 1
                                   ? $y2
                                   : _log_phi_estimate($y2,$radix)));
}

# Return an estimate of log base phi of $x, that being log($x)/log(phi),
# where phi=(1+sqrt(5))/2 the golden ratio.
#
sub _log_phi_estimate {
  my ($x) = @_;
  my ($pow,$exp) = round_down_pow ($x, 2);
  return int ($exp * (log(2) / log((1+sqrt(5))/2)));
}

1;
__END__