1 /*-
2 * Copyright (c) 2014 Colin Percival
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 */
26 #include <sys/cdefs.h>
27 __FBSDID("$FreeBSD$");
28
29 #include <stddef.h>
30 #include <stdint.h>
31
32 #include "primes.h"
33
34 /* Return a * b % n, where 0 < n. */
35 static uint64_t
mulmod(uint64_t a,uint64_t b,uint64_t n)36 mulmod(uint64_t a, uint64_t b, uint64_t n)
37 {
38 uint64_t x = 0;
39 uint64_t an = a % n;
40
41 while (b != 0) {
42 if (b & 1) {
43 x += an;
44 if ((x < an) || (x >= n))
45 x -= n;
46 }
47 if (an + an < an)
48 an = an + an - n;
49 else if (an + an >= n)
50 an = an + an - n;
51 else
52 an = an + an;
53 b >>= 1;
54 }
55
56 return (x);
57 }
58
59 /* Return a^r % n, where 0 < n. */
60 static uint64_t
powmod(uint64_t a,uint64_t r,uint64_t n)61 powmod(uint64_t a, uint64_t r, uint64_t n)
62 {
63 uint64_t x = 1;
64
65 while (r != 0) {
66 if (r & 1)
67 x = mulmod(a, x, n);
68 a = mulmod(a, a, n);
69 r >>= 1;
70 }
71
72 return (x);
73 }
74
75 /* Return non-zero if n is a strong pseudoprime to base p. */
76 static int
spsp(uint64_t n,uint64_t p)77 spsp(uint64_t n, uint64_t p)
78 {
79 uint64_t x;
80 uint64_t r = n - 1;
81 int k = 0;
82
83 /* Compute n - 1 = 2^k * r. */
84 while ((r & 1) == 0) {
85 k++;
86 r >>= 1;
87 }
88
89 /* Compute x = p^r mod n. If x = 1, n is a p-spsp. */
90 x = powmod(p, r, n);
91 if (x == 1)
92 return (1);
93
94 /* Compute x^(2^i) for 0 <= i < n. If any are -1, n is a p-spsp. */
95 while (k > 0) {
96 if (x == n - 1)
97 return (1);
98 x = powmod(x, 2, n);
99 k--;
100 }
101
102 /* Not a p-spsp. */
103 return (0);
104 }
105
106 /* Test for primality using strong pseudoprime tests. */
107 int
isprime(ubig _n)108 isprime(ubig _n)
109 {
110 uint64_t n = _n;
111
112 /*
113 * Values from:
114 * C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr.,
115 * The pseudoprimes to 25 * 10^9, Math. Comp. 35(151):1003-1026, 1980.
116 */
117
118 /* No SPSPs to base 2 less than 2047. */
119 if (!spsp(n, 2))
120 return (0);
121 if (n < 2047ULL)
122 return (1);
123
124 /* No SPSPs to bases 2,3 less than 1373653. */
125 if (!spsp(n, 3))
126 return (0);
127 if (n < 1373653ULL)
128 return (1);
129
130 /* No SPSPs to bases 2,3,5 less than 25326001. */
131 if (!spsp(n, 5))
132 return (0);
133 if (n < 25326001ULL)
134 return (1);
135
136 /* No SPSPs to bases 2,3,5,7 less than 3215031751. */
137 if (!spsp(n, 7))
138 return (0);
139 if (n < 3215031751ULL)
140 return (1);
141
142 /*
143 * Values from:
144 * G. Jaeschke, On strong pseudoprimes to several bases,
145 * Math. Comp. 61(204):915-926, 1993.
146 */
147
148 /* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */
149 if (!spsp(n, 11))
150 return (0);
151 if (n < 2152302898747ULL)
152 return (1);
153
154 /* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */
155 if (!spsp(n, 13))
156 return (0);
157 if (n < 3474749660383ULL)
158 return (1);
159
160 /* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */
161 if (!spsp(n, 17))
162 return (0);
163 if (n < 341550071728321ULL)
164 return (1);
165
166 /* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */
167 if (!spsp(n, 19))
168 return (0);
169 if (n < 341550071728321ULL)
170 return (1);
171
172 /*
173 * Value from:
174 * Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime
175 * bases, Math. Comp. 83(290):2915-2924, 2014.
176 */
177
178 /* No SPSPs to bases 2..23 less than 3825123056546413051. */
179 if (!spsp(n, 23))
180 return (0);
181 if (n < 3825123056546413051)
182 return (1);
183
184 /*
185 * Value from:
186 * J. Sorenson and J. Webster, Strong pseudoprimes to twelve prime
187 * bases, Math. Comp. 86(304):985-1003, 2017.
188 */
189
190 /* No SPSPs to bases 2..37 less than 318665857834031151167461. */
191 if (!spsp(n, 29))
192 return (0);
193 if (!spsp(n, 31))
194 return (0);
195 if (!spsp(n, 37))
196 return (0);
197
198 /* All 64-bit values are less than 318665857834031151167461. */
199 return (1);
200 }
201