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