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(2);
111 primeNumbers.add(3);
112 primeNumbers.add(5);
113 primeNumbers.add(7);
114 primeNumbers.add(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 */
170 static void boundedTrialDivision(int n,
171 int maxFactor,
172 List<Integer> factors) {
173 final int minFactor = PRIMES_LAST + 2;
174
175 /*
176 only trying integers of the form k*m + c, where k >= 0, m is the
177 product of some prime numbers which n is required not to contain
178 as prime factors, and c is an integer undivisible by all of those
179 prime numbers; in other words, skipping multiples of these primes
180 */
181 final int[] a = PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue();
182 final int m = a[a.length - 1] + 1;
183 int km = m * (minFactor / m);
184 int currentEquivalenceClassIndex = Arrays.binarySearch(a, minFactor % m);
185
186 /*
187 Since minFactor is the next smallest prime number after the
188 first 512 primes, it cannot be a multiple of one of them, therefore,
189 the index returned by the above binary search must be non-negative.
190 */
191
192 boolean done = false;
193 while (!done) {
194 // no check is done about n >= f
195 final int f = km + a[currentEquivalenceClassIndex];
196 if (f > maxFactor) {
197 done = true;
198 } else if (0 == n % f) {
199 n /= f;
200 factors.add(f);
201 done = true;
202 } else {
203 if (currentEquivalenceClassIndex == a.length - 1) {
204 km += m;
205 currentEquivalenceClassIndex = 0;
206 } else {
207 currentEquivalenceClassIndex++;
208 }
209 }
210 }
211 // Note: When this method is used after the small trial division n != 1
212 // for all n in [2, MAX_VALUE] (tested using an exhaustive enumeration).
213 factors.add(n);
214 }
215
216 /**
217 * Factorization by trial division.
218 *
219 * @param n Number to factor.
220 * @return the list of prime factors of {@code n}.
221 */
222 static List<Integer> trialDivision(int n) {
223 final List<Integer> factors = new ArrayList<>(32);
224 n = smallTrialDivision(n, factors);
225 if (1 == n) {
226 return factors;
227 }
228 // here we are sure that n is either a prime or a semi prime
229 final int bound = (int) Math.sqrt(n);
230 boundedTrialDivision(n, bound, factors);
231 return factors;
232 }
233
234 /**
235 * Miller-Rabin probabilistic primality test for int type, used in such
236 * a way that a result is always guaranteed.
237 * <p>
238 * It uses the prime numbers as successive base therefore it is guaranteed
239 * to be always correct (see Handbook of applied cryptography by Menezes,
240 * table 4.1).
241 *
242 * @param n Number to test: an odd integer ≥ 3.
243 * @return true if {@code n} is prime, false if it is definitely composite.
244 */
245 static boolean millerRabinPrimeTest(final int n) {
246 final int nMinus1 = n - 1;
247 final int s = Integer.numberOfTrailingZeros(nMinus1);
248 final int r = nMinus1 >> s;
249 // r must be odd, it is not checked here
250 int t = 1;
251 if (n >= 2047) {
252 t = 2;
253 }
254 if (n >= 1373653) {
255 t = 3;
256 }
257 if (n >= 25326001) {
258 t = 4;
259 } // works up to 3.2 billion, int range stops at 2.7 so we are safe :-)
260 final BigInteger br = BigInteger.valueOf(r);
261 final BigInteger bn = BigInteger.valueOf(n);
262
263 for (int i = 0; i < t; i++) {
264 final BigInteger a = BigInteger.valueOf(PRIMES[i]);
265 final BigInteger bPow = a.modPow(br, bn);
266 int y = bPow.intValue();
267 if (1 != y && y != nMinus1) {
268 int j = 1;
269 while (j <= s - 1 && nMinus1 != y) {
270 final long square = ((long) y) * y;
271 y = (int) (square % n);
272 if (1 == y) {
273 return false;
274 } // definitely composite
275 j++;
276 }
277 if (nMinus1 != y) {
278 return false;
279 } // definitely composite
280 }
281 }
282 return true; // definitely prime
283 }
284 }