File Coverage

blib/lib/Math/NumSeq/Fibbinary.pm

Criterion	Covered	Total	%
statement	131	135	97.0
branch	28	32	87.5
condition	6	8	75.0
subroutine	25	25	100.0
pod	11	11	100.0
total	201	211	95.2

# Copyright 2011, 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-NumSeq.
#
# Math-NumSeq 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-NumSeq 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-NumSeq.  If not, see .


# ZOrderCurve, ImaginaryBase tree shape
# DragonCurve repeating runs
#
# cf fxtbook ch38 p756
#
# cf visualizing
# http://cs-people.bu.edu/ilir/zecko/


package Math::NumSeq::Fibbinary;
use 5.004;
use strict;
use Carp;

use vars '$VERSION', '@ISA';
$VERSION = 71;
use Math::NumSeq;
@ISA = ('Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;
*_to_bigint = \&Math::NumSeq::_to_bigint;

use Math::NumSeq::Fibonacci;
*_bit_split_hightolow = \&Math::NumSeq::Fibonacci::_bit_split_hightolow;
*_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;

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


# use constant name => Math::NumSeq::__('Fibbinary Numbers');
use constant description => Math::NumSeq::__('Fibbinary numbers 0,1,2,4,5,8,9,etc, integers without adjacent 1-bits.');
use constant default_i_start => 0;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;

sub values_min {
  my ($self) = @_;
  return $self->ith($self->i_start);
}

#------------------------------------------------------------------------------
# cf A000119 - number of fibonacci sums forms
#    A003622 - n with odd Zeckendorf,  cf golden seq
#    A037011 - baum-sweet cubic, might be 1 iff i is in the fibbinary seq
#    A014417 - n in fibonacci base, the fibbinaries written out in binary
#    A139764 - smallest Zeckendorf term
#    A054204 - using only even Fibs
#
use constant oeis_anum => 'A003714';  # Fibbinary, OFFSET=0 start value=0

#------------------------------------------------------------------------------
# $self->{'i'}, $self->{'value'} are the next $i,$value to return.
# next() increments 'i' and steps 'value'.
# So the next value is calculated ahead of its actually being needed,
# but doing so

sub rewind {
  my ($self) = @_;
  $self->{'i'} = $self->i_start;
  $self->{'value'} = 0;
}
sub seek_to_i {
  my ($self, $i) = @_;
  if ($i < 0) {
    croak "Cannot seek to ",$i,", sequence begins at i=0";
  }
  $self->{'i'} = $i;
  $self->{'value'} = $self->ith($i);
}
sub seek_to_value {
  my ($self, $value) = @_;
  $self->seek_to_i($self->value_to_i_ceil($value));
}

sub next {
  my ($self) = @_;
  ### Fibbinary next() ...

  my $v = $self->{'value'};
  $self->{'value'} = _value_next($self,$v);
  return ($self->{'i'}++, $v);
}

sub _value_next {
  my ($self, $value) = @_;
  my $filled = ($value >> 1) | $value;
  my $mask = (($filled+1) ^ $filled) >> 1;

  ### value : sprintf('0b %6b',$value)
  ### filled: sprintf('0b %6b',$filled)
  ### mask  : sprintf('0b %6b',$mask)
  ### bit   : sprintf('0b %6b',$mask+1)
  ### newv  : sprintf('0b %6b',($value | $mask))

  return ($value | $mask) + 1;
}

sub ith {
  my ($self, $i) = @_;
  ### Fibbinary ith(): $i

  if (_is_infinite($i)) {
    return $i;
  }

  # 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 = ($i * 0);  # inherit bignum 0
    my $f1 = $f0 + 1;   # inherit bignum 1
    @fibs = ($f0);
    while ($f0 <= $i-$f1) {
      ($f1,$f0) = ($f1+$f0,$f1);
      push @fibs, $f1;
    }
  }
  ### @fibs

  my $value = 0;
  while (my $f = pop @fibs) {
    ### at: "$f i=$i  value=$value"
    $value *= 2;
    if ($i >= $f) {
      $value += 1;
      $i -= $f;
      ### sub: "$f to i=$i value=$value"

      # never consecutive fibs, so pop without comparing to i
      pop @fibs || last;
      $value *= 2;
    }
  }
  return $value;
}

sub pred {
  my ($self, $value) = @_;
  ### Fibbinary pred(): $value

  my $int;
  unless ($value >= 0
          && $value == ($int = int($value))) {
    return 0;
  }

  # go to BigInt if NV floating point integer bigger than UV, since "&"
  # operator will cast to a UV and lose bits
  if ($int > ~0 && ! ref $int) {
    $int = _to_bigint(sprintf('%.0f',$int));
    ### use BigInt: $int
    ### str: sprintf('%.0f',$int)
  }

  ### and: ($int & ($int >> 1)).''
  return ! ($int & ($int >> 1));
}

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

sub value_to_i_floor {
  my ($self, $value) = @_;
  ### Fibbinary value_to_i_floor(): $value
  if ($value < 0) { return 0; }
  my ($i) = _value_to_i_and_floor($value);
  return $i;
}
sub value_to_i_ceil {
  my ($self, $value) = @_;
  ### Fibbinary value_to_i_ceil(): $value
  if ($value < 0) { return 0; }
  my ($i,$floor) = _value_to_i_and_floor($value);
  return $i + $floor;
}
sub value_to_i {
  my ($self, $value) = @_;
  ### Fibbinary value_to_i(): $value
  if ($value < 0) { return undef; }
  my ($i,$floor) = _value_to_i_and_floor($value);
  return ($floor ? undef : $i);
}

# return ($i, $floor)
sub _value_to_i_and_floor {
  my ($value) = @_;

  if (_is_infinite($value)) {
    return ($value,
            0); # reckon infinite as not rounded
  }

  my $floor;
  {
    my $int = int($value);
    $floor = ($value == $int ? 0 : 1);
    $value = $int
      || return (0, $floor);  # i=0 not handled below
  }

  my @bits = _bit_split_hightolow($value);
  my @fibs;
  {
    my $f0 = ($value * 0);  # inherit bignum 0
    my $f1 = $f0 + 1;       # inherit bignum 1
    foreach (@bits) {
      ($f1,$f0) = ($f1+$f0,$f1);
      push @fibs, $f1;
    }
  }
  ### @fibs

  my $prev_bit = shift @bits; # high 1-bit
  my $i = pop @fibs;

  ### initial i: $i

  foreach my $bit (@bits) {  # high to low
    my $fib = pop @fibs;
    ### $bit
    ### $fib

    if ($bit) {
      if ($prev_bit) {
        ### consecutive bits 11xxx, round down to 10xxx with xxx=1010 ...
        while (@fibs) {
          $i += pop @fibs;
          pop @fibs;
        }
        return ($i,
                1);  # rounded down
      }
      $i += $fib;
      ### add i to: $i
    }
    $prev_bit = $bit;
  }
  ### exact i: "$i"
  return ($i,
          $floor);  # not rounded, unless $value was fractional
}

#------------------------------------------------------------------------------
# value_to_i_estimate()

use constant 1.02 _PHI => (1 + sqrt(5)) / 2;

# (phi-beta) = phi+1/phi = 2phi-1
#
# value=2^k
# log(value) = k*log(2)
# k = log(value)/log(2)
# i = F(k)
#   = phi^(k+1) / (phi-beta)
#   = phi^k * C where C=phi/(phi-beta) ~= 0.72
# log(i/C) = k*log(phi)
# k = log(i/C)/log(phi)
#
# log(i/C)/log(phi) = log(value)/log(2)
# log(i/C) = log(value) * log(phi)/log(2)
# i/C = e^ (log(value) * log(phi)/log(2))
# i/C = (e^log(value)) ^ (log(phi)/log(2)))
# i = C * value ^ (log(phi)/log(2)))
#
# log(phi)/log(2) ~= 0.694
#

sub value_to_i_estimate {
  my ($self, $value) = @_;

  if ($value <= 0) {
    return 0;
  }

  $value = int($value);
  if (my $blog2 = Math::NumSeq::Fibonacci::_blog2_estimate($value)) {
    my $shift = int ((1 - log(_PHI)/log(2))
                     * Math::NumSeq::Fibonacci::_blog2_estimate($value));
    return $value >> $shift;
  }

  return int ((((_PHI + 1/_PHI)/_PHI))
              * $value ** (log(_PHI)/log(2)));
}

# Can get close taking bits low to high and tweaking for consecutive 1s.
# But the high to low of the full value_to_i_floor() is only a little extra
# work.
#
# sub value_to_i_estimate {
#   my ($self, $value) = @_;
#   ### Fibbinary value_to_i_estimate(): $value
#
#   if (_is_infinite($value)) {
#     return $value;
#   }
#
#   my $f0 = my $f1 = ($value * 0)+1;  # inherit bignum 1
#   my $i = 0;
#
#   my $prev_bit = 0;
#   while ($value) {
#     my $bit = $value % 2;
#     if ($bit) {
#       if ($prev_bit) {
#         $i += $f0;
#       } else {
#         $i += $f1;
#       }
#     }
#     $prev_bit = $bit;
#     ($f1,$f0) = ($f1+$f0,$f1);
#     $value = int($value/2);
#   }
#   return $i;
# }

1;
__END__