diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/primes/PollardRho.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/primes/PollardRho.java | 165 |
1 files changed, 165 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/primes/PollardRho.java b/src/main/java/org/apache/commons/math3/primes/PollardRho.java new file mode 100644 index 0000000..4cbc064 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/primes/PollardRho.java @@ -0,0 +1,165 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.primes; + +import org.apache.commons.math3.util.FastMath; + +import java.util.ArrayList; +import java.util.List; + +/** + * Implementation of the Pollard's rho factorization algorithm. + * + * @since 3.2 + */ +class PollardRho { + + /** Hide utility class. */ + private PollardRho() {} + + /** + * Factorization using Pollard's rho algorithm. + * + * @param n number to factors, must be > 0 + * @return the list of prime factors of n. + */ + public static List<Integer> primeFactors(int n) { + final List<Integer> factors = new ArrayList<Integer>(); + + n = SmallPrimes.smallTrialDivision(n, factors); + if (1 == n) { + return factors; + } + + if (SmallPrimes.millerRabinPrimeTest(n)) { + factors.add(n); + return factors; + } + + int divisor = rhoBrent(n); + factors.add(divisor); + factors.add(n / divisor); + return factors; + } + + /** + * Implementation of the Pollard's rho factorization algorithm. + * + * <p>This implementation follows the paper "An improved Monte Carlo factorization algorithm" by + * Richard P. Brent. This avoids the triple computation of f(x) typically found in Pollard's rho + * implementations. It also batches several gcd computation into 1. + * + * <p>The backtracking is not implemented as we deal only with semi-primes. + * + * @param n number to factor, must be semi-prime. + * @return a prime factor of n. + */ + static int rhoBrent(final int n) { + final int x0 = 2; + final int m = 25; + int cst = SmallPrimes.PRIMES_LAST; + int y = x0; + int r = 1; + do { + int x = y; + for (int i = 0; i < r; i++) { + final long y2 = ((long) y) * y; + y = (int) ((y2 + cst) % n); + } + int k = 0; + do { + final int bound = FastMath.min(m, r - k); + int q = 1; + for (int i = -3; + i < bound; + i++) { // start at -3 to ensure we enter this loop at least 3 times + final long y2 = ((long) y) * y; + y = (int) ((y2 + cst) % n); + final long divisor = FastMath.abs(x - y); + if (0 == divisor) { + cst += SmallPrimes.PRIMES_LAST; + k = -m; + y = x0; + r = 1; + break; + } + final long prod = divisor * q; + q = (int) (prod % n); + if (0 == q) { + return gcdPositive(FastMath.abs((int) divisor), n); + } + } + final int out = gcdPositive(FastMath.abs(q), n); + if (1 != out) { + return out; + } + k += m; + } while (k < r); + r = 2 * r; + } while (true); + } + + /** + * Gcd between two positive numbers. + * + * <p>Gets the greatest common divisor of two numbers, using the "binary gcd" method, which + * avoids division and modulo operations. See Knuth 4.5.2 algorithm B. This algorithm is due to + * Josef Stein (1961). Special cases: + * + * <ul> + * <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and {@code gcd(x, 0)} is the value + * of {@code x}. + * <li>The invocation {@code gcd(0, 0)} is the only one which returns {@code 0}. + * </ul> + * + * @param a first number, must be ≥ 0 + * @param b second number, must be ≥ 0 + * @return gcd(a,b) + */ + static int gcdPositive(int a, int b) { + // both a and b must be positive, it is not checked here + // gdc(a,0) = a + if (a == 0) { + return b; + } else if (b == 0) { + return a; + } + + // make a and b odd, keep in mind the common power of twos + final int aTwos = Integer.numberOfTrailingZeros(a); + a >>= aTwos; + final int bTwos = Integer.numberOfTrailingZeros(b); + b >>= bTwos; + final int shift = FastMath.min(aTwos, bTwos); + + // a and b >0 + // if a > b then gdc(a,b) = gcd(a-b,b) + // if a < b then gcd(a,b) = gcd(b-a,a) + // so next a is the absolute difference and next b is the minimum of current values + while (a != b) { + final int delta = a - b; + b = FastMath.min(a, b); + a = FastMath.abs(delta); + // for speed optimization: + // remove any power of two in a as b is guaranteed to be odd throughout all iterations + a >>= Integer.numberOfTrailingZeros(a); + } + + // gcd(a,a) = a, just "add" the common power of twos + return a << shift; + } +} |