version 1.1.1.1, 2000/09/09 14:12:51 |
version 1.1.1.2, 2003/08/25 16:06:33 |
|
|
/* mpz_fib_ui(result, n) -- Set RESULT to the Nth Fibonacci number. |
/* mpz_fib_ui -- calculate Fibonacci numbers. |
|
|
Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc. |
Copyright 2000, 2001 Free Software Foundation, Inc. |
|
|
This file is part of the GNU MP Library. |
This file is part of the GNU MP Library. |
|
|
Line 19 along with the GNU MP Library; see the file COPYING.LI |
|
Line 19 along with the GNU MP Library; see the file COPYING.LI |
|
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
MA 02111-1307, USA. */ |
MA 02111-1307, USA. */ |
|
|
|
#include <stdio.h> |
#include "gmp.h" |
#include "gmp.h" |
#include "gmp-impl.h" |
#include "gmp-impl.h" |
|
#include "longlong.h" |
|
|
/* This is fast, but could be made somewhat faster and neater. |
|
The timing is somewhat fluctuating for even/odd sizes because |
|
of the extra hair used to save variables and operations. Here |
|
are a few things one might want to address: |
|
1. Avoid using 4 intermediate variables in mpz_fib_bigcase. |
|
2. Call mpn functions directly. Straightforward for these functions. |
|
3. Merge the three functions into one. |
|
|
|
Said by Kevin: |
/* change to "#define TRACE(x) x" to get some traces */ |
Consider using the Lucas numbers L[n] as an auxiliary sequence, making |
#define TRACE(x) |
it possible to do the "doubling" operation in mpz_fib_bigcase with two |
|
squares rather than two multiplies. The formulas are a little more |
|
complicated, something like the following (untested). |
|
|
|
F[2n] = ((F[n]+L[n])^2 - 6*F[n]^2 - 4*(-1)^n) / 2 |
|
L[2n] = 5*F[n]^2 + 2*(-1)^n |
|
|
|
F[2n+1] = (F[2n] + L[2n]) / 2 |
/* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low |
L[2n+1] = (5*F[2n] + L[2n]) / 2 |
limb because the result F[2k+1] is an F[4m+3] and such numbers are always |
|
== 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7. (This is |
|
the same as in mpn_fib2_ui.) |
|
|
The Lucas number that comes for free here could even be returned. |
In the F[2k+1] for k even, the +2 won't give a carry out of the low limb |
|
in normal circumstances. This is an F[4m+1] and we claim that F[3*2^b+1] |
|
== 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence |
|
if n < 2^BITS_PER_MP_LIMB then F[n] cannot have a low limb of 0 or 1. No |
|
proof for this claim, but it's been verified up to b==32 and has such a |
|
nice pattern it must be true :-). Of interest is that F[3*2^b] == 0 mod |
|
2^(b+1) seems to hold too. |
|
|
Maybe there's formulas with two squares using just F[n], but I don't |
When n >= 2^BITS_PER_MP_LIMB, which can arise in test setups with a small |
know of any. |
limb, then the low limb of F[4m+1] can certainly be 1, and an mpn_add_1 |
*/ |
must be used. */ |
|
|
/* Determine the needed storage for Fib(n). */ |
|
#define FIB_SIZE(n) (((mp_size_t) ((n)*0.695)) / BITS_PER_MP_LIMB + 2) |
|
|
|
static void mpz_fib_bigcase _PROTO ((mpz_t, mpz_t, unsigned long int)); |
|
static void mpz_fib_basecase _PROTO ((mpz_t, mpz_t, unsigned long int)); |
|
|
|
|
|
#ifndef FIB_THRESHOLD |
|
#define FIB_THRESHOLD 60 |
|
#endif |
|
|
|
void |
void |
#if __STDC__ |
mpz_fib_ui (mpz_ptr fn, unsigned long n) |
mpz_fib_ui (mpz_t r, unsigned long int n) |
|
#else |
|
mpz_fib_ui (r, n) |
|
mpz_t r; |
|
unsigned long int n; |
|
#endif |
|
{ |
{ |
if (n == 0) |
mp_ptr fp, xp, yp; |
mpz_set_ui (r, 0); |
mp_size_t size, xalloc; |
else |
unsigned long n2; |
|
mp_limb_t c, c2; |
|
TMP_DECL (marker); |
|
|
|
if (n <= FIB_TABLE_LIMIT) |
{ |
{ |
mpz_t t1; |
PTR(fn)[0] = FIB_TABLE (n); |
mpz_init (t1); |
SIZ(fn) = (n != 0); /* F[0]==0, others are !=0 */ |
if (n < FIB_THRESHOLD) |
return; |
mpz_fib_basecase (t1, r, n); |
|
else |
|
mpz_fib_bigcase (t1, r, n); |
|
mpz_clear (t1); |
|
} |
} |
} |
|
|
|
static void |
n2 = n/2; |
#if __STDC__ |
xalloc = MPN_FIB2_SIZE (n2) + 1; |
mpz_fib_basecase (mpz_t t1, mpz_t t2, unsigned long int n) |
MPZ_REALLOC (fn, 2*xalloc+1); |
#else |
fp = PTR (fn); |
mpz_fib_basecase (t1, t2, n) |
|
mpz_t t1; |
|
mpz_t t2; |
|
unsigned long int n; |
|
#endif |
|
{ |
|
unsigned long int m, i; |
|
|
|
mpz_set_ui (t1, 0); |
TMP_MARK (marker); |
mpz_set_ui (t2, 1); |
TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc); |
m = n/2; |
size = mpn_fib2_ui (xp, yp, n2); |
for (i = 0; i < m; i++) |
|
|
TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n", |
|
n >> 1, size, n&1); |
|
mpn_trace ("xp", xp, size); |
|
mpn_trace ("yp", yp, size)); |
|
|
|
if (n & 1) |
{ |
{ |
mpz_add (t1, t1, t2); |
/* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k */ |
mpz_add (t2, t1, t2); |
mp_size_t xsize, ysize; |
} |
|
if ((n & 1) == 0) |
|
{ |
|
mpz_sub (t1, t2, t1); |
|
mpz_sub (t2, t2, t1); /* trick: recover t1 value just overwritten */ |
|
} |
|
} |
|
|
|
static void |
#if HAVE_NATIVE_mpn_addsub_n |
#if __STDC__ |
xp[size] = mpn_lshift (xp, xp, size, 1); |
mpz_fib_bigcase (mpz_t t1, mpz_t t2, unsigned long int n) |
yp[size] = 0; |
|
ASSERT_NOCARRY (mpn_addsub_n (xp, yp, xp, yp, size+1)); |
|
xsize = size + (xp[size] != 0); |
|
ysize = size + (yp[size] != 0); |
#else |
#else |
mpz_fib_bigcase (t1, t2, n) |
c2 = mpn_lshift (fp, xp, size, 1); |
mpz_t t1; |
c = c2 + mpn_add_n (xp, fp, yp, size); |
mpz_t t2; |
xp[size] = c; |
unsigned long int n; |
xsize = size + (c != 0); |
|
c2 -= mpn_sub_n (yp, fp, yp, size); |
|
yp[size] = c2; |
|
ASSERT (c2 <= 1); |
|
ysize = size + c2; |
#endif |
#endif |
{ |
|
unsigned long int n2; |
|
int ni, i; |
|
mpz_t x1, x2, u1, u2; |
|
|
|
ni = 0; |
size = xsize + ysize; |
for (n2 = n; n2 >= FIB_THRESHOLD; n2 /= 2) |
c = mpn_mul (fp, xp, xsize, yp, ysize); |
ni++; |
|
|
|
mpz_fib_basecase (t1, t2, n2); |
#if BITS_PER_MP_LIMB >= BITS_PER_ULONG |
|
/* no overflow, see comments above */ |
mpz_init (x1); |
ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= MP_LIMB_T_MAX-2); |
mpz_init (x2); |
fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2)); |
mpz_init (u1); |
#else |
mpz_init (u2); |
/* this code only for testing with small limbs, limb<ulong is unusual */ |
|
if (n & 2) |
for (i = ni - 1; i >= 0; i--) |
{ |
|
ASSERT (fp[0] >= 2); |
|
fp[0] -= 2; |
|
} |
|
else |
|
{ |
|
ASSERT (c != MP_LIMB_T_MAX); /* because it's the high of a mul */ |
|
c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2)); |
|
fp[size-1] = c; |
|
} |
|
#endif |
|
} |
|
else |
{ |
{ |
mpz_mul_2exp (x1, t1, 1); |
/* F[2k] = F[k]*(F[k]+2F[k-1]) */ |
mpz_mul_2exp (x2, t2, 1); |
|
|
|
mpz_add (x1, x1, t2); |
mp_size_t xsize, ysize; |
mpz_sub (x2, x2, t1); |
c = mpn_lshift (yp, yp, size, 1); |
|
c += mpn_add_n (yp, yp, xp, size); |
|
yp[size] = c; |
|
xsize = size; |
|
ysize = size + (c != 0); |
|
size += ysize; |
|
c = mpn_mul (fp, yp, ysize, xp, xsize); |
|
} |
|
|
mpz_mul (u1, t2, x1); |
/* one or two high zeros */ |
mpz_mul (u2, t1, x2); |
size -= (c == 0); |
|
size -= (fp[size-1] == 0); |
|
SIZ(fn) = size; |
|
|
if (((n >> i) & 1) == 0) |
TRACE (printf ("done special, size=%ld\n", size); |
{ |
mpn_trace ("fp ", fp, size)); |
mpz_sub (t1, u1, u2); |
|
mpz_set (t2, u1); |
|
} |
|
else |
|
{ |
|
mpz_set (t1, u1); |
|
mpz_mul_2exp (t2, u1, 1); |
|
mpz_sub (t2, t2, u2); |
|
} |
|
} |
|
|
|
mpz_clear (x1); |
TMP_FREE (marker); |
mpz_clear (x2); |
|
mpz_clear (u1); |
|
mpz_clear (u2); |
|
} |
} |