File Coverage

blib/lib/Math/PlanePath/WythoffArray.pm

Criterion	Covered	Total	%
statement	43	122	35.2
branch	0	22	0.0
condition	4	16	25.0
subroutine	15	19	78.9
pod	7	7	100.0
total	69	186	37.1

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


# Classic Sequences
# http://oeis.org/classic.html
#
# Clark Kimberling
# http://faculty.evansville.edu/ck6/integer/intersp.html
#
# cf A175004 similar to wythoff but rows recurrence
#            r(n-1)+r(n-2)+1 extra +1 in each step
#            floor(n*phi+2/phi)
#
# cf Stolarsky round_nearest(n*phi)
#    A035506 stolarsky by diagonals
#    A035507   inverse
#    A007067 stolarsky first column

# Maybe:
# my ($x,$y) = $path->pair_to_xy($a,$b);
# Return the $x,$y which has ($a,$b).
# Advance $a,$b if before start of row.

# Carlitz and Hoggatt "Fibonacci Representations", Fibonacci Quarterly,
# volume 10, number 1, January 1972
# http://www.fq.math.ca/10-1.html
# http://www.fq.math.ca/Scanned/10-1/carlitz1.pdf


package Math::PlanePath::WythoffArray;
use 5.004;
use strict;
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

use vars '$VERSION', '@ISA';
$VERSION = 128;
use Math::PlanePath;
*_sqrtint = \&Math::PlanePath::_sqrtint;
@ISA = ('Math::PlanePath');

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'bit_split_lowtohigh';

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


use constant parameter_info_array =>
  [ { name        => 'x_start',
      display     => 'X start',
      type        => 'integer',
      default     => 0,
      width       => 3,
      description => 'Starting X coordinate.',
    },
    { name        => 'y_start',
      display     => 'Y start',
      type        => 'integer',
      default     => 0,
      width       => 3,
      description => 'Starting Y coordinate.',
    },
  ];

use constant default_n_start => 1;
use constant class_x_negative => 0;
use constant class_y_negative => 0;

sub x_minimum {
  my ($self) = @_;
  return $self->{'x_start'};
}
sub y_minimum {
  my ($self) = @_;
  return $self->{'y_start'};
}
use constant absdx_minimum => 1;
use constant dir_maximum_dxdy => (3,-1);  # N=4 to N=5 dX=3,dY=-1

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

sub new {
  my $self = shift->SUPER::new(@_);
  $self->{'x_start'} ||= 0;
  $self->{'y_start'} ||= 0;
  return $self;
}

sub xy_is_visited {
  my ($self, $x, $y) = @_;
  return ((round_nearest($x) >= $self->{'x_start'})
          && (round_nearest($y) >= $self->{'y_start'}));
}

#------------------------------------------------------------------------------
#   4  |  12   20   32   52   84  136  220  356  576  932 1508
#   3  |   9   15   24   39   63  102  165  267  432  699 1131
#   2  |   6   10   16   26   42   68  110  178  288  466  754
#   1  |   4    7   11   18   29   47   76  123  199  322  521
# Y=0  |   1    2    3    5    8   13   21   34   55   89  144
#      +-------------------------------------------------------
#        X=0    1    2    3    4    5    6    7    8    9   10
# 13,8,5,3,2,1
# 4 = 3+1     -> 1
# 6 = 5+1     -> 2
# 9 = 8+1     -> 3
# 12 = 8+3+1  -> 3+1=4
# 14 = 13+1   -> 5

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

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

  {
    # fractions on straight line between integer points
    my $int = int($n);
    if ($n != $int) {
      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;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int;
  }

  # f1+f0 > i
  # f0 > i-f1
  # check i-f1 as the stopping point, so that if i=UV_MAX then won't
  # overflow a UV trying to get to f1>=i
  #
  my @fibs;
  {
    my $f0 = ($n * 0);  # inherit bignum 0
    my $f1 = $f0 + 1;   # inherit bignum 1
    while ($f0 <= $n-$f1) {
      ($f1,$f0) = ($f1+$f0,$f1);
      push @fibs, $f1;      # starting $fibs[0]=1
    }
  }
  ### @fibs

  # indices into fib[] which are the Fibonaccis adding up to $n
  my @indices;
  for (my $i = $#fibs; $i >= 0; $i--) {
    ### at: "n=$n f=".$fibs[$i]
    if ($n >= $fibs[$i]) {
      push @indices, $i;
      $n -= $fibs[$i];
      ### sub: "$fibs[$i] to n=$n"
      --$i;
    }
  }
  ### @indices

  # X is low index, ie. how many low 0 bits in Zeckendorf form
  my $x = pop @indices;
  ### $x

  # Y is indices shifted down by $x and 2 more
  my $y = 0;
  my $shift = $x+2;
  foreach my $i (@indices) {
    ### y add: "ishift=".($i-$shift)." fib=".$fibs[$i-$shift]
    $y += $fibs[$i-$shift];
  }
  ### $shift
  ### $y

  return ($x+$self->{'x_start'},$y+$self->{'y_start'});
}

# phi = (sqrt(5)+1)/2
# (y+1)*phi = (y+1)*(sqrt(5)+1)/2
#           = ((y+1)*sqrt(5)+(y+1))/2
#           = (sqrt(5*(y+1)^2)+(y+1))/2
#
# from x=0,y=0
# N = floor((y+1)*Phi) * Fib(x+2) + y*Fib(x+1)
#
sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### WythoffArray xy_to_n(): "$x, $y"

  $x = round_nearest($x) - $self->{'x_start'};
  $y = round_nearest($y) - $self->{'y_start'};
  if ($x < 0 || $y < 0) {
    return undef;
  }

  my $zero = $x * 0 * $y;
  $x += 2;
  if (is_infinite($x)) { return $x; }
  if (is_infinite($y)) { return $y; }

  my @bits = bit_split_lowtohigh($x);
  ### @bits
  pop @bits; # discard high 1-bit

  my $yplus1 = $zero + $y+1;   # inherit bigint from $x perhaps

  # spectrum(Y+1) so Y,Ybefore are notional two values at X=-2 and X=-1
  my $ybefore = int((_sqrtint(5*$yplus1*$yplus1) + $yplus1) / 2);
  ### $ybefore

  # k=1, Fk1=F[k-1]=0, Fk=F[k]=1
  my $Fk1 = $zero;
  my $Fk  = 1 + $zero;

  my $add = -2;
  while (@bits) {
    ### remaining bits: @bits
    ### Fk1: "$Fk1"
    ### Fk: "$Fk"

    # two squares and some adds
    # F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
    # F[2k-1] =   F[k]^2 + F[k-1]^2
    # F[2k] = F[2k+1] - F[2k-1]
    #
    $Fk *= $Fk;
    $Fk1 *= $Fk1;
    my $F2kplus1 = 4*$Fk - $Fk1 + $add;
    $Fk1 += $Fk; # F[2k-1]
    my $F2k = $F2kplus1 - $Fk1;

    if (pop @bits) {  # high to low
      $Fk1 = $F2k;     # F[2k]
      $Fk = $F2kplus1; # F[2k+1]
      $add = -2;
    } else {
      # $Fk1 is F[2k-1] already
      $Fk = $F2k;  # F[2k]
      $add = 2;
    }
  }

  ### final pair ...
  ### Fk1: "$Fk1"
  ### Fk: "$Fk"
  ### @bits

  return ($Fk*$ybefore + $Fk1*$y);
}

# exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### WythoffArray rect_to_n_range(): "$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;

  if ($x2 < $self->{'x_start'} || $y2 < $self->{'y_start'}) {
    ### all outside first quadrant ...
    return (1, 0);
  }

  # bottom left into first quadrant
  $x1 = max($x1, $self->{'x_start'});
  $y1 = max($y1, $self->{'y_start'});

  return ($self->xy_to_n($x1,$y1),    # bottom left
          $self->xy_to_n($x2,$y2));   # top right
}

1;
__END__