Coverage Report

Created: 2025-03-06 07:58

/src/gmp/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 2463*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 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 millerrabin (mpz_srcptr,
61
      mpz_ptr, mpz_ptr,
62
      mpz_srcptr, unsigned long int);
63
64
int
65
mpz_millerrabin (mpz_srcptr n, int reps)
66
0
{
67
0
  mpz_t nm, x, y, q;
68
0
  mp_bitcnt_t k;
69
0
  int is_prime;
70
0
  TMP_DECL;
71
0
  TMP_MARK;
72
73
0
  ASSERT (SIZ (n) > 0);
74
0
  MPZ_TMP_INIT (nm, SIZ (n) + 1);
75
0
  mpz_tdiv_q_2exp (nm, n, 1);
76
77
0
  MPZ_TMP_INIT (x, SIZ (n) + 1);
78
0
  MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
79
0
  MPZ_TMP_INIT (q, SIZ (n));
80
81
  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
82
0
  k = mpn_scan1 (PTR (nm), 0);
83
0
  mpz_tdiv_q_2exp (q, nm, k);
84
0
  ++k;
85
86
  /* BPSW test */
87
0
  mpz_set_ui (x, 2);
88
0
  is_prime = millerrabin (n, x, y, q, k) && mpz_stronglucas (n, x, y);
89
90
0
  if (is_prime)
91
0
    {
92
0
      if (
93
#if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
94
    /* Consider numbers up to 2^64 that pass the BPSW test as primes. */
95
#if GMP_NUMB_BITS <= 64
96
    SIZ (n) <= 64 / GMP_NUMB_BITS
97
#else
98
    0
99
#endif
100
#if 64 % GMP_NUMB_BITS != 0
101
    || SIZ (n) - 64 / GMP_NUMB_BITS == (PTR (n) [64 / GMP_NUMB_BITS] < CNST_LIMB(1) << 64 % GMP_NUMB_BITS)
102
#endif
103
#else
104
    /* Consider numbers up to 35*2^46 that pass the BPSW test as primes.
105
       This implementation was tested up to 2463*10^12 > 2^51+2^47+2^46 */
106
    /* 2^5 < 35 = 0b100011 < 2^6 */
107
0
#define GMP_BPSW_LIMB_CONST CNST_LIMB(35)
108
0
#define GMP_BPSW_BITS_CONST (LOG2C(35) - 1)
109
0
#define GMP_BPSW_BITS_LIMIT (46 + GMP_BPSW_BITS_CONST)
110
111
0
#define GMP_BPSW_LIMBS_LIMIT (GMP_BPSW_BITS_LIMIT / GMP_NUMB_BITS)
112
0
#define GMP_BPSW_BITS_MOD (GMP_BPSW_BITS_LIMIT % GMP_NUMB_BITS)
113
114
#if GMP_NUMB_BITS <=  GMP_BPSW_BITS_LIMIT
115
    SIZ (n) <= GMP_BPSW_LIMBS_LIMIT
116
#else
117
0
    0
118
0
#endif
119
0
#if GMP_BPSW_BITS_MOD >=  GMP_BPSW_BITS_CONST
120
0
    || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST << (GMP_BPSW_BITS_MOD - GMP_BPSW_BITS_CONST))
121
#else
122
#if GMP_BPSW_BITS_MOD != 0
123
    || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST >> (GMP_BPSW_BITS_CONST -  GMP_BPSW_BITS_MOD))
124
#else
125
#if GMP_NUMB_BITS > GMP_BPSW_BITS_CONST
126
    || 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))
127
#endif
128
#endif
129
#endif
130
131
0
#undef GMP_BPSW_BITS_LIMIT
132
0
#undef GMP_BPSW_LIMB_CONST
133
0
#undef GMP_BPSW_BITS_CONST
134
0
#undef GMP_BPSW_LIMBS_LIMIT
135
0
#undef GMP_BPSW_BITS_MOD
136
137
0
#endif
138
0
    )
139
0
  is_prime = 2;
140
0
      else
141
0
  {
142
0
    reps -= 24;
143
0
    if (reps > 0)
144
0
      {
145
0
        gmp_randstate_t rstate;
146
        /* (n-5)/2 */
147
0
        mpz_sub_ui (nm, nm, 2L);
148
0
        ASSERT (mpz_cmp_ui (nm, 1L) >= 0);
149
150
0
        gmp_randinit_default (rstate);
151
152
0
        do
153
0
    {
154
      /* 3 to (n-1)/2 inclusive, don't want 1, 0 or 2 */
155
0
      mpz_urandomm (x, rstate, nm);
156
0
      mpz_add_ui (x, x, 3L);
157
158
0
      is_prime = millerrabin (n, x, y, q, k);
159
0
    } while (--reps > 0 && is_prime);
160
161
0
        gmp_randclear (rstate);
162
0
      }
163
0
  }
164
0
    }
165
0
  TMP_FREE;
166
0
  return is_prime;
167
0
}
168
169
static int
170
mod_eq_m1 (mpz_srcptr x, mpz_srcptr m)
171
0
{
172
0
  mp_size_t ms;
173
0
  mp_srcptr mp, xp;
174
175
0
  ms = SIZ (m);
176
0
  if (SIZ (x) != ms)
177
0
    return 0;
178
0
  ASSERT (ms > 0);
179
180
0
  mp = PTR (m);
181
0
  xp = PTR (x);
182
0
  ASSERT ((mp[0] - 1) == (mp[0] ^ 1)); /* n is odd */
183
184
0
  if ((*xp ^ CNST_LIMB(1) ^ *mp) != CNST_LIMB(0)) /* xp[0] != mp[0] - 1 */
185
0
    return 0;
186
0
  else
187
0
    {
188
0
      int cmp;
189
190
0
      --ms;
191
0
      ++xp;
192
0
      ++mp;
193
194
0
      MPN_CMP (cmp, xp, mp, ms);
195
196
0
      return cmp == 0;
197
0
    }
198
0
}
199
200
static int
201
millerrabin (mpz_srcptr n, mpz_ptr x, mpz_ptr y,
202
       mpz_srcptr q, mp_bitcnt_t k)
203
0
{
204
0
  mpz_powm (y, x, q, n);
205
206
0
  if (mpz_cmp_ui (y, 1L) == 0 || mod_eq_m1 (y, n))
207
0
    return 1;
208
209
0
  for (mp_bitcnt_t i = 1; i < k; ++i)
210
0
    {
211
0
      mpz_powm_ui (y, y, 2L, n);
212
0
      if (mod_eq_m1 (y, n))
213
0
  return 1;
214
0
    }
215
0
  return 0;
216
0
}