1/* mpz_lucnum_ui -- calculate Lucas number. 2 3Copyright 2001, 2003, 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 24 25/* change this to "#define TRACE(x) x" for diagnostics */ 26#define TRACE(x) 27 28 29/* Notes: 30 31 For the +4 in L[2k+1] when k is even, all L[4m+3] == 4, 5 or 7 mod 8, so 32 there can't be an overflow applying +4 to just the low limb (since that 33 would leave 0, 1, 2 or 3 mod 8). 34 35 For the -4 in L[2k+1] when k is even, it seems (no proof) that 36 L[3*2^(b-2)-3] == -4 mod 2^b, so for instance with a 32-bit limb 37 L[0xBFFFFFFD] == 0xFFFFFFFC mod 2^32, and this implies a borrow from the 38 low limb. Obviously L[0xBFFFFFFD] is a huge number, but it's at least 39 conceivable to calculate it, so it probably should be handled. 40 41 For the -2 in L[2k] with k even, it seems (no proof) L[2^(b-1)] == -1 mod 42 2^b, so for instance in 32-bits L[0x80000000] has a low limb of 43 0xFFFFFFFF so there would have been a borrow. Again L[0x80000000] is 44 obviously huge, but probably should be made to work. */ 45 46void 47mpz_lucnum_ui (mpz_ptr ln, unsigned long n) 48{ 49 mp_size_t lalloc, xalloc, lsize, xsize; 50 mp_ptr lp, xp; 51 mp_limb_t c; 52 int zeros; 53 TMP_DECL; 54 55 TRACE (printf ("mpn_lucnum_ui n=%lu\n", n)); 56 57 if (n <= FIB_TABLE_LUCNUM_LIMIT) 58 { 59 /* L[n] = F[n] + 2F[n-1] */ 60 PTR(ln)[0] = FIB_TABLE(n) + 2 * FIB_TABLE ((int) n - 1); 61 SIZ(ln) = 1; 62 return; 63 } 64 65 /* +1 since L[n]=F[n]+2F[n-1] might be 1 limb bigger than F[n], further +1 66 since square or mul used below might need an extra limb over the true 67 size */ 68 lalloc = MPN_FIB2_SIZE (n) + 2; 69 MPZ_REALLOC (ln, lalloc); 70 lp = PTR (ln); 71 72 TMP_MARK; 73 xalloc = lalloc; 74 xp = TMP_ALLOC_LIMBS (xalloc); 75 76 /* Strip trailing zeros from n, until either an odd number is reached 77 where the L[2k+1] formula can be used, or until n fits within the 78 FIB_TABLE data. The table is preferred of course. */ 79 zeros = 0; 80 for (;;) 81 { 82 if (n & 1) 83 { 84 /* L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k */ 85 86 mp_size_t yalloc, ysize; 87 mp_ptr yp; 88 89 TRACE (printf (" initial odd n=%lu\n", n)); 90 91 yalloc = MPN_FIB2_SIZE (n/2); 92 yp = TMP_ALLOC_LIMBS (yalloc); 93 ASSERT (xalloc >= yalloc); 94 95 xsize = mpn_fib2_ui (xp, yp, n/2); 96 97 /* possible high zero on F[k-1] */ 98 ysize = xsize; 99 ysize -= (yp[ysize-1] == 0); 100 ASSERT (yp[ysize-1] != 0); 101 102 /* xp = 2*F[k] + F[k-1] */ 103#if HAVE_NATIVE_mpn_addlsh1_n 104 c = mpn_addlsh1_n (xp, yp, xp, xsize); 105#else 106 c = mpn_lshift (xp, xp, xsize, 1); 107 c += mpn_add_n (xp, xp, yp, xsize); 108#endif 109 ASSERT (xalloc >= xsize+1); 110 xp[xsize] = c; 111 xsize += (c != 0); 112 ASSERT (xp[xsize-1] != 0); 113 114 ASSERT (lalloc >= xsize + ysize); 115 c = mpn_mul (lp, xp, xsize, yp, ysize); 116 lsize = xsize + ysize; 117 lsize -= (c == 0); 118 119 /* lp = 5*lp */ 120#if HAVE_NATIVE_mpn_addlshift 121 c = mpn_addlshift (lp, lp, lsize, 2); 122#else 123 c = mpn_lshift (xp, lp, lsize, 2); 124 c += mpn_add_n (lp, lp, xp, lsize); 125#endif 126 ASSERT (lalloc >= lsize+1); 127 lp[lsize] = c; 128 lsize += (c != 0); 129 130 /* lp = lp - 4*(-1)^k */ 131 if (n & 2) 132 { 133 /* no overflow, see comments above */ 134 ASSERT (lp[0] <= MP_LIMB_T_MAX-4); 135 lp[0] += 4; 136 } 137 else 138 { 139 /* won't go negative */ 140 MPN_DECR_U (lp, lsize, CNST_LIMB(4)); 141 } 142 143 TRACE (mpn_trace (" l",lp, lsize)); 144 break; 145 } 146 147 MP_PTR_SWAP (xp, lp); /* balance the swaps wanted in the L[2k] below */ 148 zeros++; 149 n /= 2; 150 151 if (n <= FIB_TABLE_LUCNUM_LIMIT) 152 { 153 /* L[n] = F[n] + 2F[n-1] */ 154 lp[0] = FIB_TABLE (n) + 2 * FIB_TABLE ((int) n - 1); 155 lsize = 1; 156 157 TRACE (printf (" initial small n=%lu\n", n); 158 mpn_trace (" l",lp, lsize)); 159 break; 160 } 161 } 162 163 for ( ; zeros != 0; zeros--) 164 { 165 /* L[2k] = L[k]^2 + 2*(-1)^k */ 166 167 TRACE (printf (" zeros=%d\n", zeros)); 168 169 ASSERT (xalloc >= 2*lsize); 170 mpn_sqr (xp, lp, lsize); 171 lsize *= 2; 172 lsize -= (xp[lsize-1] == 0); 173 174 /* First time around the loop k==n determines (-1)^k, after that k is 175 always even and we set n=0 to indicate that. */ 176 if (n & 1) 177 { 178 /* L[n]^2 == 0 or 1 mod 4, like all squares, so +2 gives no carry */ 179 ASSERT (xp[0] <= MP_LIMB_T_MAX-2); 180 xp[0] += 2; 181 n = 0; 182 } 183 else 184 { 185 /* won't go negative */ 186 MPN_DECR_U (xp, lsize, CNST_LIMB(2)); 187 } 188 189 MP_PTR_SWAP (xp, lp); 190 ASSERT (lp[lsize-1] != 0); 191 } 192 193 /* should end up in the right spot after all the xp/lp swaps */ 194 ASSERT (lp == PTR(ln)); 195 SIZ(ln) = lsize; 196 197 TMP_FREE; 198} 199