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 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 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
Indexes created Sat Feb 28 05:31:39 UTC 2026