/* mpz_fib_ui -- calculate Fibonacci numbers. Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc. This file is part of the GNU MP Library. The GNU MP Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The GNU MP Library 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 Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ #include #include "mpir.h" #include "gmp-impl.h" #include "longlong.h" /* change to "#define TRACE(x) x" to get some traces */ #define TRACE(x) /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low 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.) 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^GMP_NUMB_BITS 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. When n >= 2^GMP_NUMB_BITS, which can arise in a nails build, then the low limb of F[4m+1] can certainly be 1, and an mpn_add_1 must be used. */ void mpz_fib_ui (mpz_ptr fn, unsigned long n) { mp_ptr fp, xp, yp; mp_size_t size, xalloc; unsigned long n2; mp_limb_t c, c2; TMP_DECL; if (n <= FIB_TABLE_LIMIT) { PTR(fn)[0] = FIB_TABLE (n); SIZ(fn) = (n != 0); /* F[0]==0, others are !=0 */ return; } n2 = n/2; xalloc = MPN_FIB2_SIZE (n2) + 1; MPZ_REALLOC (fn, 2*xalloc+1); fp = PTR (fn); TMP_MARK; TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc); size = mpn_fib2_ui (xp, yp, n2); 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) { /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k */ mp_size_t xsize, ysize; #if HAVE_NATIVE_mpn_sumdiff_n xp[size] = mpn_lshift1 (xp, xp, size); yp[size] = 0; ASSERT_NOCARRY (mpn_sumdiff_n (xp, yp, xp, yp, size+1)); xsize = size + (xp[size] != 0); ysize = size + (yp[size] != 0); #else c2 = mpn_lshift1 (fp, xp, size); c = c2 + mpn_add_n (xp, fp, yp, size); xp[size] = c; xsize = size + (c != 0); c2 -= mpn_sub_n (yp, fp, yp, size); yp[size] = c2; ASSERT (c2 <= 1); ysize = size + c2; #endif size = xsize + ysize; c = mpn_mul (fp, xp, xsize, yp, ysize); #if GMP_NUMB_BITS >= BITS_PER_ULONG /* no overflow, see comments above */ ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= GMP_NUMB_MAX-2); fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2)); #else if (n & 2) { ASSERT (fp[0] >= 2); fp[0] -= 2; } else { ASSERT (c != GMP_NUMB_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 { /* F[2k] = F[k]*(F[k]+2F[k-1]) */ mp_size_t xsize, ysize; c = mpn_lshift1 (yp, yp, size); 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); } /* one or two high zeros */ size -= (c == 0); size -= (fp[size-1] == 0); SIZ(fn) = size; TRACE (printf ("done special, size=%ld\n", size); mpn_trace ("fp ", fp, size)); TMP_FREE; }