1 /* mpn_fib2m -- calculate Fibonacci numbers, modulo m.
2
3 Contributed to the GNU project by Marco Bodrato.
4
5 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
6 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
7 FUTURE GNU MP RELEASES.
8
9 Copyright 2001, 2002, 2005, 2009, 2018, 2022 Free Software Foundation, Inc.
10
11 This file is part of the GNU MP Library.
12
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of either:
15
16 * the GNU Lesser General Public License as published by the Free
17 Software Foundation; either version 3 of the License, or (at your
18 option) any later version.
19
20 or
21
22 * the GNU General Public License as published by the Free Software
23 Foundation; either version 2 of the License, or (at your option) any
24 later version.
25
26 or both in parallel, as here.
27
28 The GNU MP Library is distributed in the hope that it will be useful, but
29 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
30 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
31 for more details.
32
33 You should have received copies of the GNU General Public License and the
34 GNU Lesser General Public License along with the GNU MP Library. If not,
35 see https://www.gnu.org/licenses/. */
36
37 #include <stdio.h>
38 #include "gmp-impl.h"
39
40
41 #if ! HAVE_NATIVE_mpn_rsblsh1_n && ! HAVE_NATIVE_mpn_sublsh1_n
42 /* Stores |{ap,n}-{bp,n}| in {rp,n},
43 returns the sign of {ap,n}-{bp,n}. */
44 static int
45 abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
46 {
47 mp_limb_t x, y;
48 while (--n >= 0)
49 {
50 x = ap[n];
51 y = bp[n];
52 if (x != y)
53 {
54 ++n;
55 if (x > y)
56 {
57 ASSERT_NOCARRY (mpn_sub_n (rp, ap, bp, n));
58 return 1;
59 }
60 else
61 {
62 ASSERT_NOCARRY (mpn_sub_n (rp, bp, ap, n));
63 return -1;
64 }
65 }
66 rp[n] = 0;
67 }
68 return 0;
69 }
70 #endif
71
72 /* Computes at most count terms of the sequence needed by the
73 Lucas-Lehmer-Riesel test, indexing backward:
74 L_i = L_{i+1}^2 - 2
75
76 The sequence is computed modulo M = {mp, mn}.
77 The starting point is given in L_{count+1} = {lp, mn}.
78 The scratch pointed by sp, needs a space of at least 3 * mn + 1 limbs.
79
80 Returns the index i>0 if L_i = 0 (mod M) is found within the
81 computed count terms of the sequence. Otherwise it returns zero.
82
83 Note: (+/-2)^2-2=2, (+/-1)^2-2=-1, 0^2-2=-2
84 */
85
86 static mp_bitcnt_t
87 mpn_llriter (mp_ptr lp, mp_srcptr mp, mp_size_t mn, mp_bitcnt_t count, mp_ptr sp)
88 {
89 do
90 {
91 mpn_sqr (sp, lp, mn);
92 mpn_tdiv_qr (sp + 2 * mn, lp, 0, sp, 2 * mn, mp, mn);
93 if (lp[0] < 5)
94 {
95 /* If L^2 % M < 5, |L^2 % M - 2| <= 2 */
96 if (mn == 1 || mpn_zero_p (lp + 1, mn - 1))
97 return (lp[0] == 2) ? count : 0;
98 else
99 MPN_DECR_U (lp, mn, 2);
100 }
101 else
102 lp[0] -= 2;
103 } while (--count != 0);
104 return 0;
105 }
106
107 /* Store the Lucas' number L[n] at lp (maybe), computed modulo m. lp
108 and scratch should have room for mn*2+1 limbs.
109
110 Returns the size of L[n] normally.
111
112 If F[n] is zero modulo m, or L[n] is, returns 0 and lp is
113 undefined.
114 */
115
116 static mp_size_t
117 mpn_lucm (mp_ptr lp, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn, mp_ptr scratch)
118 {
119 int neg;
120 mp_limb_t cy;
121
122 ASSERT (! MPN_OVERLAP_P (lp, MAX(2*mn+1,5), scratch, MAX(2*mn+1,5)));
123 ASSERT (nn > 0);
124
125 neg = mpn_fib2m (lp, scratch, np, nn, mp, mn);
126
127 /* F[n] = +/-{lp, mn}, F[n-1] = +/-{scratch, mn} */
128 if (mpn_zero_p (lp, mn))
129 return 0;
130
131 if (neg) /* One sign is opposite, use sub instead of add. */
132 {
133 #if HAVE_NATIVE_mpn_rsblsh1_n || HAVE_NATIVE_mpn_sublsh1_n
134 #if HAVE_NATIVE_mpn_rsblsh1_n
135 cy = mpn_rsblsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */
136 #else
137 cy = mpn_sublsh1_n (lp, lp, scratch, mn); /* L[n] = -/+(F[n]-(-2F[n-1])) */
138 if (cy != 0)
139 cy = mpn_add_n (lp, lp, mp, mn) - cy;
140 #endif
141 if (cy > 1)
142 cy += mpn_add_n (lp, lp, mp, mn);
143 #else
144 cy = mpn_lshift (scratch, scratch, mn, 1); /* 2F[n-1] */
145 if (UNLIKELY (cy))
146 cy -= mpn_sub_n (lp, scratch, lp, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */
147 else
148 abs_sub_n (lp, lp, scratch, mn);
149 #endif
150 ASSERT (cy <= 1);
151 }
152 else
153 {
154 #if HAVE_NATIVE_mpn_addlsh1_n
155 cy = mpn_addlsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]+F[n])) */
156 #else
157 cy = mpn_lshift (scratch, scratch, mn, 1);
158 cy+= mpn_add_n (lp, lp, scratch, mn);
159 #endif
160 ASSERT (cy <= 2);
161 }
162 while (cy || mpn_cmp (lp, mp, mn) >= 0)
163 cy -= mpn_sub_n (lp, lp, mp, mn);
164 MPN_NORMALIZE (lp, mn);
165 return mn;
166 }
167
168 int
169 mpn_strongfibo (mp_srcptr mp, mp_size_t mn, mp_ptr scratch)
170 {
171 mp_ptr lp, sp;
172 mp_size_t en;
173 mp_bitcnt_t b0;
174 TMP_DECL;
175
176 #if GMP_NUMB_BITS % 4 == 0
177 b0 = mpn_scan0 (mp, 0);
178 #else
179 {
180 mpz_t m = MPZ_ROINIT_N(mp, mn);
181 b0 = mpz_scan0 (m, 0);
182 }
183 if (UNLIKELY (b0 == mn * GMP_NUMB_BITS))
184 {
185 en = 1;
186 scratch [0] = 1;
187 }
188 else
189 #endif
190 {
191 int cnt = b0 % GMP_NUMB_BITS;
192 en = b0 / GMP_NUMB_BITS;
193 if (LIKELY (cnt != 0))
194 mpn_rshift (scratch, mp + en, mn - en, cnt);
195 else
196 MPN_COPY (scratch, mp + en, mn - en);
197 en = mn - en;
198 scratch [0] |= 1;
199 en -= scratch [en - 1] == 0;
200 }
201 TMP_MARK;
202
203 lp = TMP_ALLOC_LIMBS (4 * mn + 6);
204 sp = lp + 2 * mn + 3;
205 en = mpn_lucm (sp, scratch, en, mp, mn, lp);
206 if (en != 0 && LIKELY (--b0 != 0))
207 {
208 mpn_sqr (lp, sp, en);
209 lp [0] |= 2; /* V^2 + 2 */
210 if (LIKELY (2 * en >= mn))
211 mpn_tdiv_qr (sp, lp, 0, lp, 2 * en, mp, mn);
212 else
213 MPN_ZERO (lp + 2 * en, mn - 2 * en);
214 if (! mpn_zero_p (lp, mn) && LIKELY (--b0 != 0))
215 b0 = mpn_llriter (lp, mp, mn, b0, lp + mn + 1);
216 }
217 TMP_FREE;
218 return (b0 != 0);
219 }