1 /* mpz_lucnum_ui -- calculate Lucas number. 2 3 Copyright 2001, 2003, 2005 Free Software Foundation, Inc. 4 5 This file is part of the GNU MP Library. 6 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 3 of the License, or (at your 10 option) any later version. 11 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. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along 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 46 void 47 mpz_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