1 /* mpz_fib_ui -- calculate Fibonacci numbers.
3 Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
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.
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.
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/. */
26 /* change to "#define TRACE(x) x" to get some traces */
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.)
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.
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. */
47 mpz_fib_ui (mpz_ptr fn
, unsigned long n
)
50 mp_size_t size
, xalloc
;
55 if (n
<= FIB_TABLE_LIMIT
)
57 PTR(fn
)[0] = FIB_TABLE (n
);
58 SIZ(fn
) = (n
!= 0); /* F[0]==0, others are !=0 */
63 xalloc
= MPN_FIB2_SIZE (n2
) + 1;
64 MPZ_REALLOC (fn
, 2*xalloc
+1);
68 TMP_ALLOC_LIMBS_2 (xp
,xalloc
, yp
,xalloc
);
69 size
= mpn_fib2_ui (xp
, yp
, n2
);
71 TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n",
73 mpn_trace ("xp", xp
, size
);
74 mpn_trace ("yp", yp
, size
));
78 /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k */
79 mp_size_t xsize
, ysize
;
81 #if HAVE_NATIVE_mpn_add_n_sub_n
82 xp
[size
] = mpn_lshift (xp
, xp
, size
, 1);
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);
88 c2
= mpn_lshift (fp
, xp
, size
, 1);
89 c
= c2
+ mpn_add_n (xp
, fp
, yp
, size
);
91 xsize
= size
+ (c
!= 0);
92 c2
-= mpn_sub_n (yp
, fp
, yp
, size
);
99 c
= mpn_mul (fp
, xp
, xsize
, yp
, ysize
);
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));
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));
121 /* F[2k] = F[k]*(F[k]+2F[k-1]) */
123 mp_size_t xsize
, ysize
;
124 c
= mpn_lshift (yp
, yp
, size
, 1);
125 c
+= mpn_add_n (yp
, yp
, xp
, size
);
128 ysize
= size
+ (c
!= 0);
130 c
= mpn_mul (fp
, yp
, ysize
, xp
, xsize
);
133 /* one or two high zeros */
135 size
-= (fp
[size
-1] == 0);
138 TRACE (printf ("done special, size=%ld\n", size
);
139 mpn_trace ("fp ", fp
, size
));