1/* mpz_fib_ui -- calculate Fibonacci numbers.
2
3Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library.
6
7The GNU MP Library is free software; you can redistribute it and/or modify
8it under the terms of the GNU Lesser General Public License as published by
9the Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
15License for more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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
46void
47mpz_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}
143