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