[seqfan] Re: A1595 mentioned in Medium article
M. F. Hasler
oeis at hasler.fr
Mon Nov 9 23:47:45 CET 2020
Marc,
It certainly does not harm to recall how matrix and other symbolic powers
(e.g. also for "a+b sqrt(5)" in our case) are computed,
through binary exponentiation in log2(n) operations!
And since we're at it, let's mention that thi is also one of the
best methods to compute an arbitrary term in any linear recurrence,
not only for Fibonacci or generalized Fibonacci but also for higher order:
in that case the matrix M will have the signature in its last row, 1's just
above the diagonal, and we'd multiply with the vector of initial values
rather than pick an element of the matrix M^n as we can do here because the
initial values are (0,1).
At the heart is still the computation of M^n in O(log2(n)) operations.
PS: I write to also fix a possibly confusing typo at the same occasion: my
phone changed the intended "naïve" to "native", in the first line.
- Maximilian
On Mon, Nov 9, 2020 at 6:00 PM Marc LeBrun <mlb at well.com> wrote:
> Maximilian, all great points, thanks for the additional comments!
>
> Perhaps I shouldn't have included that pseudocode; I was just following
> the pattern of the original email and thought it might be helpful.
>
> Be well and prosper! --MLB
>
>
> > On Nov 9, 2020, at 1:50 PM, M. F. Hasler <oeis at hasler.fr> wrote:
> >
> > Marc,
> >
> > I think Alonso is aware of that, his article gives the "native" method
> > explicitly as bad example to produce a timeout.
> > The matrix computation method, Fib(n) = ( M^n )[1,2] = ( M^n )[2,1] = (
> > M^(n-1) )[2,2] with M = [0,1; 1,1]
> > is already given in some of the programs of oeis.org/A45, e.g., Julia.
> >
> > But two more remarks should be made, given you mentioned this method:
> >
> > - first, there is an explicit formula,
> > Fibonacci(n) = (Phi^n - Phi^-n)/sqrt(5) = floor(Phi^n / sqrt(5))
> > with the golden ratio Phi = (sqrt(5)+1)/2
> > which is usable "out of the box" with standard (double) precision [38
> > digits] up to almost Fib(200) (~ 10^41)
> > or beyond either with increased floating point precision,
> > or using exact arithmetic in (some implementation of) Z[sqrt(5)], see,
> > e.g., the PARI code with "quadgen(5)".
> >
> > - The latter would probably, like your code for matrix multiplication,
> use
> > what is called binary exponentiation.
> > It's a bit confusing to make a procedure called matfibpow(), because
> > this is just the standard way of computing integer powers of matrices or
> > any other objects.
> > All sensible programming languages do implement M^n in that way.
> > I would find it very strange if Java didn't.
> > Let all SeqFans be assured that PARI/GP does that when you type
> > [0,1;1,1]^n, without need for a custom procedure.
> >
> > Regards,
> > Maximilian
> >
> > On Mon, 9 Nov 2020, 16:47 Marc LeBrun, <mlb at well.com> wrote:
> >
> > It's sometimes useful to know that -- if you just want to get the N-th
> > Fibonacci number without computing so many prior ones -- you can do much
> > better: there's an approach where you can compute it with fixed storage
> in
> > O(log2 N) operations, by computing the N-th power of a 2x2 matrix the
> same
> > way numeric N-th powers are efficiently calculated:
> >
> > To do this, let the identity matrix be
> >
> > 1 0
> > I :=
> > 0 1
> >
> > and also define a "multiplier" matrix
> >
> > 0 1
> > M :=
> > 1 1
> >
> > and then evaluate the following recursive procedure
> >
> > matfibpow(N):=
> > if N==0
> > then return I
> > else
> > let F = matfibpow(floor(N/2))
> > set F = F * F
> > if N is odd
> > then return M * F
> > else return F
> >
> > Finally, the N-th Fibonacci number will be the [2,1] element of the
> > returned matrix.
> >
> > Basically M encodes the the linear recurrence relationship. This
> approach
> > can be applied for different initial values by modifying the final step,
> > and for different recurrences by using a different multiplier matrix M.
> >
> >
> >> On Nov 9, 2020, at 10:22 AM, Alonso Del Arte <alonso.delarte at gmail.com>
> > wrote:
> >>
> >> Thank you for finding A1610, David. I think I made a mistake when I
> tried
> >> to calculate that sequence, so it didn't matter if I omitted the initial
> >> terms. I decided I did not actually need the sequence for the article,
> so
> > I
> >> set it aside. Now that I've had time to sleep on it, I've verified that
> >> A1610 with a 1 joined at the beginning is the sequence I was looking
> for.
> >>
> >> I should mention that I omitted one line from the article:
> >>
> >> public static BigInteger fibonacci(int n) {
> >> *counter++;*
> >> switch (n) {
> >> case 0: return BigInteger.ZERO;
> >> case 1: return BigInteger.ONE;
> >> default:
> >> return fibonacci(n - 2).add(fibonacci(n - 1));
> >> }
> >> }
> >>
> >> counter is a static variable that keeps track of how many times
> > fibonacci()
> >> is called.
> >>
> >> Fred wrote that
> >>
> >>> There is nothing wrong with using the recursion directly for small
> >> values: simply compute each F(m) in turn for 2 <= m <= n ,
> >> requiring order n additions.
> >>
> >> I think that would be the case if the program remembers previously
> > computed
> >> values.
> >>
> >> I remember a JavaScript whiteboarding meet-up where the moderator used a
> >> JavaScript implementation like the Java one quoted above. She used it to
> >> solve a Project Euler problem regarding even Fibonacci numbers. It was
> > only
> >> necessary to go up to Fibonacci(50), if I recall correctly. It took
> >> something like a full minute to give a result. No one said anything.
> >>
> >> The simple act of saving previously computed values cuts down
> computation
> >> times: Fibonacci(50) in a millisecond or two instead of a full minute,
> and
> >> Fibonacci(100) in a few milliseconds instead of the theoretical 36
> million
> >> years.
> >>
> >> Al
> >>
> >> On Mon, Nov 9, 2020 at 8:47 AM David Seal <david.j.seal at gwynmop.com>
> > wrote:
> >>
> >>>> On 09/11/2020 05:34 Alonso Del Arte <alonso.delarte at gmail.com> wrote:
> >>>> ...
> >>>> I've been wondering about Fibonacci(*n*) − A001595(*n*). That's
> > probably
> >>>> already in the OEIS, though maybe without signs. It does change sign,
> >>> right?
> >>>
> >>> From the definition of the Fibonacci numbers (A000045), which is "F(n)
> =
> >>> F(n-1) + F(n-2) with F(0) = 0 and F(1) = 1", and from the definition of
> >>> A001595, which is "a(n) = a(n-1) + a(n-2) + 1, with a(0) = a(1) = 1",
> it
> > is
> >>> very easy to prove by induction that a(n) >= F(n) for all n >= 0 (with
> >>> equality if and only if n = 1), and so the difference of the two
> > sequences
> >>> does not change sign.
> >>>
> >>> Is there some subtlety about this question whose significance I'm
> > missing?
> >>> - for instance, a meaning of "*n*"?
> >>>
> >>> A search for the first ten values of A001595(n) - A000045(n), which are
> >>> 1,0,2,3,6,10,17,28,46,75, says that it doesn't match any sequence in
> the
> >>> OEIS. However, leaving off its initial 1 finds A001610, described as
> > "a(n)
> >>> = a(n-1) + a(n-2) + 1" - an incomplete definition, so I have submitted
> a
> >>> change to add "with a(0) = 0 and a(1) = 2". With that completed
> > definition,
> >>> it's also very easy to prove by induction that A001610(n) =
> A001595(n+1)
> > -
> >>> A000045(n+1) for all n >= 0.
> >>>
> >>> David
> >>>
> >>> --
> >>> Seqfan Mailing list - http://list.seqfan.eu/
> >>>
> >>
> >>
> >> --
> >> Alonso del Arte
> >> Author at SmashWords.com
> >> <https://www.smashwords.com/profile/view/AlonsoDelarte>
> >> Musician at ReverbNation.com <http://www.reverbnation.com/alonsodelarte
> >
> >>
> >> --
> >> Seqfan Mailing list - http://list.seqfan.eu/
> >
> >
> > --
> > Seqfan Mailing list - http://list.seqfan.eu/
> >
> >
> > On Mon, 9 Nov 2020, 16:47 Marc LeBrun, <mlb at well.com> wrote:
> >
> >> It's sometimes useful to know that -- if you just want to get the N-th
> >> Fibonacci number without computing so many prior ones -- you can do much
> >> better: there's an approach where you can compute it with fixed storage
> in
> >> O(log2 N) operations, by computing the N-th power of a 2x2 matrix the
> same
> >> way numeric N-th powers are efficiently calculated:
> >>
> >> To do this, let the identity matrix be
> >>
> >> 1 0
> >> I :=
> >> 0 1
> >>
> >> and also define a "multiplier" matrix
> >>
> >> 0 1
> >> M :=
> >> 1 1
> >>
> >> and then evaluate the following recursive procedure
> >>
> >> matfibpow(N):=
> >> if N==0
> >> then return I
> >> else
> >> let F = matfibpow(floor(N/2))
> >> set F = F * F
> >> if N is odd
> >> then return M * F
> >> else return F
> >>
> >> Finally, the N-th Fibonacci number will be the [2,1] element of the
> >> returned matrix.
> >>
> >> Basically M encodes the the linear recurrence relationship. This
> approach
> >> can be applied for different initial values by modifying the final step,
> >> and for different recurrences by using a different multiplier matrix M.
> >>
> >>
> >>> On Nov 9, 2020, at 10:22 AM, Alonso Del Arte <alonso.delarte at gmail.com
> >
> >> wrote:
> >>>
> >>> Thank you for finding A1610, David. I think I made a mistake when I
> tried
> >>> to calculate that sequence, so it didn't matter if I omitted the
> initial
> >>> terms. I decided I did not actually need the sequence for the article,
> >> so I
> >>> set it aside. Now that I've had time to sleep on it, I've verified that
> >>> A1610 with a 1 joined at the beginning is the sequence I was looking
> for.
> >>>
> >>> I should mention that I omitted one line from the article:
> >>>
> >>> public static BigInteger fibonacci(int n) {
> >>> *counter++;*
> >>> switch (n) {
> >>> case 0: return BigInteger.ZERO;
> >>> case 1: return BigInteger.ONE;
> >>> default:
> >>> return fibonacci(n - 2).add(fibonacci(n - 1));
> >>> }
> >>> }
> >>>
> >>> counter is a static variable that keeps track of how many times
> >> fibonacci()
> >>> is called.
> >>>
> >>> Fred wrote that
> >>>
> >>>> There is nothing wrong with using the recursion directly for small
> >>> values: simply compute each F(m) in turn for 2 <= m <= n ,
> >>> requiring order n additions.
> >>>
> >>> I think that would be the case if the program remembers previously
> >> computed
> >>> values.
> >>>
> >>> I remember a JavaScript whiteboarding meet-up where the moderator used
> a
> >>> JavaScript implementation like the Java one quoted above. She used it
> to
> >>> solve a Project Euler problem regarding even Fibonacci numbers. It was
> >> only
> >>> necessary to go up to Fibonacci(50), if I recall correctly. It took
> >>> something like a full minute to give a result. No one said anything.
> >>>
> >>> The simple act of saving previously computed values cuts down
> computation
> >>> times: Fibonacci(50) in a millisecond or two instead of a full minute,
> >> and
> >>> Fibonacci(100) in a few milliseconds instead of the theoretical 36
> >> million
> >>> years.
> >>>
> >>> Al
> >>>
> >>> On Mon, Nov 9, 2020 at 8:47 AM David Seal <david.j.seal at gwynmop.com>
> >> wrote:
> >>>
> >>>>> On 09/11/2020 05:34 Alonso Del Arte <alonso.delarte at gmail.com>
> wrote:
> >>>>> ...
> >>>>> I've been wondering about Fibonacci(*n*) − A001595(*n*). That's
> >> probably
> >>>>> already in the OEIS, though maybe without signs. It does change sign,
> >>>> right?
> >>>>
> >>>> From the definition of the Fibonacci numbers (A000045), which is
> "F(n) =
> >>>> F(n-1) + F(n-2) with F(0) = 0 and F(1) = 1", and from the definition
> of
> >>>> A001595, which is "a(n) = a(n-1) + a(n-2) + 1, with a(0) = a(1) = 1",
> >> it is
> >>>> very easy to prove by induction that a(n) >= F(n) for all n >= 0 (with
> >>>> equality if and only if n = 1), and so the difference of the two
> >> sequences
> >>>> does not change sign.
> >>>>
> >>>> Is there some subtlety about this question whose significance I'm
> >> missing?
> >>>> - for instance, a meaning of "*n*"?
> >>>>
> >>>> A search for the first ten values of A001595(n) - A000045(n), which
> are
> >>>> 1,0,2,3,6,10,17,28,46,75, says that it doesn't match any sequence in
> the
> >>>> OEIS. However, leaving off its initial 1 finds A001610, described as
> >> "a(n)
> >>>> = a(n-1) + a(n-2) + 1" - an incomplete definition, so I have
> submitted a
> >>>> change to add "with a(0) = 0 and a(1) = 2". With that completed
> >> definition,
> >>>> it's also very easy to prove by induction that A001610(n) =
> >> A001595(n+1) -
> >>>> A000045(n+1) for all n >= 0.
> >>>>
> >>>> David
> >>>>
> >>>> --
> >>>> Seqfan Mailing list - http://list.seqfan.eu/
> >>>>
> >>>
> >>>
> >>> --
> >>> Alonso del Arte
> >>> Author at SmashWords.com
> >>> <https://www.smashwords.com/profile/view/AlonsoDelarte>
> >>> Musician at ReverbNation.com <
> http://www.reverbnation.com/alonsodelarte>
> >>>
> >>> --
> >>> Seqfan Mailing list - http://list.seqfan.eu/
> >>
> >>
> >> --
> >> Seqfan Mailing list - http://list.seqfan.eu/
> >>
> >
> > --
> > Seqfan Mailing list - http://list.seqfan.eu/
>
>
> --
> Seqfan Mailing list - http://list.seqfan.eu/
>
More information about the SeqFan
mailing list