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