lucnum2_ui.c revision 1.1.1.2
1/* mpz_lucnum2_ui -- calculate Lucas numbers. 2 3Copyright 2001, 2003, 2005, 2012 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 25void 26mpz_lucnum2_ui (mpz_ptr ln, mpz_ptr lnsub1, unsigned long n) 27{ 28 mp_ptr lp, l1p, f1p; 29 mp_size_t size; 30 mp_limb_t c; 31 TMP_DECL; 32 33 ASSERT (ln != lnsub1); 34 35 /* handle small n quickly, and hide the special case for L[-1]=-1 */ 36 if (n <= FIB_TABLE_LUCNUM_LIMIT) 37 { 38 mp_limb_t f = FIB_TABLE (n); 39 mp_limb_t f1 = FIB_TABLE ((int) n - 1); 40 41 /* L[n] = F[n] + 2F[n-1] */ 42 PTR(ln)[0] = f + 2*f1; 43 SIZ(ln) = 1; 44 45 /* L[n-1] = 2F[n] - F[n-1], but allow for L[-1]=-1 */ 46 PTR(lnsub1)[0] = (n == 0 ? 1 : 2*f - f1); 47 SIZ(lnsub1) = (n == 0 ? -1 : 1); 48 49 return; 50 } 51 52 TMP_MARK; 53 size = MPN_FIB2_SIZE (n); 54 f1p = TMP_ALLOC_LIMBS (size); 55 56 lp = MPZ_REALLOC (ln, size+1); 57 l1p = MPZ_REALLOC (lnsub1, size+1); 58 59 size = mpn_fib2_ui (l1p, f1p, n); 60 61 /* L[n] = F[n] + 2F[n-1] */ 62#if HAVE_NATIVE_mpn_addlsh1_n 63 c = mpn_addlsh1_n (lp, l1p, f1p, size); 64#else 65 c = mpn_lshift (lp, f1p, size, 1); 66 c += mpn_add_n (lp, lp, l1p, size); 67#endif 68 lp[size] = c; 69 SIZ(ln) = size + (c != 0); 70 71 /* L[n-1] = 2F[n] - F[n-1] */ 72 c = mpn_lshift (l1p, l1p, size, 1); 73 c -= mpn_sub_n (l1p, l1p, f1p, size); 74 ASSERT ((mp_limb_signed_t) c >= 0); 75 l1p[size] = c; 76 SIZ(lnsub1) = size + (c != 0); 77 78 TMP_FREE; 79} 80