1102782Skan/*- 2102782Skan * Copyright (c) 2014 Colin Percival 3169691Skan * All rights reserved. 4102782Skan * 5102782Skan * Redistribution and use in source and binary forms, with or without 6102782Skan * modification, are permitted provided that the following conditions 7102782Skan * are met: 8102782Skan * 1. Redistributions of source code must retain the above copyright 9102782Skan * notice, this list of conditions and the following disclaimer. 10102782Skan * 2. Redistributions in binary form must reproduce the above copyright 11102782Skan * notice, this list of conditions and the following disclaimer in the 12102782Skan * documentation and/or other materials provided with the distribution. 13102782Skan * 14102782Skan * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15102782Skan * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16102782Skan * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17102782Skan * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18102782Skan * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19169691Skan * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20102782Skan * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21102782Skan * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22102782Skan * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23102782Skan * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24102782Skan * SUCH DAMAGE. 25102782Skan */ 26102782Skan#include <sys/cdefs.h> 27102782Skan__FBSDID("$FreeBSD: releng/11.0/usr.bin/primes/spsp.c 272166 2014-09-26 09:40:48Z cperciva $"); 28102782Skan 29102782Skan#include <assert.h> 30102782Skan#include <stddef.h> 31102782Skan#include <stdint.h> 32102782Skan 33102782Skan#include "primes.h" 34102782Skan 35102782Skan/* Return a * b % n, where 0 <= a, b < 2^63, 0 < n < 2^63. */ 36102782Skanstatic uint64_t 37102782Skanmulmod(uint64_t a, uint64_t b, uint64_t n) 38169691Skan{ 39169691Skan uint64_t x = 0; 40102782Skan 41102782Skan while (b != 0) { 42102782Skan if (b & 1) 43102782Skan x = (x + a) % n; 44102782Skan a = (a + a) % n; 45102782Skan b >>= 1; 46102782Skan } 47102782Skan 48102782Skan return (x); 49102782Skan} 50102782Skan 51102782Skan/* Return a^r % n, where 0 <= a < 2^63, 0 < n < 2^63. */ 52102782Skanstatic uint64_t 53102782Skanpowmod(uint64_t a, uint64_t r, uint64_t n) 54102782Skan{ 55102782Skan uint64_t x = 1; 56102782Skan 57102782Skan while (r != 0) { 58169691Skan if (r & 1) 59102782Skan x = mulmod(a, x, n); 60169691Skan a = mulmod(a, a, n); 61169691Skan r >>= 1; 62169691Skan } 63169691Skan 64169691Skan return (x); 65169691Skan} 66169691Skan 67169691Skan/* Return non-zero if n is a strong pseudoprime to base p. */ 68169691Skanstatic int 69169691Skanspsp(uint64_t n, uint64_t p) 70169691Skan{ 71132720Skan uint64_t x; 72102782Skan uint64_t r = n - 1; 73102782Skan int k = 0; 74102782Skan 75102782Skan /* Compute n - 1 = 2^k * r. */ 76102782Skan while ((r & 1) == 0) { 77102782Skan k++; 78102782Skan r >>= 1; 79102782Skan } 80102782Skan 81102782Skan /* Compute x = p^r mod n. If x = 1, n is a p-spsp. */ 82102782Skan x = powmod(p, r, n); 83102782Skan if (x == 1) 84102782Skan return (1); 85102782Skan 86102782Skan /* Compute x^(2^i) for 0 <= i < n. If any are -1, n is a p-spsp. */ 87169691Skan while (k > 0) { 88169691Skan if (x == n - 1) 89169691Skan return (1); 90169691Skan x = powmod(x, 2, n); 91169691Skan k--; 92169691Skan } 93169691Skan 94169691Skan /* Not a p-spsp. */ 95102782Skan return (0); 96169691Skan} 97169691Skan 98169691Skan/* Test for primality using strong pseudoprime tests. */ 99169691Skanint 100169691Skanisprime(ubig _n) 101169691Skan{ 102169691Skan uint64_t n = _n; 103169691Skan 104169691Skan /* 105169691Skan * Values from: 106169691Skan * C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr., 107169691Skan * The pseudoprimes to 25 * 10^9, Math. Comp. 35(151):1003-1026, 1980. 108169691Skan */ 109169691Skan 110169691Skan /* No SPSPs to base 2 less than 2047. */ 111169691Skan if (!spsp(n, 2)) 112169691Skan return (0); 113169691Skan if (n < 2047ULL) 114169691Skan return (1); 115169691Skan 116169691Skan /* No SPSPs to bases 2,3 less than 1373653. */ 117 if (!spsp(n, 3)) 118 return (0); 119 if (n < 1373653ULL) 120 return (1); 121 122 /* No SPSPs to bases 2,3,5 less than 25326001. */ 123 if (!spsp(n, 5)) 124 return (0); 125 if (n < 25326001ULL) 126 return (1); 127 128 /* No SPSPs to bases 2,3,5,7 less than 3215031751. */ 129 if (!spsp(n, 7)) 130 return (0); 131 if (n < 3215031751ULL) 132 return (1); 133 134 /* 135 * Values from: 136 * G. Jaeschke, On strong pseudoprimes to several bases, 137 * Math. Comp. 61(204):915-926, 1993. 138 */ 139 140 /* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */ 141 if (!spsp(n, 11)) 142 return (0); 143 if (n < 2152302898747ULL) 144 return (1); 145 146 /* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */ 147 if (!spsp(n, 13)) 148 return (0); 149 if (n < 3474749660383ULL) 150 return (1); 151 152 /* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */ 153 if (!spsp(n, 17)) 154 return (0); 155 if (n < 341550071728321ULL) 156 return (1); 157 158 /* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */ 159 if (!spsp(n, 19)) 160 return (0); 161 if (n < 341550071728321ULL) 162 return (1); 163 164 /* 165 * Value from: 166 * Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime 167 * bases, Math. Comp. 83(290):2915-2924, 2014. 168 */ 169 170 /* No SPSPs to bases 2..23 less than 3825123056546413051. */ 171 if (!spsp(n, 23)) 172 return (0); 173 if (n < 3825123056546413051) 174 return (1); 175 176 /* We can't handle values larger than this. */ 177 assert(n <= SPSPMAX); 178 179 /* UNREACHABLE */ 180 return (0); 181} 182