summaryrefslogtreecommitdiff
path: root/gmp/mpz/lucnum_ui.c
diff options
context:
space:
mode:
Diffstat (limited to 'gmp/mpz/lucnum_ui.c')
-rw-r--r--gmp/mpz/lucnum_ui.c177
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);
View Lorry Controller Status