[NUMBERS-133] Use iteration algorithm from bounded trial division in Primes.nextPrime(int)

commons-numbers-primes/src/main/java/org/apache/commons/numbers/primes/Primes.java

@@ -17,7 +17,10 @@
 package org.apache.commons.numbers.primes;
 
 import java.text.MessageFormat;
+
+import java.util.Arrays;
 import java.util.List;
+import java.util.PrimitiveIterator;
 
 
 /**
@@ -70,37 +73,21 @@ public static boolean isPrime(int n) {
  public static int nextPrime(int n) {
  if (n < 0) {
  throw new IllegalArgumentException(MessageFormat.format(NUMBER_TOO_SMALL, n, 0));
- }
- if (n == 2) {
- return 2;
- }
- n |= 1; // make sure n is odd
- if (n == 1) {
- return 2;
- }
-
- if (isPrime(n)) {
- return n;
- }
-
- // prepare entry in the +2, +4 loop:
- // n should not be a multiple of 3
- final int rem = n % 3;
- if (0 == rem) { // if n % 3 == 0
- n += 2; // n % 3 == 2
- } else if (1 == rem) { // if n % 3 == 1
- // if (isPrime(n)) return n;
- n += 4; // n % 3 == 2
- }
- while (true) { // this loop skips all multiple of 3
- if (isPrime(n)) {
- return n;
+ } else if (n <= SmallPrimes.PRIMES_LAST) {
+ int index = Arrays.binarySearch(SmallPrimes.PRIMES, n);
+ if (index < 0) {
+ index = - (index + 1);
  }
- n += 2; // n % 3 == 1
- if (isPrime(n)) {
- return n;
+ return SmallPrimes.PRIMES[index];
+ } else {
+ PrimitiveIterator.OfInt potentialPrimesIterator = SmallPrimes.potentialPrimesGTE(n);
+ while (true) {
+ // Integer.MAX_VALUE is a prime number, so no risk of overflow
+ int candidate = potentialPrimesIterator.next();
+ if (isPrime(candidate)) {
+ return candidate;
+ }
  }
- n += 4; // n % 3 == 2
  }
  }
 

commons-numbers-primes/src/main/java/org/apache/commons/numbers/primes/SmallPrimes.java

@@ -24,14 +24,15 @@
 import java.util.HashSet;
 import java.util.List;
 import java.util.Map.Entry;
+import java.util.PrimitiveIterator;
 import java.util.Set;
 
 /**
  * Utility methods to work on primes within the <code>int</code> range.
  */
 class SmallPrimes {
  /**
- * The first 512 prime numbers.
+ * The first 512 prime numbers, in ascending order.
  * <p>
  * It contains all primes smaller or equal to the cubic square of Integer.MAX_VALUE.
  * As a result, <code>int</code> numbers which are not reduced by those primes are guaranteed
@@ -73,38 +74,42 @@ class SmallPrimes {
  * 0 (inclusive) and the least common multiple, i.e. the product, of those
  * prime numbers (exclusive) that are not divisible by any of these prime
  * numbers. The prime numbers in the set are among the first 512 prime
- * numbers, and the {@code int} array containing the numbers undivisible by
+ * numbers, and the {@code int} array containing the numbers indivisible by
  * these prime numbers is sorted in ascending order.
  *
- * <p>The purpose of this field is to speed up trial division by skipping
- * multiples of individual prime numbers, specifically those contained
- * in the key of this {@code Entry}, by only trying integers that are equivalent
- * to one of the integers contained in the value of this {@code Entry} modulo
- * the least common multiple of the prime numbers in the set.</p>
+ * <p>The purpose of this field is to facilitate iterating over potential
+ * prime numbers by skipping multiples of individual primes, specifically
+ * those contained in the key of this {@code Entry}, by only selecting
+ * integers that are equivalent to one of the integers contained in the
+ * value of this {@code Entry} modulo the least common multiple of the prime
+ * numbers in the set.</p>
  *
- * <p>Note that, if {@code product} is the product of the prime numbers,
+ * <p>Note that, if \(\mathit{product}\) is the product of the prime numbers,
  * the last number in the array of coprime integers is necessarily
- * {@code product - 1}, because if {@code product ≡ 0 mod p}, then
- * {@code product - 1 ≡ -1 mod p}, and {@code 0 ≢ -1 mod p} for any prime number p.</p>
+ * \(\mathit{product} - 1\), because if \(\mathit{product} \equiv 0 \pmod{p}\),
+ * then \(\mathit{product} - 1 \equiv -1 \pmod{p}\), and
+ * \(0 \not\equiv -1 \pmod{p}\) for any prime number \(p\).</p>
+ *
+ * @see #potentialPrimesGTE(int)
  */
  static final Entry<Set<Integer>, int[]> PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES;
 
  static {
  /*
- According to the Chinese Remainder Theorem, for every combination of
- congruence classes modulo distinct, pairwise coprime moduli, there
- exists exactly one congruence class modulo the product of these
- moduli that is contained in every one of the former congruence
- classes. Since the number of congruence classes coprime to a prime
- number p is p-1, the number of congruence classes coprime to all
- prime numbers p_1, p_2, p_3 … is (p_1 - 1) * (p_2 - 1) * (p_3 - 1) …
-
- Therefore, when using the first five prime numbers as those whose multiples
- are to be skipped in trial division, the array containing the coprime
- equivalence classes will have to hold (2-1)*(3-1)*(5-1)*(7-1)*(11-1) = 480
- values. As a consequence, the amount of integers to be tried in
- trial division is reduced to 480/(2*3*5*7*11), which is about 20.78%,
- of all integers.
+  * According to the Chinese Remainder Theorem, for every combination of
+  * congruence classes modulo distinct, pairwise coprime moduli, there
+  * exists exactly one congruence class modulo the product of these
+  * moduli that is contained in every one of the former congruence
+  * classes. Since the number of congruence classes coprime to a prime
+  * number p is p-1, the number of congruence classes coprime to all
+  * prime numbers p_1, p_2, p_3 … is (p_1 - 1) * (p_2 - 1) * (p_3 - 1) …
+ *
+  * Therefore, when choosing the first five prime numbers as those whose
+  * multiples should be skipped during an iteration, the array containing
+  * the coprime equivalence classes will have to hold
+  * (2-1)*(3-1)*(5-1)*(7-1)*(11-1) = 480 values. As a consequence, the
+  * amount of potential prime numbers is reduced to 480/(2*3*5*7*11),
+  * which is about 20.78%, of all integers.
  */
  Set<Integer> primeNumbers = new HashSet<>();
  primeNumbers.add(Integer.valueOf(2));
@@ -171,43 +176,18 @@ static int smallTrialDivision(int n,
  static int boundedTrialDivision(int n,
  int maxFactor,
  List<Integer> factors) {
- int minFactor = PRIMES_LAST + 2;
-
- /*
- only trying integers of the form k*m + c, where k >= 0, m is the
- product of some prime numbers which n is required not to contain
- as prime factors, and c is an integer undivisible by all of those
- prime numbers; in other words, skipping multiples of these primes
- */
- int m = PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue()[PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue().length - 1] + 1;
- int km = m * (minFactor / m);
- int currentEquivalenceClassIndex = Arrays.binarySearch(
- PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue(),
- minFactor % m);
-
- /*
- Since minFactor is the next smallest prime number after the
- first 512 primes, it cannot be a multiple of one of them, therefore,
- the index returned by the above binary search must be non-negative.
- */
+ PrimitiveIterator.OfInt potentialPrimesIterator = potentialPrimesGTE(PRIMES_LAST + 2);
 
  boolean done = false;
  while (!done) {
  // no check is done about n >= f
- int f = km + PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue()[currentEquivalenceClassIndex];
+ int f = potentialPrimesIterator.next();
  if (f > maxFactor) {
  done = true;
  } else if (0 == n % f) {
  n /= f;
  factors.add(f);
  done = true;
- } else {
- if (currentEquivalenceClassIndex == PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue().length - 1) {
- km += m;
- currentEquivalenceClassIndex = 0;
- } else {
- currentEquivalenceClassIndex++;
- }
  }
  }
  if (n != 1) {
@@ -284,4 +264,71 @@ static boolean millerRabinPrimeTest(final int n) {
  }
  return true; // definitely prime
  }
+
+ /**
+ * Returns an iterator that iterates, in ascending oder, over all positive
+ * integers greater than or equal to ("GTE") the passed lower bound, skipping
+ * numbers that are a multiple of any of the prime numbers contained in the
+ * key of {@link #PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES}.
+ *
+ * <p>
+ * NOTE: The iterator returned by this method will be completely oblivious
+ * to the inherent limitation of the {@code int} range and will simply
+ * continue iterating past {@code Integer.MAX_VALUE}, producing an infinite
+ * sequence of invalid results from then on. It is therefore the
+ * responsibility of the caller to account for the possibility of overflow
+ * (the first invalid result will be negative).
+ * </p>
+ *
+ * @param lowerBound the lower bound for the iteration results
+ * @return an iterator as described above
+ */
+ static PrimitiveIterator.OfInt potentialPrimesGTE(final int lowerBound) {
+ /*
+ * The numbers that should be iterated over are of the form k*m + c,
+ * where k >= 0, m is the least common multiple, that is, the product of
+ * the prime numbers whose multiples should be skipped, and c is a
+ * positive integer between 0 and m that is indivisible by any of these
+ * prime numbers.
+ */
+ return new PrimitiveIterator.OfInt() {
+ private final int m;
+ private int kTimesM;
+ private int indexForNextC;
+
+ {
+ m = PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue()[PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue().length - 1] + 1;
+
+ int initialValue = Math.max(lowerBound, 1);
+ int remainder = initialValue % m;
+ kTimesM = initialValue - remainder;
+
+ indexForNextC = Arrays.binarySearch(
+ PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue(),
+ remainder
+ );
+ if (indexForNextC < 0) {
+ // the last element in the array is m-1, so the remainder
+ // cannot be greater than the last element in the array
+ indexForNextC = - (indexForNextC + 1);
+ }
+ }
+
+ @Override
+ public int nextInt() {
+ int next = kTimesM + PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue()[indexForNextC];
+ indexForNextC++;
+ if (indexForNextC == PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getValue().length) {
+ indexForNextC = 0;
+ kTimesM += m;
+ }
+ return next;
+ }
+
+ @Override
+ public boolean hasNext() {
+ return true;
+ }
+ };
+ }
 }

commons-numbers-primes/src/test/java/org/apache/commons/numbers/primes/SmallPrimesTest.java

@@ -20,6 +20,7 @@
 import java.util.Arrays;
 import java.util.Collections;
 import java.util.List;
+import java.util.PrimitiveIterator;
 
 import org.junit.jupiter.api.Assertions;
 import org.junit.jupiter.api.Test;
@@ -141,4 +142,46 @@ public void millerRabinPrimeTest_composites() {
  }
  }
  }
+
+ @Test
+ public void testPotentialPrimes() {
+ testPotentialPrimes(-23456, 1500);
+ testPotentialPrimes(34567, 1500);
+ testPotentialPrimes(Integer.MAX_VALUE - 97, 1500);
+ }
+
+ private void testPotentialPrimes(int lowerBound, int maxIterations) {
+ PrimitiveIterator.OfInt potentialPrimesIterator = SmallPrimes.potentialPrimesGTE(lowerBound);
+
+ int iterationCount = 0;
+ boolean overflow = false;
+ int previous = Math.max(1, lowerBound) - 1;
+ while (iterationCount != maxIterations && !overflow) {
+ Assertions.assertTrue(potentialPrimesIterator.hasNext());
+ int next = potentialPrimesIterator.nextInt();
+ if (next < 0) {
+ overflow = true;
+ }
+
+ for (int i = previous + 1;
+ i >= 0 && i <= (overflow ? Integer.MAX_VALUE : next - 1);
+ i++) {
+ boolean hasPrimeFactor = false;
+ for (Integer prime : SmallPrimes.PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getKey()) {
+ if (i % prime == 0) {
+ hasPrimeFactor = true;
+ break;
+ }
+ }
+ Assertions.assertTrue(hasPrimeFactor);
+ }
+ if (!overflow) {
+ for (Integer prime : SmallPrimes.PRIME_NUMBERS_AND_COPRIME_EQUIVALENCE_CLASSES.getKey()) {
+ Assertions.assertNotEquals(0, next % prime);
+ }
+ }
+ previous = next;
+ iterationCount++;
+ }
+ }
 }