Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpz/millerrabin.c
Line
Count
Source (jump to first uncovered line)
1
/* mpz_millerrabin(n,reps) -- An implementation of the probabilistic primality
2
   test found in Knuth's Seminumerical Algorithms book.  If the function
3
   mpz_millerrabin() returns 0 then n is not prime.  If it returns 1, then n is
4
   'probably' prime.  The probability of a false positive is (1/4)**reps, where
5
   reps is the number of internal passes of the probabilistic algorithm.  Knuth
6
   indicates that 25 passes are reasonable.
7
8
   With the current implementation, the first 24 MR-tests are substituted by a
9
   Baillie-PSW probable prime test.
10
11
   This implementation of the Baillie-PSW test was checked up to 2640*10^12,
12
   for smaller values no MR-test is performed, regardless of reps, and
13
   2 ("surely prime") is returned if the number was not proved composite.
14
15
   If GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS is defined as non-zero,
16
   the code assumes that the Baillie-PSW test was checked up to 2^64.
17
18
   THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
19
   CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
20
   FUTURE GNU MP RELEASES.
21
22
Copyright 1991, 1993, 1994, 1996-2002, 2005, 2014, 2018-2022, 2024 Free
23
Software Foundation, Inc.
24
25
Contributed by John Amanatides.
26
Changed to "BPSW, then Miller Rabin if required" by Marco Bodrato.
27
28
This file is part of the GNU MP Library.
29
30
The GNU MP Library is free software; you can redistribute it and/or modify
31
it under the terms of either:
32
33
  * the GNU Lesser General Public License as published by the Free
34
    Software Foundation; either version 3 of the License, or (at your
35
    option) any later version.
36
37
or
38
39
  * the GNU General Public License as published by the Free Software
40
    Foundation; either version 2 of the License, or (at your option) any
41
    later version.
42
43
or both in parallel, as here.
44
45
The GNU MP Library is distributed in the hope that it will be useful, but
46
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
47
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
48
for more details.
49
50
You should have received copies of the GNU General Public License and the
51
GNU Lesser General Public License along with the GNU MP Library.  If not,
52
see https://www.gnu.org/licenses/.  */
53
54
#include "gmp-impl.h"
55
56
#ifndef GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
57
#define GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS 0
58
#endif
59
60
static int
61
mod_eq_m1 (mpz_srcptr x, mpz_srcptr m)
62
46.7k
{
63
46.7k
  mp_size_t ms;
64
46.7k
  mp_srcptr mp, xp;
65
66
46.7k
  ms = SIZ (m);
67
46.7k
  if (SIZ (x) != ms)
68
6.06k
    return 0;
69
40.6k
  ASSERT (ms > 0);
70
71
40.6k
  mp = PTR (m);
72
40.6k
  xp = PTR (x);
73
40.6k
  ASSERT ((mp[0] - 1) == (mp[0] ^ 1)); /* n is odd */
74
75
40.6k
  if ((*xp ^ CNST_LIMB(1) ^ *mp) != CNST_LIMB(0)) /* xp[0] != mp[0] - 1 */
76
37.9k
    return 0;
77
2.67k
  else
78
2.67k
    {
79
2.67k
      int cmp;
80
81
2.67k
      --ms;
82
2.67k
      ++xp;
83
2.67k
      ++mp;
84
85
2.67k
      MPN_CMP (cmp, xp, mp, ms);
86
87
2.67k
      return cmp == 0;
88
2.67k
    }
89
40.6k
}
90
91
static int
92
millerrabin (mpz_srcptr n, mpz_ptr x, mpz_ptr y,
93
       mpz_srcptr q, mp_bitcnt_t k)
94
24.4k
{
95
24.4k
  mpz_powm (y, x, q, n);
96
97
24.4k
  if (mpz_cmp_ui (y, 1L) == 0 || mod_eq_m1 (y, n))
98
2.41k
    return 1;
99
100
44.0k
  for (mp_bitcnt_t i = 1; i < k; ++i)
101
23.4k
    {
102
23.4k
      mpz_powm_ui (y, y, 2L, n);
103
23.4k
      if (mod_eq_m1 (y, n))
104
1.50k
  return 1;
105
23.4k
    }
106
20.5k
  return 0;
107
22.0k
}
108
109
int
110
mpz_millerrabin (mpz_srcptr n, int reps)
111
22.7k
{
112
22.7k
  mpz_t nm, x, y, q;
113
22.7k
  mp_bitcnt_t k;
114
22.7k
  int is_prime;
115
22.7k
  TMP_DECL;
116
22.7k
  TMP_MARK;
117
118
22.7k
  ASSERT (SIZ (n) > 0);
119
22.7k
  MPZ_TMP_INIT (nm, SIZ (n) + 1);
120
22.7k
  mpz_tdiv_q_2exp (nm, n, 1);
121
122
22.7k
  MPZ_TMP_INIT (x, SIZ (n) + 1);
123
22.7k
  MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
124
22.7k
  MPZ_TMP_INIT (q, SIZ (n));
125
126
  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
127
22.7k
  k = mpn_scan1 (PTR (nm), 0);
128
22.7k
  mpz_tdiv_q_2exp (q, nm, k);
129
22.7k
  ++k;
130
131
  /* BPSW test */
132
22.7k
  mpz_set_ui (x, 2);
133
22.7k
  is_prime = millerrabin (n, x, y, q, k) && mpz_stronglucas (n, x, y);
134
135
22.7k
  if (is_prime)
136
2.13k
    {
137
2.13k
      if (
138
#if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
139
    /* Consider numbers up to 2^64 that pass the BPSW test as primes. */
140
#if GMP_NUMB_BITS <= 64
141
    SIZ (n) <= 64 / GMP_NUMB_BITS
142
#else
143
    0
144
#endif
145
#if 64 % GMP_NUMB_BITS != 0
146
    || SIZ (n) - 64 / GMP_NUMB_BITS == (PTR (n) [64 / GMP_NUMB_BITS] < CNST_LIMB(1) << 64 % GMP_NUMB_BITS)
147
#endif
148
#else
149
    /* Consider numbers up to 75*2^45 that pass the BPSW test as primes.
150
       This implementation was tested up to 264*10^13 > 2^51+2^48+2^46+2^45 */
151
    /* 2^6 < 75 = 0b1001011 < 2^7 */
152
2.13k
#define GMP_BPSW_LIMB_CONST CNST_LIMB(75)
153
8.54k
#define GMP_BPSW_BITS_CONST (LOG2C(75) - 1)
154
6.40k
#define GMP_BPSW_BITS_LIMIT (45 + GMP_BPSW_BITS_CONST)
155
156
4.27k
#define GMP_BPSW_LIMBS_LIMIT (GMP_BPSW_BITS_LIMIT / GMP_NUMB_BITS)
157
2.13k
#define GMP_BPSW_BITS_MOD (GMP_BPSW_BITS_LIMIT % GMP_NUMB_BITS)
158
159
#if GMP_NUMB_BITS <=  GMP_BPSW_BITS_LIMIT
160
    SIZ (n) <= GMP_BPSW_LIMBS_LIMIT
161
#else
162
2.13k
    0
163
2.13k
#endif
164
2.13k
#if GMP_BPSW_BITS_MOD >=  GMP_BPSW_BITS_CONST
165
2.13k
    || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST << (GMP_BPSW_BITS_MOD - GMP_BPSW_BITS_CONST))
166
#else
167
#if GMP_BPSW_BITS_MOD != 0
168
    || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST >> (GMP_BPSW_BITS_CONST -  GMP_BPSW_BITS_MOD))
169
#else
170
#if GMP_NUMB_BITS > GMP_BPSW_BITS_CONST
171
    || SIZ (nm) - GMP_BPSW_LIMBS_LIMIT + 1 == (PTR (nm) [GMP_BPSW_LIMBS_LIMIT - 1] < GMP_BPSW_LIMB_CONST << (GMP_NUMB_BITS - 1 - GMP_BPSW_BITS_CONST))
172
#endif
173
#endif
174
#endif
175
176
2.13k
#undef GMP_BPSW_BITS_LIMIT
177
2.13k
#undef GMP_BPSW_LIMB_CONST
178
2.13k
#undef GMP_BPSW_BITS_CONST
179
2.13k
#undef GMP_BPSW_LIMBS_LIMIT
180
2.13k
#undef GMP_BPSW_BITS_MOD
181
182
2.13k
#endif
183
2.13k
    )
184
259
  is_prime = 2;
185
1.87k
      else
186
1.87k
  {
187
1.87k
    reps -= 24;
188
1.87k
    if (reps > 0)
189
1.78k
      {
190
1.78k
        gmp_randstate_t rstate;
191
        /* (n-5)/2 */
192
1.78k
        mpz_sub_ui (nm, nm, 2L);
193
1.78k
        ASSERT (mpz_cmp_ui (nm, 1L) >= 0);
194
195
1.78k
        gmp_randinit_default (rstate);
196
197
1.78k
        do
198
1.78k
    {
199
      /* 3 to (n-1)/2 inclusive, don't want 1, 0 or 2 */
200
1.78k
      mpz_urandomm (x, rstate, nm);
201
1.78k
      mpz_add_ui (x, x, 3L);
202
203
1.78k
      is_prime = millerrabin (n, x, y, q, k);
204
1.78k
    } while (--reps > 0 && is_prime);
205
206
1.78k
        gmp_randclear (rstate);
207
1.78k
      }
208
1.87k
  }
209
2.13k
    }
210
22.7k
  TMP_FREE;
211
22.7k
  return is_prime;
212
22.7k
}