1/* mpz_fib_ui -- calculate Fibonacci numbers. 2 3Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc. 4 5This file is part of the GNU MP Library. 6 7The GNU MP Library is free software; you can redistribute it and/or modify 8it under the terms of the GNU Lesser General Public License as published by 9the Free Software Foundation; either version 3 of the License, or (at your 10option) any later version. 11 12The GNU MP Library is distributed in the hope that it will be useful, but 13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15License for more details. 16 17You should have received a copy of the GNU Lesser General Public License 18along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 19 20#include <stdio.h> 21#include "gmp.h" 22#include "gmp-impl.h" 23#include "longlong.h" 24 25 26/* change to "#define TRACE(x) x" to get some traces */ 27#define TRACE(x) 28 29 30/* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low 31 limb because the result F[2k+1] is an F[4m+3] and such numbers are always 32 == 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7. (This is 33 the same as in mpn_fib2_ui.) 34 35 In the F[2k+1] for k even, the +2 won't give a carry out of the low limb 36 in normal circumstances. This is an F[4m+1] and we claim that F[3*2^b+1] 37 == 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence 38 if n < 2^GMP_NUMB_BITS then F[n] cannot have a low limb of 0 or 1. No 39 proof for this claim, but it's been verified up to b==32 and has such a 40 nice pattern it must be true :-). Of interest is that F[3*2^b] == 0 mod 41 2^(b+1) seems to hold too. 42 43 When n >= 2^GMP_NUMB_BITS, which can arise in a nails build, then the low 44 limb of F[4m+1] can certainly be 1, and an mpn_add_1 must be used. */ 45 46void 47mpz_fib_ui (mpz_ptr fn, unsigned long n) 48{ 49 mp_ptr fp, xp, yp; 50 mp_size_t size, xalloc; 51 unsigned long n2; 52 mp_limb_t c, c2; 53 TMP_DECL; 54 55 if (n <= FIB_TABLE_LIMIT) 56 { 57 PTR(fn)[0] = FIB_TABLE (n); 58 SIZ(fn) = (n != 0); /* F[0]==0, others are !=0 */ 59 return; 60 } 61 62 n2 = n/2; 63 xalloc = MPN_FIB2_SIZE (n2) + 1; 64 MPZ_REALLOC (fn, 2*xalloc+1); 65 fp = PTR (fn); 66 67 TMP_MARK; 68 TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc); 69 size = mpn_fib2_ui (xp, yp, n2); 70 71 TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n", 72 n >> 1, size, n&1); 73 mpn_trace ("xp", xp, size); 74 mpn_trace ("yp", yp, size)); 75 76 if (n & 1) 77 { 78 /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k */ 79 mp_size_t xsize, ysize; 80 81#if HAVE_NATIVE_mpn_add_n_sub_n 82 xp[size] = mpn_lshift (xp, xp, size, 1); 83 yp[size] = 0; 84 ASSERT_NOCARRY (mpn_add_n_sub_n (xp, yp, xp, yp, size+1)); 85 xsize = size + (xp[size] != 0); 86 ysize = size + (yp[size] != 0); 87#else 88 c2 = mpn_lshift (fp, xp, size, 1); 89 c = c2 + mpn_add_n (xp, fp, yp, size); 90 xp[size] = c; 91 xsize = size + (c != 0); 92 c2 -= mpn_sub_n (yp, fp, yp, size); 93 yp[size] = c2; 94 ASSERT (c2 <= 1); 95 ysize = size + c2; 96#endif 97 98 size = xsize + ysize; 99 c = mpn_mul (fp, xp, xsize, yp, ysize); 100 101#if GMP_NUMB_BITS >= BITS_PER_ULONG 102 /* no overflow, see comments above */ 103 ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= GMP_NUMB_MAX-2); 104 fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2)); 105#else 106 if (n & 2) 107 { 108 ASSERT (fp[0] >= 2); 109 fp[0] -= 2; 110 } 111 else 112 { 113 ASSERT (c != GMP_NUMB_MAX); /* because it's the high of a mul */ 114 c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2)); 115 fp[size-1] = c; 116 } 117#endif 118 } 119 else 120 { 121 /* F[2k] = F[k]*(F[k]+2F[k-1]) */ 122 123 mp_size_t xsize, ysize; 124 c = mpn_lshift (yp, yp, size, 1); 125 c += mpn_add_n (yp, yp, xp, size); 126 yp[size] = c; 127 xsize = size; 128 ysize = size + (c != 0); 129 size += ysize; 130 c = mpn_mul (fp, yp, ysize, xp, xsize); 131 } 132 133 /* one or two high zeros */ 134 size -= (c == 0); 135 size -= (fp[size-1] == 0); 136 SIZ(fn) = size; 137 138 TRACE (printf ("done special, size=%ld\n", size); 139 mpn_trace ("fp ", fp, size)); 140 141 TMP_FREE; 142} 143