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 <assert.h> 30 #include <stddef.h> 31 #include <stdint.h> 32 33 #include "primes.h" 34 35 /* Return a * b % n, where 0 < n. */ 36 static uint64_t 37 mulmod(uint64_t a, uint64_t b, uint64_t n) 38 { 39 uint64_t x = 0; 40 uint64_t an = a % n; 41 42 while (b != 0) { 43 if (b & 1) { 44 x += an; 45 if ((x < an) || (x >= n)) 46 x -= n; 47 } 48 if (an + an < an) 49 an = an + an - n; 50 else if (an + an >= n) 51 an = an + an - n; 52 else 53 an = an + an; 54 b >>= 1; 55 } 56 57 return (x); 58 } 59 60 /* Return a^r % n, where 0 < n. */ 61 static uint64_t 62 powmod(uint64_t a, uint64_t r, uint64_t n) 63 { 64 uint64_t x = 1; 65 66 while (r != 0) { 67 if (r & 1) 68 x = mulmod(a, x, n); 69 a = mulmod(a, a, n); 70 r >>= 1; 71 } 72 73 return (x); 74 } 75 76 /* Return non-zero if n is a strong pseudoprime to base p. */ 77 static int 78 spsp(uint64_t n, uint64_t p) 79 { 80 uint64_t x; 81 uint64_t r = n - 1; 82 int k = 0; 83 84 /* Compute n - 1 = 2^k * r. */ 85 while ((r & 1) == 0) { 86 k++; 87 r >>= 1; 88 } 89 90 /* Compute x = p^r mod n. If x = 1, n is a p-spsp. */ 91 x = powmod(p, r, n); 92 if (x == 1) 93 return (1); 94 95 /* Compute x^(2^i) for 0 <= i < n. If any are -1, n is a p-spsp. */ 96 while (k > 0) { 97 if (x == n - 1) 98 return (1); 99 x = powmod(x, 2, n); 100 k--; 101 } 102 103 /* Not a p-spsp. */ 104 return (0); 105 } 106 107 /* Test for primality using strong pseudoprime tests. */ 108 int 109 isprime(ubig _n) 110 { 111 uint64_t n = _n; 112 113 /* 114 * Values from: 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. 117 */ 118 119 /* No SPSPs to base 2 less than 2047. */ 120 if (!spsp(n, 2)) 121 return (0); 122 if (n < 2047ULL) 123 return (1); 124 125 /* No SPSPs to bases 2,3 less than 1373653. */ 126 if (!spsp(n, 3)) 127 return (0); 128 if (n < 1373653ULL) 129 return (1); 130 131 /* No SPSPs to bases 2,3,5 less than 25326001. */ 132 if (!spsp(n, 5)) 133 return (0); 134 if (n < 25326001ULL) 135 return (1); 136 137 /* No SPSPs to bases 2,3,5,7 less than 3215031751. */ 138 if (!spsp(n, 7)) 139 return (0); 140 if (n < 3215031751ULL) 141 return (1); 142 143 /* 144 * Values from: 145 * G. Jaeschke, On strong pseudoprimes to several bases, 146 * Math. Comp. 61(204):915-926, 1993. 147 */ 148 149 /* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */ 150 if (!spsp(n, 11)) 151 return (0); 152 if (n < 2152302898747ULL) 153 return (1); 154 155 /* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */ 156 if (!spsp(n, 13)) 157 return (0); 158 if (n < 3474749660383ULL) 159 return (1); 160 161 /* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */ 162 if (!spsp(n, 17)) 163 return (0); 164 if (n < 341550071728321ULL) 165 return (1); 166 167 /* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */ 168 if (!spsp(n, 19)) 169 return (0); 170 if (n < 341550071728321ULL) 171 return (1); 172 173 /* 174 * Value from: 175 * Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime 176 * bases, Math. Comp. 83(290):2915-2924, 2014. 177 */ 178 179 /* No SPSPs to bases 2..23 less than 3825123056546413051. */ 180 if (!spsp(n, 23)) 181 return (0); 182 if (n < 3825123056546413051) 183 return (1); 184 185 /* 186 * Value from: 187 * J. Sorenson and J. Webster, Strong pseudoprimes to twelve prime 188 * bases, Math. Comp. 86(304):985-1003, 2017. 189 */ 190 191 /* No SPSPs to bases 2..37 less than 318665857834031151167461. */ 192 if (!spsp(n, 29)) 193 return (0); 194 if (!spsp(n, 31)) 195 return (0); 196 if (!spsp(n, 37)) 197 return (0); 198 199 /* All 64-bit values are less than 318665857834031151167461. */ 200 return (1); 201 } 202