2 * Copyright (c) 2014 Colin Percival
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
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.
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
26 #include <sys/cdefs.h>
27 __FBSDID("$FreeBSD$");
35 /* Return a * b % n, where 0 < n. */
37 mulmod(uint64_t a, uint64_t b, uint64_t n)
45 if ((x < an) || (x >= n))
50 else if (an + an >= n)
60 /* Return a^r % n, where 0 < n. */
62 powmod(uint64_t a, uint64_t r, uint64_t n)
76 /* Return non-zero if n is a strong pseudoprime to base p. */
78 spsp(uint64_t n, uint64_t p)
84 /* Compute n - 1 = 2^k * r. */
85 while ((r & 1) == 0) {
90 /* Compute x = p^r mod n. If x = 1, n is a p-spsp. */
95 /* Compute x^(2^i) for 0 <= i < n. If any are -1, n is a p-spsp. */
107 /* Test for primality using strong pseudoprime tests. */
115 * C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr.,
116 * The pseudoprimes to 25 * 10^9, Math. Comp. 35(151):1003-1026, 1980.
119 /* No SPSPs to base 2 less than 2047. */
125 /* No SPSPs to bases 2,3 less than 1373653. */
131 /* No SPSPs to bases 2,3,5 less than 25326001. */
137 /* No SPSPs to bases 2,3,5,7 less than 3215031751. */
140 if (n < 3215031751ULL)
145 * G. Jaeschke, On strong pseudoprimes to several bases,
146 * Math. Comp. 61(204):915-926, 1993.
149 /* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */
152 if (n < 2152302898747ULL)
155 /* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */
158 if (n < 3474749660383ULL)
161 /* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */
164 if (n < 341550071728321ULL)
167 /* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */
170 if (n < 341550071728321ULL)
175 * Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime
176 * bases, Math. Comp. 83(290):2915-2924, 2014.
179 /* No SPSPs to bases 2..23 less than 3825123056546413051. */
182 if (n < 3825123056546413051)
187 * J. Sorenson and J. Webster, Strong pseudoprimes to twelve prime
188 * bases, Math. Comp. 86(304):985-1003, 2017.
191 /* No SPSPs to bases 2..37 less than 318665857834031151167461. */
199 /* All 64-bit values are less than 318665857834031151167461. */