1 /* mpz_lucnum_ui -- calculate Lucas number.
2
3 Copyright 2001, 2003, 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
24
25 /* change this to "#define TRACE(x) x" for diagnostics */
26 #define TRACE(x)
27
28
29 /* Notes:
30
31 For the +4 in L[2k+1] when k is even, all L[4m+3] == 4, 5 or 7 mod 8, so
32 there can't be an overflow applying +4 to just the low limb (since that
33 would leave 0, 1, 2 or 3 mod 8).
34
35 For the -4 in L[2k+1] when k is even, it seems (no proof) that
36 L[3*2^(b-2)-3] == -4 mod 2^b, so for instance with a 32-bit limb
37 L[0xBFFFFFFD] == 0xFFFFFFFC mod 2^32, and this implies a borrow from the
38 low limb. Obviously L[0xBFFFFFFD] is a huge number, but it's at least
39 conceivable to calculate it, so it probably should be handled.
40
41 For the -2 in L[2k] with k even, it seems (no proof) L[2^(b-1)] == -1 mod
42 2^b, so for instance in 32-bits L[0x80000000] has a low limb of
43 0xFFFFFFFF so there would have been a borrow. Again L[0x80000000] is
44 obviously huge, but probably should be made to work. */
45
46 void
mpz_lucnum_ui(mpz_ptr ln,unsigned long n)47 mpz_lucnum_ui (mpz_ptr ln, unsigned long n)
48 {
49 mp_size_t lalloc, xalloc, lsize, xsize;
50 mp_ptr lp, xp;
51 mp_limb_t c;
52 int zeros;
53 TMP_DECL;
54
55 TRACE (printf ("mpn_lucnum_ui n=%lu\n", n));
56
57 if (n <= FIB_TABLE_LUCNUM_LIMIT)
58 {
59 /* L[n] = F[n] + 2F[n-1] */
60 PTR(ln)[0] = FIB_TABLE(n) + 2 * FIB_TABLE ((int) n - 1);
61 SIZ(ln) = 1;
62 return;
63 }
64
65 /* +1 since L[n]=F[n]+2F[n-1] might be 1 limb bigger than F[n], further +1
66 since square or mul used below might need an extra limb over the true
67 size */
68 lalloc = MPN_FIB2_SIZE (n) + 2;
69 MPZ_REALLOC (ln, lalloc);
70 lp = PTR (ln);
71
72 TMP_MARK;
73 xalloc = lalloc;
74 xp = TMP_ALLOC_LIMBS (xalloc);
75
76 /* Strip trailing zeros from n, until either an odd number is reached
77 where the L[2k+1] formula can be used, or until n fits within the
78 FIB_TABLE data. The table is preferred of course. */
79 zeros = 0;
80 for (;;)
81 {
82 if (n & 1)
83 {
84 /* L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k */
85
86 mp_size_t yalloc, ysize;
87 mp_ptr yp;
88
89 TRACE (printf (" initial odd n=%lu\n", n));
90
91 yalloc = MPN_FIB2_SIZE (n/2);
92 yp = TMP_ALLOC_LIMBS (yalloc);
93 ASSERT (xalloc >= yalloc);
94
95 xsize = mpn_fib2_ui (xp, yp, n/2);
96
97 /* possible high zero on F[k-1] */
98 ysize = xsize;
99 ysize -= (yp[ysize-1] == 0);
100 ASSERT (yp[ysize-1] != 0);
101
102 /* xp = 2*F[k] + F[k-1] */
103 #if HAVE_NATIVE_mpn_addlsh1_n
104 c = mpn_addlsh1_n (xp, yp, xp, xsize);
105 #else
106 c = mpn_lshift (xp, xp, xsize, 1);
107 c += mpn_add_n (xp, xp, yp, xsize);
108 #endif
109 ASSERT (xalloc >= xsize+1);
110 xp[xsize] = c;
111 xsize += (c != 0);
112 ASSERT (xp[xsize-1] != 0);
113
114 ASSERT (lalloc >= xsize + ysize);
115 c = mpn_mul (lp, xp, xsize, yp, ysize);
116 lsize = xsize + ysize;
117 lsize -= (c == 0);
118
119 /* lp = 5*lp */
120 #if HAVE_NATIVE_mpn_addlshift
121 c = mpn_addlshift (lp, lp, lsize, 2);
122 #else
123 c = mpn_lshift (xp, lp, lsize, 2);
124 c += mpn_add_n (lp, lp, xp, lsize);
125 #endif
126 ASSERT (lalloc >= lsize+1);
127 lp[lsize] = c;
128 lsize += (c != 0);
129
130 /* lp = lp - 4*(-1)^k */
131 if (n & 2)
132 {
133 /* no overflow, see comments above */
134 ASSERT (lp[0] <= MP_LIMB_T_MAX-4);
135 lp[0] += 4;
136 }
137 else
138 {
139 /* won't go negative */
140 MPN_DECR_U (lp, lsize, CNST_LIMB(4));
141 }
142
143 TRACE (mpn_trace (" l",lp, lsize));
144 break;
145 }
146
147 MP_PTR_SWAP (xp, lp); /* balance the swaps wanted in the L[2k] below */
148 zeros++;
149 n /= 2;
150
151 if (n <= FIB_TABLE_LUCNUM_LIMIT)
152 {
153 /* L[n] = F[n] + 2F[n-1] */
154 lp[0] = FIB_TABLE (n) + 2 * FIB_TABLE ((int) n - 1);
155 lsize = 1;
156
157 TRACE (printf (" initial small n=%lu\n", n);
158 mpn_trace (" l",lp, lsize));
159 break;
160 }
161 }
162
163 for ( ; zeros != 0; zeros--)
164 {
165 /* L[2k] = L[k]^2 + 2*(-1)^k */
166
167 TRACE (printf (" zeros=%d\n", zeros));
168
169 ASSERT (xalloc >= 2*lsize);
170 mpn_sqr (xp, lp, lsize);
171 lsize *= 2;
172 lsize -= (xp[lsize-1] == 0);
173
174 /* First time around the loop k==n determines (-1)^k, after that k is
175 always even and we set n=0 to indicate that. */
176 if (n & 1)
177 {
178 /* L[n]^2 == 0 or 1 mod 4, like all squares, so +2 gives no carry */
179 ASSERT (xp[0] <= MP_LIMB_T_MAX-2);
180 xp[0] += 2;
181 n = 0;
182 }
183 else
184 {
185 /* won't go negative */
186 MPN_DECR_U (xp, lsize, CNST_LIMB(2));
187 }
188
189 MP_PTR_SWAP (xp, lp);
190 ASSERT (lp[lsize-1] != 0);
191 }
192
193 /* should end up in the right spot after all the xp/lp swaps */
194 ASSERT (lp == PTR(ln));
195 SIZ(ln) = lsize;
196
197 TMP_FREE;
198 }
199