Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpz/stronglucas.c
Line
Count
Source (jump to first uncovered line)
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
5
{
50
5
  int s;
51
52
5
  ASSERT (x > 2);
53
5
  count_leading_zeros (s, x);
54
5
  s = (GMP_LIMB_BITS - s) >> 1;
55
5
  return ((CNST_LIMB(1) << (s - 1)) + (x >> 1 >> s));
56
5
}
57
58
static int
59
mpz_oddjacobi_ui (mpz_t b, mp_limb_t a)
60
616
{
61
616
  mp_limb_t  b_rem;
62
616
  int        result_bit1;
63
64
616
  ASSERT (a & 1);
65
616
  ASSERT (a > 1);
66
616
  ASSERT (SIZ (b) > 0);
67
616
  ASSERT ((*PTR (b) & 1) == 1);
68
69
616
  result_bit1 = 0;
70
616
  JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, PTR (b), SIZ (b), a);
71
616
  if (UNLIKELY (b_rem == 0))
72
0
    return 0;
73
616
  else
74
616
    return mpn_jacobi_base (b_rem, a, result_bit1);
75
616
}
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
2.13k
{
84
2.13k
  mp_bitcnt_t b0;
85
2.13k
  mpz_t n;
86
2.13k
  mp_limb_t D; /* The absolute value is stored. */
87
2.13k
  mp_limb_t g;
88
2.13k
  long Q;
89
2.13k
  mpz_t T1, T2;
90
91
  /* Test on the absolute value. */
92
2.13k
  mpz_roinit_n (n, PTR (x), ABSIZ (x));
93
94
2.13k
  ASSERT (mpz_odd_p (n));
95
  /* ASSERT (mpz_gcd_ui (NULL, n, 6) == 1); */
96
2.13k
#if GMP_NUMB_BITS % 16 == 0
97
  /* (2^12 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1) */
98
2.13k
  g = mpn_mod_34lsub1 (PTR (n), SIZ (n));
99
  /* (2^12 - 1) = 3^2 * 5 * 7 * 13    */
100
2.13k
  ASSERT (g % 3 != 0 && g % 5 != 0 && g % 7 != 0);
101
2.13k
  if ((g % 5 & 2) != 0)
102
    /* (5/n) = -1, iff n = 2 or 3 (mod 5) */
103
    /* D = 5; Q = -1 */
104
1.04k
    return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
105
1.09k
  else if (! POW2_P (g % 7))
106
    /* (-7/n) = -1, iff n = 3,5 or 6 (mod 7)  */
107
567
    D = 7; /* Q = 2 */
108
    /* (9/n) = -1, never: 9 = 3^2 */
109
524
  else if (mpz_oddjacobi_ui (n, 11) == -1)
110
    /* (-11/n) = (n/11) */
111
278
    D = 11; /* Q = 3 */
112
246
  else if ((((g % 13 - (g % 13 >> 3)) & 7) > 4) ||
113
246
     (((g % 13 - (g % 13 >> 3)) & 7) == 2))
114
    /* (13/n) = -1, iff n = 2,5,6,7,8 or 11 (mod 13)  */
115
116
    D = 13; /* Q = -3 */
116
130
  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
63
    D = 15; /* Q = 4 */
121
67
#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
67
  else if (! POW2_P (g % 17) && ! POW2_P (17 - g % 17))
125
    /* (17/n) = -1, iff n != +-1,+-2,+-4,+-8 (mod 17) */
126
30
    D = 17; /* Q = -4 */
127
37
#endif
128
#else
129
  if (mpz_oddjacobi_ui (n, 5) == -1)
130
    return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
131
#endif
132
37
  else
133
37
  {
134
37
    mp_limb_t maxD;
135
37
    int jac;
136
137
    /* n is odd, to possibly be a square, n % 8 = 1 is needed. */
138
37
    if (((*PTR (n) & 6) == 0) && UNLIKELY (mpz_perfect_square_p (n)))
139
0
      return 0; /* A square is composite. */
140
141
    /* Check Ds up to square root (in case, n is prime)
142
       or avoid overflows */
143
37
    if (SIZ (n) == 1)
144
5
      maxD = limb_apprsqrt (* PTR (n));
145
32
    else if (BITS_PER_ULONG >= GMP_NUMB_BITS && SIZ (n) == 2)
146
9
      mpn_sqrtrem (&maxD, (mp_ptr) NULL, PTR (n), 2);
147
23
    else
148
23
      maxD = GMP_NUMB_MAX;
149
37
    maxD = MIN (maxD, ULONG_MAX);
150
151
37
    unsigned Ddiff = 2;
152
37
#if GMP_NUMB_BITS % 16 == 0
153
37
    const unsigned D2 = 6;
154
37
#if GMP_NUMB_BITS % 32 == 0
155
37
    D = 19;
156
37
    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
37
    for (;;)
170
92
      {
171
92
  jac = mpz_oddjacobi_ui (n, D);
172
92
  if (jac != 1)
173
37
    break;
174
55
  if (UNLIKELY (D >= maxD))
175
0
    return 1;
176
55
  D += Ddiff;
177
55
  Ddiff = D2 - Ddiff;
178
55
      }
179
180
37
    if (UNLIKELY (jac == 0))
181
0
      return 0;
182
37
  }
183
184
  /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
185
1.09k
  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
1.09k
  b0 = mpz_scan0 (n, 0);
190
191
1.09k
  mpz_init (T1);
192
1.09k
  mpz_init (T2);
193
194
  /* If Ud != 0 && Vd != 0 */
195
1.09k
  if (mpz_lucas_mod (V, Qk, Q, b0, n, T1, T2) == 0)
196
398
    if (LIKELY (--b0 != 0))
197
398
      for (;;)
198
645
  {
199
    /* V_{2k} <- V_k ^ 2 - 2Q^k */
200
645
    mpz_mul (T2, V, V);
201
645
    mpz_submul_ui (T2, Qk, 2);
202
645
    mpz_tdiv_r (V, T2, n);
203
645
    if (SIZ (V) == 0 || UNLIKELY (--b0 == 0))
204
398
      break;
205
    /* Q^{2k} = (Q^k)^2 */
206
247
    mpz_mul (T2, Qk, Qk);
207
247
    mpz_tdiv_r (Qk, T2, n);
208
247
  }
209
210
1.09k
  mpz_clear (T1);
211
1.09k
  mpz_clear (T2);
212
213
1.09k
  return (b0 != 0);
214
2.13k
}