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