``````# 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
# 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 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]
} else {
\$Fk = \$F2k;  # F[2k]
}
}

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

Math::NumSeq::Fibonacci -- Fibonacci numbers

use Math::NumSeq::Fibonacci;
my \$seq = Math::NumSeq::Fibonacci->new;
my (\$i, \$value) = \$seq->next;

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

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

=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

=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

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

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

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)

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>

L<http://user42.tuxfamily.org/math-numseq/index.html>

Copyright 2010, 2011, 2012, 2013, 2014, 2016, 2019, 2020 Kevin Ryde

Math-NumSeq is free software; you can redistribute it and/or modify it