diff options
author | Pedro Alvarez <pedro.alvarez@codethink.co.uk> | 2016-05-27 17:39:31 +0100 |
---|---|---|
committer | Pedro Alvarez <pedro.alvarez@codethink.co.uk> | 2016-05-27 17:53:32 +0100 |
commit | 26c75cf8267919f81a1759c9c965a52c660233f9 (patch) | |
tree | cf2a39cf56c2c8ac45760854413ab233e6263974 /gmp/mpz/lucnum_ui.c | |
parent | 56892c1d217baea02092b51a09bbc924130ca84c (diff) | |
download | gcc-tarball-baserock/pedroalvarez/gcc-5.3.0-gmp432.tar.gz |
go to gmp 4.3.2baserock/pedroalvarez/gcc-5.3.0-gmp432
Diffstat (limited to 'gmp/mpz/lucnum_ui.c')
-rw-r--r-- | gmp/mpz/lucnum_ui.c | 177 |
1 files changed, 83 insertions, 94 deletions
diff --git a/gmp/mpz/lucnum_ui.c b/gmp/mpz/lucnum_ui.c index d1fe6b54ce..1215a04b4d 100644 --- a/gmp/mpz/lucnum_ui.c +++ b/gmp/mpz/lucnum_ui.c @@ -1,32 +1,21 @@ /* mpz_lucnum_ui -- calculate Lucas number. -Copyright 2001, 2003, 2005, 2011, 2012 Free Software Foundation, Inc. +Copyright 2001, 2003, 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 either: - - * the GNU Lesser General Public License as published by the Free - Software Foundation; either version 3 of the License, or (at your - option) any later version. - -or - - * the GNU General Public License as published by the Free Software - Foundation; either version 2 of the License, or (at your option) any - later version. - -or both in parallel, as here. +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 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 General Public License -for more details. +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. -You should have received copies of the GNU General Public License and the -GNU Lesser General Public License along with the GNU MP Library. If not, -see https://www.gnu.org/licenses/. */ +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ #include <stdio.h> #include "gmp.h" @@ -77,7 +66,8 @@ mpz_lucnum_ui (mpz_ptr ln, unsigned long n) since square or mul used below might need an extra limb over the true size */ lalloc = MPN_FIB2_SIZE (n) + 2; - lp = MPZ_REALLOC (ln, lalloc); + MPZ_REALLOC (ln, lalloc); + lp = PTR (ln); TMP_MARK; xalloc = lalloc; @@ -90,85 +80,84 @@ mpz_lucnum_ui (mpz_ptr ln, unsigned long n) for (;;) { if (n & 1) - { - /* L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k */ + { + /* L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k */ - mp_size_t yalloc, ysize; - mp_ptr yp; + mp_size_t yalloc, ysize; + mp_ptr yp; - TRACE (printf (" initial odd n=%lu\n", n)); + TRACE (printf (" initial odd n=%lu\n", n)); - yalloc = MPN_FIB2_SIZE (n/2); - yp = TMP_ALLOC_LIMBS (yalloc); - ASSERT (xalloc >= yalloc); + yalloc = MPN_FIB2_SIZE (n/2); + yp = TMP_ALLOC_LIMBS (yalloc); + ASSERT (xalloc >= yalloc); - xsize = mpn_fib2_ui (xp, yp, n/2); + xsize = mpn_fib2_ui (xp, yp, n/2); - /* possible high zero on F[k-1] */ - ysize = xsize; - ysize -= (yp[ysize-1] == 0); - ASSERT (yp[ysize-1] != 0); + /* possible high zero on F[k-1] */ + ysize = xsize; + ysize -= (yp[ysize-1] == 0); + ASSERT (yp[ysize-1] != 0); - /* xp = 2*F[k] + F[k-1] */ + /* xp = 2*F[k] + F[k-1] */ #if HAVE_NATIVE_mpn_addlsh1_n - c = mpn_addlsh1_n (xp, yp, xp, xsize); + c = mpn_addlsh1_n (xp, yp, xp, xsize); #else - c = mpn_lshift (xp, xp, xsize, 1); - c += mpn_add_n (xp, xp, yp, xsize); + c = mpn_lshift (xp, xp, xsize, 1); + c += mpn_add_n (xp, xp, yp, xsize); #endif - ASSERT (xalloc >= xsize+1); - xp[xsize] = c; - xsize += (c != 0); - ASSERT (xp[xsize-1] != 0); - - ASSERT (lalloc >= xsize + ysize); - c = mpn_mul (lp, xp, xsize, yp, ysize); - lsize = xsize + ysize; - lsize -= (c == 0); - - /* lp = 5*lp */ -#if HAVE_NATIVE_mpn_addlsh2_n - c = mpn_addlsh2_n (lp, lp, lp, lsize); + ASSERT (xalloc >= xsize+1); + xp[xsize] = c; + xsize += (c != 0); + ASSERT (xp[xsize-1] != 0); + + ASSERT (lalloc >= xsize + ysize); + c = mpn_mul (lp, xp, xsize, yp, ysize); + lsize = xsize + ysize; + lsize -= (c == 0); + + /* lp = 5*lp */ +#if HAVE_NATIVE_mpn_addlshift + c = mpn_addlshift (lp, lp, lsize, 2); #else - /* FIXME: Is this faster than mpn_mul_1 ? */ - c = mpn_lshift (xp, lp, lsize, 2); - c += mpn_add_n (lp, lp, xp, lsize); + c = mpn_lshift (xp, lp, lsize, 2); + c += mpn_add_n (lp, lp, xp, lsize); #endif - ASSERT (lalloc >= lsize+1); - lp[lsize] = c; - lsize += (c != 0); - - /* lp = lp - 4*(-1)^k */ - if (n & 2) - { - /* no overflow, see comments above */ - ASSERT (lp[0] <= MP_LIMB_T_MAX-4); - lp[0] += 4; - } - else - { - /* won't go negative */ - MPN_DECR_U (lp, lsize, CNST_LIMB(4)); - } - - TRACE (mpn_trace (" l",lp, lsize)); - break; - } + ASSERT (lalloc >= lsize+1); + lp[lsize] = c; + lsize += (c != 0); + + /* lp = lp - 4*(-1)^k */ + if (n & 2) + { + /* no overflow, see comments above */ + ASSERT (lp[0] <= MP_LIMB_T_MAX-4); + lp[0] += 4; + } + else + { + /* won't go negative */ + MPN_DECR_U (lp, lsize, CNST_LIMB(4)); + } + + TRACE (mpn_trace (" l",lp, lsize)); + break; + } MP_PTR_SWAP (xp, lp); /* balance the swaps wanted in the L[2k] below */ zeros++; n /= 2; if (n <= FIB_TABLE_LUCNUM_LIMIT) - { - /* L[n] = F[n] + 2F[n-1] */ - lp[0] = FIB_TABLE (n) + 2 * FIB_TABLE ((int) n - 1); - lsize = 1; - - TRACE (printf (" initial small n=%lu\n", n); - mpn_trace (" l",lp, lsize)); - break; - } + { + /* L[n] = F[n] + 2F[n-1] */ + lp[0] = FIB_TABLE (n) + 2 * FIB_TABLE ((int) n - 1); + lsize = 1; + + TRACE (printf (" initial small n=%lu\n", n); + mpn_trace (" l",lp, lsize)); + break; + } } for ( ; zeros != 0; zeros--) @@ -178,24 +167,24 @@ mpz_lucnum_ui (mpz_ptr ln, unsigned long n) TRACE (printf (" zeros=%d\n", zeros)); ASSERT (xalloc >= 2*lsize); - mpn_sqr (xp, lp, lsize); + mpn_sqr_n (xp, lp, lsize); lsize *= 2; lsize -= (xp[lsize-1] == 0); /* First time around the loop k==n determines (-1)^k, after that k is - always even and we set n=0 to indicate that. */ + always even and we set n=0 to indicate that. */ if (n & 1) - { - /* L[n]^2 == 0 or 1 mod 4, like all squares, so +2 gives no carry */ - ASSERT (xp[0] <= MP_LIMB_T_MAX-2); - xp[0] += 2; - n = 0; - } + { + /* L[n]^2 == 0 or 1 mod 4, like all squares, so +2 gives no carry */ + ASSERT (xp[0] <= MP_LIMB_T_MAX-2); + xp[0] += 2; + n = 0; + } else - { - /* won't go negative */ - MPN_DECR_U (xp, lsize, CNST_LIMB(2)); - } + { + /* won't go negative */ + MPN_DECR_U (xp, lsize, CNST_LIMB(2)); + } MP_PTR_SWAP (xp, lp); ASSERT (lp[lsize-1] != 0); |