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