fib_ui.c revision 1.1.1.3 
1/* mpz_fib_ui -- calculate Fibonacci numbers.
2
3Copyright 2000-2002, 2005, 2012, 2014 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.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
57void
58mpz_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}
160