1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.numbers.primes;
18
19 import java.math.BigInteger;
20 import java.util.AbstractMap.SimpleImmutableEntry;
21 import java.util.ArrayList;
22 import java.util.Arrays;
23 import java.util.HashSet;
24 import java.util.List;
25 import java.util.Map.Entry;
26 import java.util.Set;
27
28 /**
29 * Utility methods to work on primes within the <code>int</code> range.
30 */
31 final class SmallPrimes {
32 /**
33 * The first 512 prime numbers.
34 * <p>
35 * It contains all primes smaller or equal to the cubic square of Integer.MAX_VALUE.
36 * As a result, <code>int</code> numbers which are not reduced by those primes are guaranteed
37 * to be either prime or semi prime.
38 */
39 static final int[] PRIMES = {
40 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107,
41 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
42 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359,
43 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
44 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641,
45 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787,
46 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
47 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
48 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213,
49 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
50 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481,
51 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601,
52 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
53 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
54 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017,
55 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143,
56 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
57 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
58 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593,
59 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713,
60 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
61 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
62 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181,
63 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323,
64 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467,
65 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
66 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671};
67
68 /** The last number in {@link #PRIMES}. */
69 static final int PRIMES_LAST = PRIMES[PRIMES.length - 1];
70
71 /**
72 * A set of prime numbers mapped to an array of all integers between
73 * 0 (inclusive) and the least common multiple, i.e. the product, of those
74 * prime numbers (exclusive) that are not divisible by any of these prime
75 * numbers. The prime numbers in the set are among the first 512 prime
76 * numbers, and the {@code int} array containing the numbers undivisible by
77 * these prime numbers is sorted in ascending order.
78 *
79 * <p>The purpose of this field is to speed up trial division by skipping
80 * multiples of individual prime numbers, specifically those contained
81 * in the key of this {@code Entry}, by only trying integers that are equivalent
82 * to one of the integers contained in the value of this {@code Entry} modulo
83 * the least common multiple of the prime numbers in the set.</p>
84 *
85 * <p>Note that, if {@code product} is the product of the prime numbers,
86 * the last number in the array of coprime integers is necessarily
87 * {@code product - 1}, because if {@code product ≡ 0 mod p}, then
88 * {@code product - 1 ≡ -1 mod p}, and {@code 0 ≢ -1 mod p} for any prime number p.</p>
89 */
90 static final Entry<Set<Integer>, int[]> PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES;
91
92 static {
93 /*
94 According to the Chinese Remainder Theorem, for every combination of
95 congruence classes modulo distinct, pairwise coprime moduli, there
96 exists exactly one congruence class modulo the product of these
97 moduli that is contained in every one of the former congruence
98 classes. Since the number of congruence classes coprime to a prime
99 number p is p-1, the number of congruence classes coprime to all
100 prime numbers p_1, p_2, p_3 … is (p_1 - 1) * (p_2 - 1) * (p_3 - 1) …
101
102 Therefore, when using the first five prime numbers as those whose multiples
103 are to be skipped in trial division, the array containing the coprime
104 equivalence classes will have to hold (2-1)*(3-1)*(5-1)*(7-1)*(11-1) = 480
105 values. As a consequence, the amount of integers to be tried in
106 trial division is reduced to 480/(2*3*5*7*11), which is about 20.78%,
107 of all integers.
108 */
109 final Set<Integer> primeNumbers = new HashSet<>();
110 primeNumbers.add(Integer.valueOf(2));
111 primeNumbers.add(Integer.valueOf(3));
112 primeNumbers.add(Integer.valueOf(5));
113 primeNumbers.add(Integer.valueOf(7));
114 primeNumbers.add(Integer.valueOf(11));
115
116 final int product = primeNumbers.stream().reduce(1, (a, b) -> a * b);
117 final int[] equivalenceClasses = new int[primeNumbers.stream().mapToInt(a -> a - 1).reduce(1, (a, b) -> a * b)];
118
119 int equivalenceClassIndex = 0;
120 for (int i = 0; i < product; i++) {
121 boolean foundPrimeFactor = false;
122 for (final Integer prime : primeNumbers) {
123 if (i % prime == 0) {
124 foundPrimeFactor = true;
125 break;
126 }
127 }
128 if (!foundPrimeFactor) {
129 equivalenceClasses[equivalenceClassIndex] = i;
130 equivalenceClassIndex++;
131 }
132 }
133
134 PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES = new SimpleImmutableEntry<>(primeNumbers, equivalenceClasses);
135 }
136
137 /**
138 * Utility class.
139 */
140 private SmallPrimes() {}
141
142 /**
143 * Extract small factors.
144 *
145 * @param n Number to factor, must be > 0.
146 * @param factors List where to add the factors.
147 * @return the part of {@code n} which remains to be factored, it is either
148 * a prime or a semi-prime.
149 */
150 static int smallTrialDivision(int n,
151 final List<Integer> factors) {
152 for (final int p : PRIMES) {
153 while (0 == n % p) {
154 n /= p;
155 factors.add(p);
156 }
157 }
158 return n;
159 }
160
161 /**
162 * Extract factors between {@code PRIME_LAST + 2} and {@code maxFactors}.
163 *
164 * @param n Number to factorize, must be larger than {@code PRIME_LAST + 2}
165 * and must not contain any factor below {@code PRIME_LAST + 2}.
166 * @param maxFactor Upper bound of trial division: if it is reached, the
167 * method gives up and returns {@code n}.
168 * @param factors the list where to add the factors.
169 * @return {@code n} (or 1 if factorization is completed).
170 */
171 static int boundedTrialDivision(int n,
172 int maxFactor,
173 List<Integer> factors) {
174 final int minFactor = PRIMES_LAST + 2;
175
176 /*
177 only trying integers of the form k*m + c, where k >= 0, m is the
178 product of some prime numbers which n is required not to contain
179 as prime factors, and c is an integer undivisible by all of those
180 prime numbers; in other words, skipping multiples of these primes
181 */
182 final int m = PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue()
183 [PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue().length - 1] + 1;
184 int km = m * (minFactor / m);
185 int currentEquivalenceClassIndex = Arrays.binarySearch(
186 PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue(),
187 minFactor % m);
188
189 /*
190 Since minFactor is the next smallest prime number after the
191 first 512 primes, it cannot be a multiple of one of them, therefore,
192 the index returned by the above binary search must be non-negative.
193 */
194
195 boolean done = false;
196 while (!done) {
197 // no check is done about n >= f
198 final int f = km + PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue()[currentEquivalenceClassIndex];
199 if (f > maxFactor) {
200 done = true;
201 } else if (0 == n % f) {
202 n /= f;
203 factors.add(f);
204 done = true;
205 } else {
206 if (currentEquivalenceClassIndex ==
207 PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue().length - 1) {
208 km += m;
209 currentEquivalenceClassIndex = 0;
210 } else {
211 currentEquivalenceClassIndex++;
212 }
213 }
214 }
215 if (n != 1) {
216 factors.add(n);
217 }
218 return n;
219 }
220
221 /**
222 * Factorization by trial division.
223 *
224 * @param n Number to factor.
225 * @return the list of prime factors of {@code n}.
226 */
227 static List<Integer> trialDivision(int n) {
228 final List<Integer> factors = new ArrayList<>(32);
229 n = smallTrialDivision(n, factors);
230 if (1 == n) {
231 return factors;
232 }
233 // here we are sure that n is either a prime or a semi prime
234 final int bound = (int) Math.sqrt(n);
235 boundedTrialDivision(n, bound, factors);
236 return factors;
237 }
238
239 /**
240 * Miller-Rabin probabilistic primality test for int type, used in such
241 * a way that a result is always guaranteed.
242 * <p>
243 * It uses the prime numbers as successive base therefore it is guaranteed
244 * to be always correct (see Handbook of applied cryptography by Menezes,
245 * table 4.1).
246 *
247 * @param n Number to test: an odd integer ≥ 3.
248 * @return true if {@code n} is prime, false if it is definitely composite.
249 */
250 static boolean millerRabinPrimeTest(final int n) {
251 final int nMinus1 = n - 1;
252 final int s = Integer.numberOfTrailingZeros(nMinus1);
253 final int r = nMinus1 >> s;
254 // r must be odd, it is not checked here
255 int t = 1;
256 if (n >= 2047) {
257 t = 2;
258 }
259 if (n >= 1373653) {
260 t = 3;
261 }
262 if (n >= 25326001) {
263 t = 4;
264 } // works up to 3.2 billion, int range stops at 2.7 so we are safe :-)
265 final BigInteger br = BigInteger.valueOf(r);
266 final BigInteger bn = BigInteger.valueOf(n);
267
268 for (int i = 0; i < t; i++) {
269 final BigInteger a = BigInteger.valueOf(SmallPrimes.PRIMES[i]);
270 final BigInteger bPow = a.modPow(br, bn);
271 int y = bPow.intValue();
272 if (1 != y && y != nMinus1) {
273 int j = 1;
274 while (j <= s - 1 && nMinus1 != y) {
275 final long square = ((long) y) * y;
276 y = (int) (square % n);
277 if (1 == y) {
278 return false;
279 } // definitely composite
280 j++;
281 }
282 if (nMinus1 != y) {
283 return false;
284 } // definitely composite
285 }
286 }
287 return true; // definitely prime
288 }
289 }