contrib/gmp/mpz/fib_ui.c

   1 /* mpz_fib_ui -- calculate Fibonacci numbers.
   2
   3 Copyright 2000, 2001, 2002, 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 #include "longlong.h"
  24
  25
  26 /* change to "#define TRACE(x) x" to get some traces */
  27 #define TRACE(x)
  28
  29
  30 /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low
  31    limb because the result F[2k+1] is an F[4m+3] and such numbers are always
  32    == 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7.  (This is
  33    the same as in mpn_fib2_ui.)
  34
  35    In the F[2k+1] for k even, the +2 won't give a carry out of the low limb
  36    in normal circumstances.  This is an F[4m+1] and we claim that F[3*2^b+1]
  37    == 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence
  38    if n < 2^GMP_NUMB_BITS then F[n] cannot have a low limb of 0 or 1.  No
  39    proof for this claim, but it's been verified up to b==32 and has such a
  40    nice pattern it must be true :-).  Of interest is that F[3*2^b] == 0 mod
  41    2^(b+1) seems to hold too.
  42
  43    When n >= 2^GMP_NUMB_BITS, which can arise in a nails build, then the low
  44    limb of F[4m+1] can certainly be 1, and an mpn_add_1 must be used.  */
  45
  46 void
  47 mpz_fib_ui (mpz_ptr fn, unsigned long n)
  48 {
  49   mp_ptr         fp, xp, yp;
  50   mp_size_t      size, xalloc;
  51   unsigned long  n2;
  52   mp_limb_t      c, c2;
  53   TMP_DECL;
  54
  55   if (n <= FIB_TABLE_LIMIT)
  56     {
  57       PTR(fn)[0] = FIB_TABLE (n);
  58       SIZ(fn) = (n != 0);      /* F[0]==0, others are !=0 */
  59       return;
  60     }
  61
  62   n2 = n/2;
  63   xalloc = MPN_FIB2_SIZE (n2) + 1;
  64   MPZ_REALLOC (fn, 2*xalloc+1);
  65   fp = PTR (fn);
  66
  67   TMP_MARK;
  68   TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc);
  69   size = mpn_fib2_ui (xp, yp, n2);
  70
  71   TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n",
  72                  n >> 1, size, n&1);
  73          mpn_trace ("xp", xp, size);
  74          mpn_trace ("yp", yp, size));
  75
  76   if (n & 1)
  77     {
  78       /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k  */
  79       mp_size_t  xsize, ysize;
  80
  81 #if HAVE_NATIVE_mpn_add_n_sub_n
  82       xp[size] = mpn_lshift (xp, xp, size, 1);
  83       yp[size] = 0;
  84       ASSERT_NOCARRY (mpn_add_n_sub_n (xp, yp, xp, yp, size+1));
  85       xsize = size + (xp[size] != 0);
  86       ysize = size + (yp[size] != 0);
  87 #else
  88       c2 = mpn_lshift (fp, xp, size, 1);
  89       c = c2 + mpn_add_n (xp, fp, yp, size);
  90       xp[size] = c;
  91       xsize = size + (c != 0);
  92       c2 -= mpn_sub_n (yp, fp, yp, size);
  93       yp[size] = c2;
  94       ASSERT (c2 <= 1);
  95       ysize = size + c2;
  96 #endif
  97
  98       size = xsize + ysize;
  99       c = mpn_mul (fp, xp, xsize, yp, ysize);
 100
 101 #if GMP_NUMB_BITS >= BITS_PER_ULONG
 102       /* no overflow, see comments above */
 103       ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= GMP_NUMB_MAX-2);
 104       fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2));
 105 #else
 106       if (n & 2)
 107         {
 108           ASSERT (fp[0] >= 2);
 109           fp[0] -= 2;
 110         }
 111       else
 112         {
 113           ASSERT (c != GMP_NUMB_MAX); /* because it's the high of a mul */
 114           c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2));
 115           fp[size-1] = c;
 116         }
 117 #endif
 118     }
 119   else
 120     {
 121       /* F[2k] = F[k]*(F[k]+2F[k-1]) */
 122
 123       mp_size_t  xsize, ysize;
 124       c = mpn_lshift (yp, yp, size, 1);
 125       c += mpn_add_n (yp, yp, xp, size);
 126       yp[size] = c;
 127       xsize = size;
 128       ysize = size + (c != 0);
 129       size += ysize;
 130       c = mpn_mul (fp, yp, ysize, xp, xsize);
 131     }
 132
 133   /* one or two high zeros */
 134   size -= (c == 0);
 135   size -= (fp[size-1] == 0);
 136   SIZ(fn) = size;
 137
 138   TRACE (printf ("done special, size=%ld\n", size);
 139          mpn_trace ("fp ", fp, size));
 140
 141   TMP_FREE;
 142 }