1 /* mpn_fib2m -- calculate Fibonacci numbers, modulo m.
2
3 Contributed to the GNU project by Marco Bodrato, based on the previous
4 fib2_ui.c file.
5
6 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
7 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
8 FUTURE GNU MP RELEASES.
9
10 Copyright 2001, 2002, 2005, 2009, 2018 Free Software Foundation, Inc.
11
12 This file is part of the GNU MP Library.
13
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of either:
16
17 * the GNU Lesser General Public License as published by the Free
18 Software Foundation; either version 3 of the License, or (at your
19 option) any later version.
20
21 or
22
23 * the GNU General Public License as published by the Free Software
24 Foundation; either version 2 of the License, or (at your option) any
25 later version.
26
27 or both in parallel, as here.
28
29 The GNU MP Library is distributed in the hope that it will be useful, but
30 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
32 for more details.
33
34 You should have received copies of the GNU General Public License and the
35 GNU Lesser General Public License along with the GNU MP Library. If not,
36 see https://www.gnu.org/licenses/. */
37
38 #include <stdio.h>
39 #include "gmp-impl.h"
40 #include "longlong.h"
41
42
43 /* Stores |{ap,n}-{bp,n}| in {rp,n},
44 returns the sign of {ap,n}-{bp,n}. */
45 static int
46 abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
47 {
48 mp_limb_t x, y;
49 while (--n >= 0)
50 {
51 x = ap[n];
52 y = bp[n];
53 if (x != y)
54 {
55 ++n;
56 if (x > y)
57 {
58 ASSERT_NOCARRY (mpn_sub_n (rp, ap, bp, n));
59 return 1;
60 }
61 else
62 {
63 ASSERT_NOCARRY (mpn_sub_n (rp, bp, ap, n));
64 return -1;
65 }
66 }
67 rp[n] = 0;
68 }
69 return 0;
70 }
71
72 /* Store F[n] at fp and F[n-1] at f1p. Both are computed modulo m.
73 fp and f1p should have room for mn*2+1 limbs.
74
75 The sign of one or both the values may be flipped (n-F, instead of F),
76 the return value is 0 (zero) if the signs are coherent (both positive
77 or both negative) and 1 (one) otherwise.
78
79 Notes:
80
81 In F[2k+1] with k even, +2 is applied to 4*F[k]^2 just by ORing into the
82 low limb.
83
84 In F[2k+1] with k odd, -2 is applied to F[k-1]^2 just by ORing into the
85 low limb.
86
87 TODO: Should {tp, 2 * mn} be passed as a scratch pointer?
88 Should the call to mpn_fib2_ui() obtain (up to) 2*mn limbs?
89 */
90
91 int
92 mpn_fib2m (mp_ptr fp, mp_ptr f1p, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn)
93 {
94 unsigned long nfirst;
95 mp_limb_t nh;
96 mp_bitcnt_t nbi;
97 mp_size_t sn, fn;
98 int fcnt, ncnt;
99
100 ASSERT (! MPN_OVERLAP_P (fp, MAX(2*mn+1,5), f1p, MAX(2*mn+1,5)));
101 ASSERT (nn > 0 && np[nn - 1] != 0);
102
103 /* Estimate the maximal n such that fibonacci(n) fits in mn limbs. */
104 #if GMP_NUMB_BITS % 16 == 0
105 if (UNLIKELY (ULONG_MAX / (23 * (GMP_NUMB_BITS / 16)) <= mn))
106 nfirst = ULONG_MAX;
107 else
108 nfirst = mn * (23 * (GMP_NUMB_BITS / 16));
109 #else
110 {
111 mp_bitcnt_t mbi;
112 mbi = (mp_bitcnt_t) mn * GMP_NUMB_BITS;
113
114 if (UNLIKELY (ULONG_MAX / 23 < mbi))
115 {
116 if (UNLIKELY (ULONG_MAX / 23 * 16 <= mbi))
117 nfirst = ULONG_MAX;
118 else
119 nfirst = mbi / 16 * 23;
120 }
121 else
122 nfirst = mbi * 23 / 16;
123 }
124 #endif
125
126 sn = nn - 1;
127 nh = np[sn];
128 count_leading_zeros (ncnt, nh);
129 count_leading_zeros (fcnt, nfirst);
130
131 if (fcnt >= ncnt)
132 {
133 ncnt = fcnt - ncnt;
134 nh >>= ncnt;
135 }
136 else if (sn > 0)
137 {
138 ncnt -= fcnt;
139 nh <<= ncnt;
140 ncnt = GMP_NUMB_BITS - ncnt;
141 --sn;
142 nh |= np[sn] >> ncnt;
143 }
144 else
145 ncnt = 0;
146
147 nbi = sn * GMP_NUMB_BITS + ncnt;
148 if (nh > nfirst)
149 {
150 nh >>= 1;
151 ++nbi;
152 }
153
154 ASSERT (nh <= nfirst);
155 /* Take a starting pair from mpn_fib2_ui. */
156 fn = mpn_fib2_ui (fp, f1p, nh);
157 MPN_ZERO (fp + fn, mn - fn);
158 MPN_ZERO (f1p + fn, mn - fn);
159
160 if (nbi == 0)
161 {
162 if (fn == mn)
163 {
164 mp_limb_t qp[2];
165 mpn_tdiv_qr (qp, fp, 0, fp, fn, mp, mn);
166 mpn_tdiv_qr (qp, f1p, 0, f1p, fn, mp, mn);
167 }
168
169 return 0;
170 }
171 else
172 {
173 mp_ptr tp;
174 unsigned pb = nh & 1;
175 int neg;
176 TMP_DECL;
177
178 TMP_MARK;
179
180 tp = TMP_ALLOC_LIMBS (2 * mn + (mn < 2));
181
182 do
183 {
184 mp_ptr rp;
185 /* Here fp==F[k] and f1p==F[k-1], with k being the bits of n from
186 nbi upwards.
187
188 Based on the next bit of n, we'll double to the pair
189 fp==F[2k],f1p==F[2k-1] or fp==F[2k+1],f1p==F[2k], according as
190 that bit is 0 or 1 respectively. */
191
192 mpn_sqr (tp, fp, mn);
193 mpn_sqr (fp, f1p, mn);
194
195 /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2. */
196 f1p[2 * mn] = mpn_add_n (f1p, tp, fp, 2 * mn);
197
198 /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k.
199 pb is the low bit of our implied k. */
200
201 /* fp is F[k-1]^2 == 0 or 1 mod 4, like all squares. */
202 ASSERT ((fp[0] & 2) == 0);
203 ASSERT (pb == (pb & 1));
204 ASSERT ((fp[0] + (pb ? 2 : 0)) == (fp[0] | (pb << 1)));
205 fp[0] |= pb << 1; /* possible -2 */
206 #if HAVE_NATIVE_mpn_rsblsh2_n
207 fp[2 * mn] = 1 + mpn_rsblsh2_n (fp, fp, tp, 2 * mn);
208 MPN_INCR_U(fp, 2 * mn + 1, (1 ^ pb) << 1); /* possible +2 */
209 fp[2 * mn] = (fp[2 * mn] - 1) & GMP_NUMB_MAX;
210 #else
211 {
212 mp_limb_t c;
213
214 c = mpn_lshift (tp, tp, 2 * mn, 2);
215 tp[0] |= (1 ^ pb) << 1; /* possible +2 */
216 c -= mpn_sub_n (fp, tp, fp, 2 * mn);
217 fp[2 * mn] = c & GMP_NUMB_MAX;
218 }
219 #endif
220 neg = fp[2 * mn] == GMP_NUMB_MAX;
221
222 /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2 */
223 /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k */
224
225 /* Calculate F[2k] = F[2k+1] - F[2k-1], replacing the unwanted one of
226 F[2k+1] and F[2k-1]. */
227 --nbi;
228 pb = (np [nbi / GMP_NUMB_BITS] >> (nbi % GMP_NUMB_BITS)) & 1;
229 rp = pb ? f1p : fp;
230 if (neg)
231 {
232 /* Calculate -(F[2k+1] - F[2k-1]) */
233 rp[2 * mn] = f1p[2 * mn] + 1 - mpn_sub_n (rp, f1p, fp, 2 * mn);
234 neg = ! pb;
235 if (pb) /* fp not overwritten, negate it. */
236 fp [2 * mn] = 1 ^ mpn_neg (fp, fp, 2 * mn);
237 }
238 else
239 {
240 neg = abs_sub_n (rp, fp, f1p, 2 * mn + 1) < 0;
241 }
242
243 mpn_tdiv_qr (tp, fp, 0, fp, 2 * mn + 1, mp, mn);
244 mpn_tdiv_qr (tp, f1p, 0, f1p, 2 * mn + 1, mp, mn);
245 }
246 while (nbi != 0);
247
248 TMP_FREE;
249
250 return neg;
251 }
252 }