contrib/gmp/mpz/lucnum_ui.c

   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_n (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 }