1 /* mpz_fib_ui(result, n) -- Set RESULT to the Nth Fibonacci number.
3 Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 2.1 of the License, or (at your
10 option) any later version.
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20 MA 02111-1307, USA. */
25 /* This is fast, but could be made somewhat faster and neater.
26 The timing is somewhat fluctuating for even/odd sizes because
27 of the extra hair used to save variables and operations. Here
28 are a few things one might want to address:
29 1. Avoid using 4 intermediate variables in mpz_fib_bigcase.
30 2. Call mpn functions directly. Straightforward for these functions.
31 3. Merge the three functions into one.
34 Consider using the Lucas numbers L[n] as an auxiliary sequence, making
35 it possible to do the "doubling" operation in mpz_fib_bigcase with two
36 squares rather than two multiplies. The formulas are a little more
37 complicated, something like the following (untested).
39 F[2n] = ((F[n]+L[n])^2 - 6*F[n]^2 - 4*(-1)^n) / 2
40 L[2n] = 5*F[n]^2 + 2*(-1)^n
42 F[2n+1] = (F[2n] + L[2n]) / 2
43 L[2n+1] = (5*F[2n] + L[2n]) / 2
45 The Lucas number that comes for free here could even be returned.
47 Maybe there's formulas with two squares using just F[n], but I don't
51 /* Determine the needed storage for Fib(n). */
52 #define FIB_SIZE(n) (((mp_size_t) ((n)*0.695)) / BITS_PER_MP_LIMB + 2)
54 static void mpz_fib_bigcase _PROTO ((mpz_t, mpz_t, unsigned long int));
55 static void mpz_fib_basecase _PROTO ((mpz_t, mpz_t, unsigned long int));
59 #define FIB_THRESHOLD 60
64 mpz_fib_ui (mpz_t r, unsigned long int n)
77 if (n < FIB_THRESHOLD)
78 mpz_fib_basecase (t1, r, n);
80 mpz_fib_bigcase (t1, r, n);
87 mpz_fib_basecase (mpz_t t1, mpz_t t2, unsigned long int n)
89 mpz_fib_basecase (t1, t2, n)
95 unsigned long int m, i;
100 for (i = 0; i < m; i++)
102 mpz_add (t1, t1, t2);
103 mpz_add (t2, t1, t2);
107 mpz_sub (t1, t2, t1);
108 mpz_sub (t2, t2, t1); /* trick: recover t1 value just overwritten */
114 mpz_fib_bigcase (mpz_t t1, mpz_t t2, unsigned long int n)
116 mpz_fib_bigcase (t1, t2, n)
122 unsigned long int n2;
124 mpz_t x1, x2, u1, u2;
127 for (n2 = n; n2 >= FIB_THRESHOLD; n2 /= 2)
130 mpz_fib_basecase (t1, t2, n2);
137 for (i = ni - 1; i >= 0; i--)
139 mpz_mul_2exp (x1, t1, 1);
140 mpz_mul_2exp (x2, t2, 1);
142 mpz_add (x1, x1, t2);
143 mpz_sub (x2, x2, t1);
145 mpz_mul (u1, t2, x1);
146 mpz_mul (u2, t1, x2);
148 if (((n >> i) & 1) == 0)
150 mpz_sub (t1, u1, u2);
156 mpz_mul_2exp (t2, u1, 1);
157 mpz_sub (t2, t2, u2);