source/libs/gmp/gmp-src/mpz/fib_ui.c

   1 /* mpz_fib_ui -- calculate Fibonacci numbers.
   2
   3 Copyright 2000-2002, 2005, 2012, 2014 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 either:
   9
  10   * the GNU Lesser General Public License as published by the Free
  11     Software Foundation; either version 3 of the License, or (at your
  12     option) any later version.
  13
  14 or
  15
  16   * the GNU General Public License as published by the Free Software
  17     Foundation; either version 2 of the License, or (at your option) any
  18     later version.
  19
  20 or both in parallel, as here.
  21
  22 The GNU MP Library is distributed in the hope that it will be useful, but
  23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  25 for more details.
  26
  27 You should have received copies of the GNU General Public License and the
  28 GNU Lesser General Public License along with the GNU MP Library.  If not,
  29 see https://www.gnu.org/licenses/.  */
  30
  31 #include <stdio.h>
  32 #include "gmp.h"
  33 #include "gmp-impl.h"
  34 #include "longlong.h"
  35
  36
  37 /* change to "#define TRACE(x) x" to get some traces */
  38 #define TRACE(x)
  39
  40
  41 /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low
  42    limb because the result F[2k+1] is an F[4m+3] and such numbers are always
  43    == 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7.  (This is
  44    the same as in mpn_fib2_ui.)
  45
  46    In the F[2k+1] for k even, the +2 won't give a carry out of the low limb
  47    in normal circumstances.  This is an F[4m+1] and we claim that F[3*2^b+1]
  48    == 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence
  49    if n < 2^GMP_NUMB_BITS then F[n] cannot have a low limb of 0 or 1.  No
  50    proof for this claim, but it's been verified up to b==32 and has such a
  51    nice pattern it must be true :-).  Of interest is that F[3*2^b] == 0 mod
  52    2^(b+1) seems to hold too.
  53
  54    When n >= 2^GMP_NUMB_BITS, which can arise in a nails build, then the low
  55    limb of F[4m+1] can certainly be 1, and an mpn_add_1 must be used.  */
  56
  57 void
  58 mpz_fib_ui (mpz_ptr fn, unsigned long n)
  59 {
  60   mp_ptr         fp, xp, yp;
  61   mp_size_t      size, xalloc;
  62   unsigned long  n2;
  63   mp_limb_t      c;
  64   TMP_DECL;
  65
  66   if (n <= FIB_TABLE_LIMIT)
  67     {
  68       PTR(fn)[0] = FIB_TABLE (n);
  69       SIZ(fn) = (n != 0);      /* F[0]==0, others are !=0 */
  70       return;
  71     }
  72
  73   n2 = n/2;
  74   xalloc = MPN_FIB2_SIZE (n2) + 1;
  75   fp = MPZ_NEWALLOC (fn, 2 * xalloc);
  76
  77   TMP_MARK;
  78   TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc);
  79   size = mpn_fib2_ui (xp, yp, n2);
  80
  81   TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n",
  82                  n >> 1, size, n&1);
  83          mpn_trace ("xp", xp, size);
  84          mpn_trace ("yp", yp, size));
  85
  86   if (n & 1)
  87     {
  88       /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k  */
  89       mp_size_t  xsize, ysize;
  90
  91 #if HAVE_NATIVE_mpn_add_n_sub_n
  92       xp[size] = mpn_lshift (xp, xp, size, 1);
  93       yp[size] = 0;
  94       ASSERT_NOCARRY (mpn_add_n_sub_n (xp, yp, xp, yp, size+1));
  95       xsize = size + (xp[size] != 0);
  96       ASSERT (yp[size] <= 1);
  97       ysize = size + yp[size];
  98 #else
  99       mp_limb_t  c2;
 100
 101       c2 = mpn_lshift (fp, xp, size, 1);
 102       c = c2 + mpn_add_n (xp, fp, yp, size);
 103       xp[size] = c;
 104       xsize = size + (c != 0);
 105       c2 -= mpn_sub_n (yp, fp, yp, size);
 106       yp[size] = c2;
 107       ASSERT (c2 <= 1);
 108       ysize = size + c2;
 109 #endif
 110
 111       size = xsize + ysize;
 112       c = mpn_mul (fp, xp, xsize, yp, ysize);
 113
 114 #if GMP_NUMB_BITS >= BITS_PER_ULONG
 115       /* no overflow, see comments above */
 116       ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= GMP_NUMB_MAX-2);
 117       fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2));
 118 #else
 119       if (n & 2)
 120         {
 121           ASSERT (fp[0] >= 2);
 122           fp[0] -= 2;
 123         }
 124       else
 125         {
 126           ASSERT (c != GMP_NUMB_MAX); /* because it's the high of a mul */
 127           c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2));
 128           fp[size-1] = c;
 129         }
 130 #endif
 131     }
 132   else
 133     {
 134       /* F[2k] = F[k]*(F[k]+2F[k-1]) */
 135
 136       mp_size_t  xsize, ysize;
 137 #if HAVE_NATIVE_mpn_addlsh1_n
 138       c = mpn_addlsh1_n (yp, xp, yp, size);
 139 #else
 140       c = mpn_lshift (yp, yp, size, 1);
 141       c += mpn_add_n (yp, yp, xp, size);
 142 #endif
 143       yp[size] = c;
 144       xsize = size;
 145       ysize = size + (c != 0);
 146       size += ysize;
 147       c = mpn_mul (fp, yp, ysize, xp, xsize);
 148     }
 149
 150   /* one or two high zeros */
 151   size -= (c == 0);
 152   size -= (fp[size-1] == 0);
 153   SIZ(fn) = size;
 154
 155   TRACE (printf ("done special, size=%ld\n", size);
 156          mpn_trace ("fp ", fp, size));
 157
 158   TMP_FREE;
 159 }