Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpz/pprime_p.c
Line
Count
Source (jump to first uncovered line)
1
/* mpz_probab_prime_p --
2
   An implementation of the probabilistic primality test found in Knuth's
3
   Seminumerical Algorithms book.  If the function mpz_probab_prime_p()
4
   returns 0 then n is not prime.  If it returns 1, then n is 'probably'
5
   prime.  If it returns 2, n is surely prime.  The probability of a false
6
   positive is (1/4)**reps, where reps is the number of internal passes of the
7
   probabilistic algorithm.  Knuth indicates that 25 passes are reasonable.
8
9
Copyright 1991, 1993, 1994, 1996-2002, 2005, 2015, 2016 Free Software
10
Foundation, Inc.
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
static int isprime (unsigned long int);
42
43
44
/* MPN_MOD_OR_MODEXACT_1_ODD can be used instead of mpn_mod_1 for the trial
45
   division.  It gives a result which is not the actual remainder r but a
46
   value congruent to r*2^n mod d.  Since all the primes being tested are
47
   odd, r*2^n mod p will be 0 if and only if r mod p is 0.  */
48
49
int
50
mpz_probab_prime_p (mpz_srcptr n, int reps)
51
658
{
52
658
  mp_limb_t r;
53
658
  mpz_t n2;
54
55
  /* Handle small and negative n.  */
56
658
  if (mpz_cmp_ui (n, 1000000L) <= 0)
57
40
    {
58
40
      if (mpz_cmpabs_ui (n, 1000000L) <= 0)
59
40
  {
60
40
    int is_prime;
61
40
    unsigned long n0;
62
40
    n0 = mpz_get_ui (n);
63
40
    is_prime = n0 & (n0 > 1) ? isprime (n0) : n0 == 2;
64
40
    return is_prime ? 2 : 0;
65
40
  }
66
      /* Negative number.  Negate and fall out.  */
67
0
      PTR(n2) = PTR(n);
68
0
      SIZ(n2) = -SIZ(n);
69
0
      n = n2;
70
0
    }
71
72
  /* If n is now even, it is not a prime.  */
73
618
  if (mpz_even_p (n))
74
31
    return 0;
75
76
587
#if defined (PP)
77
  /* Check if n has small factors.  */
78
587
#if defined (PP_INVERTED)
79
587
  r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP,
80
587
             (mp_limb_t) PP_INVERTED);
81
#else
82
  r = mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP);
83
#endif
84
587
  if (r % 3 == 0
85
587
#if GMP_LIMB_BITS >= 4
86
587
      || r % 5 == 0
87
587
#endif
88
587
#if GMP_LIMB_BITS >= 8
89
587
      || r % 7 == 0
90
587
#endif
91
587
#if GMP_LIMB_BITS >= 16
92
587
      || r % 11 == 0 || r % 13 == 0
93
587
#endif
94
587
#if GMP_LIMB_BITS >= 32
95
587
      || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
96
587
#endif
97
587
#if GMP_LIMB_BITS >= 64
98
587
      || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
99
587
      || r % 47 == 0 || r % 53 == 0
100
587
#endif
101
587
      )
102
205
    {
103
205
      return 0;
104
205
    }
105
382
#endif /* PP */
106
107
  /* Do more dividing.  We collect small primes, using umul_ppmm, until we
108
     overflow a single limb.  We divide our number by the small primes product,
109
     and look for factors in the remainder.  */
110
382
  {
111
382
    unsigned long int ln2;
112
382
    unsigned long int q;
113
382
    mp_limb_t p1, p0, p;
114
382
    unsigned int primes[15];
115
382
    int nprimes;
116
117
382
    nprimes = 0;
118
382
    p = 1;
119
382
    ln2 = mpz_sizeinbase (n, 2);  /* FIXME: tune this limit */
120
451k
    for (q = PP_FIRST_OMITTED; q < ln2; q += 2)
121
451k
      {
122
451k
  if (isprime (q))
123
117k
    {
124
117k
      umul_ppmm (p1, p0, p, q);
125
117k
      if (p1 != 0)
126
21.6k
        {
127
21.6k
    r = MPN_MOD_OR_MODEXACT_1_ODD (PTR(n), (mp_size_t) SIZ(n), p);
128
137k
    while (--nprimes >= 0)
129
116k
      if (r % primes[nprimes] == 0)
130
90
        {
131
90
          ASSERT_ALWAYS (mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) primes[nprimes]) == 0);
132
90
          return 0;
133
90
        }
134
21.5k
    p = q;
135
21.5k
    nprimes = 0;
136
21.5k
        }
137
96.0k
      else
138
96.0k
        {
139
96.0k
    p = p0;
140
96.0k
        }
141
117k
      primes[nprimes++] = q;
142
117k
    }
143
451k
      }
144
382
  }
145
146
  /* Perform a number of Miller-Rabin tests.  */
147
292
  return mpz_millerrabin (n, reps);
148
382
}
149
150
static int
151
isprime (unsigned long int t)
152
451k
{
153
451k
  unsigned long int q, r, d;
154
155
451k
  ASSERT (t >= 3 && (t & 1) != 0);
156
157
451k
  d = 3;
158
4.44M
  do {
159
4.44M
      q = t / d;
160
4.44M
      r = t - q * d;
161
4.44M
      if (q < d)
162
117k
  return 1;
163
4.33M
      d += 2;
164
4.33M
  } while (r != 0);
165
333k
  return 0;
166
451k
}