# Copyright 2010, 2011, 2012, 2013, 2014, 2016, 2019, 2020 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 <http://www.gnu.org/licenses/>.
package Math::NumSeq::Fibonacci;
use 5.004;
use strict;
use vars '$VERSION','@ISA';
$VERSION = 75;
use Math::NumSeq::Base::Sparse; # FIXME: implement pred() directly ...
@ISA = ('Math::NumSeq::Base::Sparse');
use Math::NumSeq;
*_is_infinite = \&Math::NumSeq::_is_infinite;
*_to_bigint = \&Math::NumSeq::_to_bigint;
# uncomment this to run the ### lines
# use Smart::Comments;
# use constant name => Math::NumSeq::__('Fibonacci Numbers');
use constant description => Math::NumSeq::__('The Fibonacci numbers 1,1,2,3,5,8,13,21, etc, each F(i) = F(i-1) + F(i-2), starting from 1,1.');
use constant values_min => 0;
use constant i_start => 0;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;
#------------------------------------------------------------------------------
# cf A105527 - index when n-nacci exceeds Fibonacci
# A020695 - Pisot 2,3,5,8,etc starting OFFSET=0
# A212804 - starting 1,0 OFFSET=0
use constant oeis_anum => 'A000045'; # Fibonacci starting at i=0 0,1,1,2,3
#------------------------------------------------------------------------------
# $uv_limit is the biggest Fibonacci number f0 for which both f0 and f1 fit
# into a UV. Upon reaching $uv_limit the next step will require BigInt.
# $uv_i_limit is the i index of $uv_limit.
#
my $uv_limit;
my $uv_i_limit = 0;
{
# Float integers too in 32 bits ?
# my $max = 1;
# for (1 .. 256) {
# my $try = $max*2 + 1;
# ### $try
# if ($try == 2*$max || $try == 2*$max+2) {
# last;
# }
# $max = $try;
# }
my $max = ~0;
# 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 $f0 = 1;
my $f1 = 1;
my $prev_f0;
while ($f0 <= $max - $f1) {
$prev_f0 = $f0;
($f1,$f0) = ($f1+$f0,$f1);
$uv_i_limit++;
}
### Fibonacci UV limit ...
### $prev_f0
### $f0
### $f1
### ~0 : ~0
$uv_limit = $prev_f0;
### $uv_limit
### $uv_i_limit
__PACKAGE__->ith($uv_i_limit) == $uv_limit
or warn "Oops, wrong uv_i_limit";
}
sub rewind {
my ($self) = @_;
### Fibonacci rewind()
$self->{'f0'} = 0;
$self->{'f1'} = 1;
$self->{'i'} = $self->i_start;
}
sub seek_to_i {
my ($self, $i) = @_;
($self->{'f0'}, $self->{'f1'}) = $self->ith_pair($i);
$self->{'i'} = $i;
}
sub next {
my ($self) = @_;
### Fibonacci next(): "f0=$self->{'f0'}, f1=$self->{'f1'}"
(my $ret,
$self->{'f0'},
$self->{'f1'})
= ($self->{'f0'},
$self->{'f1'},
$self->{'f0'} + $self->{'f1'});
### $ret
if ($ret == $uv_limit) {
### go to bigint f1 ...
$self->{'f1'} = _to_bigint($self->{'f1'});
}
return ($self->{'i'}++, $ret);
}
# F[k-1] + F[k] = F[k+1]
# F[k-1] = F[k+1] - F[k]
# F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k
# = (2F[k] + F[k+1] - F[k])*(2F[k] - (F[k+1] - F[k])) + 2*(-1)^k
# = (F[k] + F[k+1])*(2F[k] - F[k+1] + F[k]) + 2*(-1)^k
# = (F[k] + F[k+1])*(3F[k] - F[k+1]) + 2*(-1)^k
# F[2k] = F[k]*(F[k]+2F[k-1])
# = F[k]*(F[k]+2(F[k+1] - F[k]))
# = F[k]*(F[k]+2F[k+1] - 2F[k])
# = F[k]*(2F[k+1] - F[k])
sub ith {
my ($self, $i) = @_;
### Fibonacci ith(): $i
my $lowbit = ($i % 2);
my $pair_i = ($i - $lowbit) / 2;
my ($F0, $F1) = $self->ith_pair($pair_i);
if ($i > $uv_i_limit && ! ref $F0) {
### automatic BigInt as not another bignum class ...
$F0 = _to_bigint($F0);
$F1 = _to_bigint($F1);
}
# last step needing just one of F[2k] or F[2k+1] done by one multiply
# instead of two squares in the ith_pair() loop
#
if ($lowbit) {
$F0 = ($F0 + $F1) * (3*$F0 - $F1) + ($pair_i % 2 ? -2 : 2);
} else {
$F0 *= (2*$F1 - $F0);
}
return $F0;
}
sub ith_pair {
my ($self, $i) = @_;
### Fibonacci ith_pair(): $i
if (_is_infinite($i)) {
return ($i,$i);
}
my $neg = ($i < 0);
if ($neg) {
$i = -1-$i;
}
my @bits = _bit_split_hightolow($i+1);
### @bits
shift @bits; # drop high 1-bit
# k=1 which is the high bit of @bits
# $Fk1 = F[k-1] = 0
# $Fk = F[k] = 1
#
my $Fk1 = ($i * 0); # inherit bignum 0
if ($i >= $uv_i_limit && ! ref $Fk1) {
# automatic BigInt if not another number class
$Fk1 = _to_bigint(0);
}
my $Fk = $Fk1 + 1; # inherit bignum 1
my $add = -2; # (-1)^k
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 (shift @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;
}
}
if ($neg) {
($Fk1,$Fk) = ($Fk, $Fk1);
if ($i % 2) {
$Fk1 = -$Fk1;
} else {
$Fk = -$Fk;
}
}
### final ...
### Fk1: "$Fk1"
### Fk: "$Fk"
return ($Fk1, $Fk);
}
sub _bit_split_hightolow {
my ($n) = @_;
### _bit_split_hightolow(): "$n"
if (ref $n) {
if ($n->isa('Math::BigInt')
&& $n->can('as_bin')) {
### BigInt: $n->as_bin
return split //, substr($n->as_bin,2);
}
}
my @bits;
while ($n) {
push @bits, $n % 2;
$n = int($n/2);
}
return reverse @bits;
}
use constant 1.02 _PHI => (1 + sqrt(5)) / 2;
use constant 1.02 _BETA => -1/_PHI;
sub value_to_i_estimate {
my ($self, $value) = @_;
if (_is_infinite($value)) {
return $value;
}
if ($value <= 0) {
return 0;
}
if (defined (my $blog2 = _blog2_estimate($value))) {
# i ~= (log2(F(i)) + log2(phi)) / log2(phi-beta)
# with log2(x) = log(x)/log(2)
return int( ($blog2 + (log(_PHI - _BETA)/log(2)))
/ (log(_PHI)/log(2)) );
}
# i ~= (log(F(i)) + log(phi)) / log(phi-beta)
return int( (log($value) + log(_PHI - _BETA))
/ log(_PHI) );
}
sub _UNTESTED__value_to_i {
my ($self, $value) = @_;
if ($value < 0) { return undef; }
my $i = $self->value_to_i_estimate($value) - 3;
if (_is_infinite($i)) { return $i; }
if ($i < 0) { $i = 0; }
my ($f0,$f1) = $self->ith_pair($i);
foreach (1 .. 10) {
if ($f0 == $value) {
return $i;
}
last if $f0 > $value;
if ($i == $uv_i_limit && ! ref $f0) {
# automatic BigInt if not another number class
$f1 = _to_bigint($f1);
}
($f0, $f1) = ($f1, $f0+$f1);
$i += 1;
}
return undef;
}
#------------------------------------------------------------------------------
# generic, shared
# if $n is a BigInt, BigRat or BigFloat then return an estimate of log base 2
# otherwise return undef.
#
# For Math::BigInt
#
# For BigRat the calculation is just a bit count of the numerator less the
# denominator so may be off by +/-1 or +/-2 or some such. For
#
sub _blog2_estimate {
my ($n) = @_;
if (ref $n) {
### _blog2_estimate(): "$n"
if ($n->isa('Math::BigRat')) {
return ($n->numerator->copy->blog(2) - $n->denominator->copy->blog(2))->numify;
}
if ($n->isa('Math::BigFloat')) {
return $n->as_int->blog(2)->numify;
}
if ($n->isa('Math::BigInt')) {
return $n->copy->blog(2)->numify;
}
}
return undef;
}
1;
__END__
=for stopwords Ryde Math-NumSeq Ith bignum
=head1 NAME
Math::NumSeq::Fibonacci -- Fibonacci numbers
=head1 SYNOPSIS
use Math::NumSeq::Fibonacci;
my $seq = Math::NumSeq::Fibonacci->new;
my ($i, $value) = $seq->next;
=head1 DESCRIPTION
The Fibonacci numbers F(i) = F(i-1) + F(i-2) starting from 0,1,
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, ...
starting i=0
=head1 FUNCTIONS
See L<Math::NumSeq/FUNCTIONS> for behaviour common to all sequence classes.
=over 4
=item C<$seq = Math::NumSeq::Fibonacci-E<gt>new ()>
Create and return a new sequence object.
=back
=head2 Iterating
=over
=item C<($i, $value) = $seq-E<gt>next()>
Return the next index and value in the sequence.
When C<$value> exceeds the range of a Perl unsigned integer the return is a
C<Math::BigInt> to preserve precision.
=item C<$seq-E<gt>seek_to_i($i)>
Move the current sequence position to C<$i>. The next call to C<next()>
will return C<$i> and corresponding value.
=back
=head2 Random Access
=over
=item C<$value = $seq-E<gt>ith($i)>
Return the C<$i>'th Fibonacci number.
For negative <$i> the sequence is extended backwards as F[i]=F[i+2]-F[i+1].
The effect is the same Fibonacci numbers but negative at negative even i.
i F[i]
--- ----
0 0
-1 1
-2 -1 <----+ negative at even i
-3 2 |
-4 -3 <----+
When C<$value> exceeds the range of a Perl unsigned integer the return is a
C<Math::BigInt> to preserve precision.
=item C<$bool = $seq-E<gt>pred($value)>
Return true if C<$value> occurs in the sequence, so is a positive Fibonacci
number.
=item C<$i = $seq-E<gt>value_to_i_estimate($value)>
Return an estimate of the i corresponding to C<$value>. See L</Value to i
Estimate> below.
=back
=head1 FORMULAS
=head2 Ith
Fibonacci F[i] can be calculated by a powering procedure with two squares
per step. A pair of values F[k] and F[k-1] are maintained and advanced
according to bits of i from high to low
start k=1, F[k]=1, F[k-1]=0
add = -2 # 2*(-1)^k
loop
F[2k+1] = 4*F[k]^2 - F[k-1]^2 + add
F[2k-1] = F[k]^2 + F[k-1]^2
F[2k] = F[2k+1] - F[2k-1]
bit = next bit of i, high to low, skip high 1 bit
if bit == 1
take F[2k+1], F[2k] as new F[k],F[k-1]
add = -2 (for next loop)
else bit == 0
take F[2k], F[2k-1] as new F[k],F[k-1]
add = 2 (for next loop)
For the last (least significant) bit of i an optimization can be made with a
single multiple for that last step, instead of two squares.
bit = least significant bit of i
if bit == 1
F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + add
else
F[2k] = F[k]*(F[k]+2F[k-1])
The "add" amount is 2*(-1)^k which means +2 or -2 according to k odd or
even, which means whether the previous bit taken from i was 1 or 0. That
can be easily noted from each bit, to be used in the following loop
iteration or the final step F[2k+1] formula.
For small i it's usually faster to just successively add F[k+1]=F[k]+F[k-1],
but when in bignums the doubling k-E<gt>2k by two squares is faster than
doing k many individual additions for the same thing.
=head2 Value to i Estimate
F[i] increases as a power of phi, the golden ratio. The exact value is
F[i] = (phi^i - beta^i) / (phi - beta) # exactly
phi = (1+sqrt(5))/2 = 1.618
beta = -1/phi = -0.618
Since abs(beta)E<lt>1 the beta^i term quickly becomes small. So taking a
log (natural logarithm) to get i,
log(F[i]) ~= i*log(phi) - log(phi-beta)
i ~= (log(F[i]) + log(phi-beta)) / log(phi)
Or the same using log base 2 which can be estimated from the highest bit
position of a bignum,
log2(F[i]) ~= i*log2(phi) - log2(phi-beta)
i ~= (log2(F[i]) + log2(phi-beta)) / log2(phi)
=head1 SEE ALSO
L<Math::NumSeq>,
L<Math::NumSeq::LucasNumbers>,
L<Math::NumSeq::Fibbinary>,
L<Math::NumSeq::FibonacciWord>,
L<Math::NumSeq::Pell>,
L<Math::NumSeq::Tribonacci>
L<Math::Fibonacci>,
L<Math::Fibonacci::Phi>
=head1 HOME PAGE
L<http://user42.tuxfamily.org/math-numseq/index.html>
=head1 LICENSE
Copyright 2010, 2011, 2012, 2013, 2014, 2016, 2019, 2020 Kevin Ryde
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 <http://www.gnu.org/licenses/>.
=cut