1 /* mpz_stronglucas(n, t1, t2) -- An implementation of the strong Lucas
2 primality test on n, using parameters as suggested by the BPSW test.
3
4 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
5 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
6 FUTURE GNU MP RELEASES.
7
8 Copyright 2018, 2020 Free Software Foundation, Inc.
9
10 Contributed by Marco Bodrato.
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 "gmp-impl.h"
39 #include "longlong.h"
40
41 /* Returns an approximation of the sqare root of x.
42 * It gives:
43 * limb_apprsqrt (x) ^ 2 <= x < (limb_apprsqrt (x)+1) ^ 2
44 * or
45 * x <= limb_apprsqrt (x) ^ 2 <= x * 9/8
46 */
47 static mp_limb_t
48 limb_apprsqrt (mp_limb_t x)
49 {
50 int s;
51
52 ASSERT (x > 2);
53 count_leading_zeros (s, x);
54 s = (GMP_LIMB_BITS - s) >> 1;
55 return ((CNST_LIMB(1) << (s - 1)) + (x >> 1 >> s));
56 }
57
58 static int
59 mpz_oddjacobi_ui (mpz_t b, mp_limb_t a)
60 {
61 mp_limb_t b_rem;
62 int result_bit1;
63
64 ASSERT (a & 1);
65 ASSERT (a > 1);
66 ASSERT (SIZ (b) > 0);
67 ASSERT ((*PTR (b) & 1) == 1);
68
69 result_bit1 = 0;
70 JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, PTR (b), SIZ (b), a);
71 if (UNLIKELY (b_rem == 0))
72 return 0;
73 else
74 return mpn_jacobi_base (b_rem, a, result_bit1);
75 }
76
77
78 /* Performs strong Lucas' test on x, with parameters suggested */
79 /* for the BPSW test. Qk and V are passed to recycle variables. */
80 /* Requires GCD (x,6) = 1.*/
81 int
82 mpz_stronglucas (mpz_srcptr x, mpz_ptr V, mpz_ptr Qk)
83 {
84 mp_bitcnt_t b0;
85 mpz_t n;
86 mp_limb_t D; /* The absolute value is stored. */
87 mp_limb_t g;
88 long Q;
89 mpz_t T1, T2;
90
91 /* Test on the absolute value. */
92 mpz_roinit_n (n, PTR (x), ABSIZ (x));
93
94 ASSERT (mpz_odd_p (n));
95 /* ASSERT (mpz_gcd_ui (NULL, n, 6) == 1); */
96 #if GMP_NUMB_BITS % 16 == 0
97 /* (2^12 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1) */
98 g = mpn_mod_34lsub1 (PTR (n), SIZ (n));
99 /* (2^12 - 1) = 3^2 * 5 * 7 * 13 */
100 ASSERT (g % 3 != 0 && g % 5 != 0 && g % 7 != 0);
101 if ((g % 5 & 2) != 0)
102 /* (5/n) = -1, iff n = 2 or 3 (mod 5) */
103 /* D = 5; Q = -1 */
104 return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
105 else if (! POW2_P (g % 7))
106 /* (-7/n) = -1, iff n = 3,5 or 6 (mod 7) */
107 D = 7; /* Q = 2 */
108 /* (9/n) = -1, never: 9 = 3^2 */
109 else if (mpz_oddjacobi_ui (n, 11) == -1)
110 /* (-11/n) = (n/11) */
111 D = 11; /* Q = 3 */
112 else if ((((g % 13 - (g % 13 >> 3)) & 7) > 4) ||
113 (((g % 13 - (g % 13 >> 3)) & 7) == 2))
114 /* (13/n) = -1, iff n = 2,5,6,7,8 or 11 (mod 13) */
115 D = 13; /* Q = -3 */
116 else if (g % 3 == 2)
117 /* (-15/n) = (n/15) = (n/5)*(n/3) */
118 /* Here, (n/5) = 1, and */
119 /* (n/3) = -1, iff n = 2 (mod 3) */
120 D = 15; /* Q = 4 */
121 #if GMP_NUMB_BITS % 32 == 0
122 /* (2^24 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1) */
123 /* (2^24 - 1) = (2^12 - 1) * 17 * 241 */
124 else if (! POW2_P (g % 17) && ! POW2_P (17 - g % 17))
125 /* (17/n) = -1, iff n != +-1,+-2,+-4,+-8 (mod 17) */
126 D = 17; /* Q = -4 */
127 #endif
128 #else
129 if (mpz_oddjacobi_ui (n, 5) == -1)
130 return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
131 #endif
132 else
133 {
134 mp_limb_t maxD;
135 int jac;
136
137 /* n is odd, to possibly be a square, n % 8 = 1 is needed. */
138 if (((*PTR (n) & 6) == 0) && UNLIKELY (mpz_perfect_square_p (n)))
139 return 0; /* A square is composite. */
140
141 /* Check Ds up to square root (in case, n is prime)
142 or avoid overflows */
143 if (SIZ (n) == 1)
144 maxD = limb_apprsqrt (* PTR (n));
145 else if (BITS_PER_ULONG >= GMP_NUMB_BITS && SIZ (n) == 2)
146 mpn_sqrtrem (&maxD, (mp_ptr) NULL, PTR (n), 2);
147 else
148 maxD = GMP_NUMB_MAX;
149 maxD = MIN (maxD, ULONG_MAX);
150
151 unsigned Ddiff = 2;
152 #if GMP_NUMB_BITS % 16 == 0
153 const unsigned D2 = 6;
154 #if GMP_NUMB_BITS % 32 == 0
155 D = 19;
156 Ddiff = 4;
157 #else
158 D = 17;
159 #endif
160 #else
161 const unsigned D2 = 4;
162 D = 7;
163 #endif
164
165 /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */
166 /* For those Ds we have (D/n) = (n/|D|) */
167 /* FIXME: Should we loop only on prime Ds? */
168 /* The only interesting composite D is 15, because 3 is not tested. */
169 for (;;)
170 {
171 jac = mpz_oddjacobi_ui (n, D);
172 if (jac != 1)
173 break;
174 if (UNLIKELY (D >= maxD))
175 return 1;
176 D += Ddiff;
177 Ddiff = D2 - Ddiff;
178 }
179
180 if (UNLIKELY (jac == 0))
181 return 0;
182 }
183
184 /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
185 Q = (D & 2) ? (D >> 2) + 1 : -(long) (D >> 2);
186 /* ASSERT (mpz_si_kronecker ((D & 2) ? NEG_CAST (long, D) : D, n) == -1); */
187
188 /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */
189 b0 = mpz_scan0 (n, 0);
190
191 mpz_init (T1);
192 mpz_init (T2);
193
194 /* If Ud != 0 && Vd != 0 */
195 if (mpz_lucas_mod (V, Qk, Q, b0, n, T1, T2) == 0)
196 if (LIKELY (--b0 != 0))
197 for (;;)
198 {
199 /* V_{2k} <- V_k ^ 2 - 2Q^k */
200 mpz_mul (T2, V, V);
201 mpz_submul_ui (T2, Qk, 2);
202 mpz_tdiv_r (V, T2, n);
203 if (SIZ (V) == 0 || UNLIKELY (--b0 == 0))
204 break;
205 /* Q^{2k} = (Q^k)^2 */
206 mpz_mul (T2, Qk, Qk);
207 mpz_tdiv_r (Qk, T2, n);
208 }
209
210 mpz_clear (T1);
211 mpz_clear (T2);
212
213 return (b0 != 0);
214 }