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 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
limb_apprsqrt(mp_limb_t x)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) + (x >> s)) >> 1;
56 }
57
58 /* Performs strong Lucas' test on x, with parameters suggested */
59 /* for the BPSW test. Qk and V are passed to recycle variables. */
60 /* Requires GCD (x,6) = 1.*/
61 int
mpz_stronglucas(mpz_srcptr x,mpz_ptr V,mpz_ptr Qk)62 mpz_stronglucas (mpz_srcptr x, mpz_ptr V, mpz_ptr Qk)
63 {
64 mp_bitcnt_t b0;
65 mpz_t n;
66 mp_limb_t D; /* The absolute value is stored. */
67 long Q;
68 mpz_t T1, T2;
69
70 /* Test on the absolute value. */
71 mpz_roinit_n (n, PTR (x), ABSIZ (x));
72
73 ASSERT (mpz_odd_p (n));
74 /* ASSERT (mpz_gcd_ui (NULL, n, 6) == 1); */
75 #if GMP_NUMB_BITS % 16 == 0
76 /* (2^12 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1) */
77 D = mpn_mod_34lsub1 (PTR (n), SIZ (n));
78 /* (2^12 - 1) = 3^2 * 5 * 7 * 13 */
79 ASSERT (D % 3 != 0 && D % 5 != 0 && D % 7 != 0);
80 if ((D % 5 & 2) != 0)
81 /* (5/n) = -1, iff n = 2 or 3 (mod 5) */
82 /* D = 5; Q = -1 */
83 return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
84 else if (! POW2_P (D % 7))
85 /* (-7/n) = -1, iff n = 3,5 or 6 (mod 7) */
86 D = 7; /* Q = 2 */
87 /* (9/n) = -1, never: 9 = 3^2 */
88 else if (mpz_kronecker_ui (n, 11) == -1)
89 /* (-11/n) = (n/11) */
90 D = 11; /* Q = 3 */
91 else if ((((D % 13 - (D % 13 >> 3)) & 7) > 4) ||
92 (((D % 13 - (D % 13 >> 3)) & 7) == 2))
93 /* (13/n) = -1, iff n = 2,5,6,7,8 or 11 (mod 13) */
94 D = 13; /* Q = -3 */
95 else if (D % 3 == 2)
96 /* (-15/n) = (n/15) = (n/5)*(n/3) */
97 /* Here, (n/5) = 1, and */
98 /* (n/3) = -1, iff n = 2 (mod 3) */
99 D = 15; /* Q = 4 */
100 #if GMP_NUMB_BITS % 32 == 0
101 /* (2^24 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1) */
102 /* (2^24 - 1) = (2^12 - 1) * 17 * 241 */
103 else if (! POW2_P (D % 17) && ! POW2_P (17 - D % 17))
104 D = 17; /* Q = -4 */
105 #endif
106 #else
107 if (mpz_kronecker_ui (n, 5) == -1)
108 return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
109 #endif
110 else
111 {
112 mp_limb_t tl;
113 mp_limb_t maxD;
114 int jac_bit1;
115
116 if (UNLIKELY (mpz_perfect_square_p (n)))
117 return 0; /* A square is composite. */
118
119 /* Check Ds up to square root (in case, n is prime)
120 or avoid overflows */
121 if (SIZ (n) == 1)
122 maxD = limb_apprsqrt (* PTR (n));
123 else if (BITS_PER_ULONG >= GMP_NUMB_BITS && SIZ (n) == 2)
124 mpn_sqrtrem (&maxD, (mp_ptr) NULL, PTR (n), 2);
125 else
126 maxD = GMP_NUMB_MAX;
127 maxD = MIN (maxD, ULONG_MAX);
128
129 D = GMP_NUMB_BITS % 16 == 0 ? (GMP_NUMB_BITS % 32 == 0 ? 17 : 15) : 5;
130
131 /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */
132 /* For those Ds we have (D/n) = (n/|D|) */
133 /* FIXME: Should we loop only on prime Ds? */
134 /* The only interesting composite D is 15. */
135 do
136 {
137 if (UNLIKELY (D >= maxD))
138 return 1;
139 D += 2;
140 jac_bit1 = 0;
141 JACOBI_MOD_OR_MODEXACT_1_ODD (jac_bit1, tl, PTR (n), SIZ (n), D);
142 if (UNLIKELY (tl == 0))
143 return 0;
144 }
145 while (mpn_jacobi_base (tl, D, jac_bit1) == 1);
146 }
147
148 /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
149 Q = (D & 2) ? (D >> 2) + 1 : -(long) (D >> 2);
150 /* ASSERT (mpz_si_kronecker ((D & 2) ? NEG_CAST (long, D) : D, n) == -1); */
151
152 /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */
153 b0 = mpz_scan0 (n, 0);
154
155 mpz_init (T1);
156 mpz_init (T2);
157
158 /* If Ud != 0 && Vd != 0 */
159 if (mpz_lucas_mod (V, Qk, Q, b0, n, T1, T2) == 0)
160 if (LIKELY (--b0 != 0))
161 do
162 {
163 /* V_{2k} <- V_k ^ 2 - 2Q^k */
164 mpz_mul (T2, V, V);
165 mpz_submul_ui (T2, Qk, 2);
166 mpz_tdiv_r (V, T2, n);
167 if (SIZ (V) == 0 || UNLIKELY (--b0 == 0))
168 break;
169 /* Q^{2k} = (Q^k)^2 */
170 mpz_mul (T2, Qk, Qk);
171 mpz_tdiv_r (Qk, T2, n);
172 } while (1);
173
174 mpz_clear (T1);
175 mpz_clear (T2);
176
177 return (b0 != 0);
178 }
179