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 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 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 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 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