/symja_android_library/matheclipse-gpl/src/main/java/de/tilman_neumann/jml/factor/cfrac/CFrac63.java
Java | 378 lines | 223 code | 34 blank | 121 comment | 36 complexity | b98caa56fbfc8f770dce8ddc9f1511a4 MD5 | raw file
- /*
- * java-math-library is a Java library focused on number theory, but not necessarily limited to it. It is based on the PSIQS 4.0 factoring project.
- * Copyright (C) 2018 Tilman Neumann (www.tilman-neumann.de)
- *
- * This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
- * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License along with this program;
- * if not, see <http://www.gnu.org/licenses/>.
- */
- package de.tilman_neumann.jml.factor.cfrac;
- import java.math.BigInteger;
- import java.util.HashSet;
- import java.util.Iterator;
- import java.util.TreeMap;
- import org.apache.log4j.Logger;
- import de.tilman_neumann.jml.factor.FactorAlgorithm;
- import de.tilman_neumann.jml.factor.FactorException;
- import de.tilman_neumann.jml.factor.base.PrimeBaseGenerator;
- import de.tilman_neumann.jml.factor.base.congruence.AQPair;
- import de.tilman_neumann.jml.factor.base.congruence.CongruenceCollector;
- import de.tilman_neumann.jml.factor.base.matrixSolver.FactorTest;
- import de.tilman_neumann.jml.factor.base.matrixSolver.FactorTest01;
- import de.tilman_neumann.jml.factor.base.matrixSolver.MatrixSolver;
- import de.tilman_neumann.jml.factor.cfrac.tdiv.TDiv_CF63;
- import static de.tilman_neumann.jml.base.BigIntConstants.*;
- import static de.tilman_neumann.jml.factor.base.GlobalFactoringOptions.ANALYZE;
- /**
- * 63 bit CFrac with Knuth-Schroeppel multiplier.
- *
- * @author Tilman Neumann
- */
- public class CFrac63 extends FactorAlgorithm {
- private static final Logger LOG = Logger.getLogger(CFrac63.class);
- private static final boolean DEBUG = false;
- // input
- private BigInteger N, kN;
- private long floor_sqrt_kN;
- /**
- * Test all Q or only those Q_i+1 with odd i ? "All Q" boosts performance for N > 45 bits
- * approximately
- */
- private boolean use_all_i; //
- // maximum number of SquFoF iterations for a single k
- private int stopRoot;
- private float stopMult;
- private long maxI;
- // multiplier
- private int ks_adjust;
- private KnuthSchroeppel_CFrac ks = new KnuthSchroeppel_CFrac();
- // prime base
- private float C;
- private int primeBaseSize;
- private PrimeBaseGenerator primeBaseBuilder = new PrimeBaseGenerator();
- private int[] primesArray;
- /**
- * The union of all reduced prime bases -- this is the number of variables in the equation system
- */
- private HashSet<Integer> combinedPrimesSet;
- // factorizer for Q
- private TDiv_CF63 auxFactorizer;
- private float maxQRestExponent;
- // collects the congruences we find
- private CongruenceCollector congruenceCollector;
- // the number of congruences we need to find before we try to solve the smooth congruence equation
- // system
- private int requiredSmoothCongruenceCount;
- // extra congruences to have a bigger chance that the equation system solves. the likelihood is >=
- // 1-2^(extraCongruences+1)
- private int extraCongruences;
- /** The solver used for smooth congruence equation systems. */
- private MatrixSolver matrixSolver;
- // time capturing
- private long startTime, linAlgStartTime;
- /**
- * Standard constructor.
- *
- * @param use_all_i
- * @param stopRoot order of the root to compute the maximum number of iterations
- * @param stopMult multiplier to compute the maximum number of iterations
- * @param C multiplier for prime base size
- * @param maxQRestExponent
- * @param auxFactorizer the algorithm to find smooth Q
- * @param extraCongruences the number of surplus congruences we collect to have a greater chance
- * that the equation system solves.
- * @param matrixSolver matrix solver for the smooth congruence equation system
- * @param ks_adjust
- */
- public CFrac63(
- boolean use_all_i,
- int stopRoot,
- float stopMult,
- float C,
- float maxQRestExponent,
- TDiv_CF63 auxFactorizer,
- int extraCongruences,
- MatrixSolver matrixSolver,
- int ks_adjust) {
- this.use_all_i = use_all_i;
- this.stopRoot = stopRoot;
- this.stopMult = stopMult;
- this.C = C;
- this.maxQRestExponent = maxQRestExponent;
- this.auxFactorizer = auxFactorizer;
- this.congruenceCollector = new CongruenceCollector();
- this.extraCongruences = extraCongruences;
- this.matrixSolver = matrixSolver;
- this.ks_adjust = ks_adjust;
- }
- @Override
- public String getName() {
- return "CFrac63(all_i="
- + use_all_i
- + ", ks_adjust="
- + ks_adjust
- + ", stop=("
- + stopRoot
- + ", "
- + stopMult
- + "), C="
- + C
- + ", maxSuSmoothExp="
- + maxQRestExponent
- + ", "
- + auxFactorizer.getName()
- + ")";
- }
- /**
- * Test the current N.
- *
- * @return factor, or null if no factor was found.
- */
- public BigInteger findSingleFactor(BigInteger N) {
- if (ANALYZE) startTime = System.currentTimeMillis();
- this.N = N;
- // compute prime base size
- double N_dbl = N.doubleValue();
- double lnN = Math.log(N_dbl);
- double lnlnN = Math.log(lnN);
- double lnNPow = 0.666667; // heuristics for CFrac
- // we want that the exponents of lnN and lnlnN sum to 1
- this.primeBaseSize =
- 25 + (int) (Math.exp(Math.pow(lnN, lnNPow) * Math.pow(lnlnN, 1 - lnNPow) * C));
- if (DEBUG) LOG.debug("N = " + N + ": primeBaseSize = " + primeBaseSize);
- this.primesArray = new int[primeBaseSize];
- // compute the biggest unfactored rest where some Q is still considered "sufficiently smooth".
- double maxQRest = Math.pow(N_dbl, maxQRestExponent);
- // initialize sub-algorithms for N
- this.auxFactorizer.initialize(N, maxQRest);
- FactorTest factorTest = new FactorTest01(N);
- this.congruenceCollector.initialize(N, factorTest);
- this.combinedPrimesSet = new HashSet<Integer>();
- this.matrixSolver.initialize(N, factorTest);
- // max iterations: there is no need to account for k, because expansions of smooth kN are
- // typically not longer than those for N.
- // long maxI is sufficient to hold any practicable number of iteration steps (~9.22*10^18);
- // stopRoot must be chosen appropriately such that there is no overflow of long values.
- // with stopRoot=5, the overflow of longs starts at N~6.67 * 10^94...
- this.maxI = (long) (stopMult * Math.pow(N_dbl, 1.0 / stopRoot));
- // compute multiplier k: though k=1 is better than Knuth-Schroeppel for N<47 bits,
- // we can ignore that here because that is far out of the optimal CFrac range
- TreeMap<Double, Integer> kMap = ks.computeMultiplier(N, ks_adjust);
- Iterator<Integer> kIter = kMap.values().iterator();
- while (kIter.hasNext()) {
- // get a new k, return immediately if kN is square
- this.kN = BigInteger.valueOf(kIter.next()).multiply(N);
- floor_sqrt_kN =
- (long)
- Math.sqrt(
- kN
- .doubleValue()); // faster than BigInteger sqrt; but the type conversion may
- // give ceil(sqrt(kN)) for some N >= 54 bit
- BigInteger sqrt_kN_big = BigInteger.valueOf(floor_sqrt_kN);
- long diff = kN.subtract(sqrt_kN_big.multiply(sqrt_kN_big)).longValue();
- if (diff == 0) return N.gcd(sqrt_kN_big);
- if (diff < 0) {
- // floor_sqrt_kN was too big, diff too small -> compute correction
- floor_sqrt_kN--;
- diff += (floor_sqrt_kN << 1) + 1;
- }
- // Create the reduced prime base for kN:
- primeBaseBuilder.computeReducedPrimeBase(kN, primeBaseSize, primesArray);
- // add new reduced prime base to the combined prime base
- for (int i = 0; i < primeBaseSize; i++) combinedPrimesSet.add(primesArray[i]);
- // we want: #equations = #variables + some extra congruences
- this.requiredSmoothCongruenceCount = combinedPrimesSet.size() + extraCongruences;
- try {
- // initialize the Q-factorizer with new prime base
- this.auxFactorizer.initialize(kN, primeBaseSize, primesArray); // may throw FactorException
- // search square Q_i
- test(diff);
- } catch (FactorException fe) {
- if (ANALYZE) {
- long endTime = System.currentTimeMillis();
- LOG.info(
- "Found factor of N="
- + N
- + " in "
- + (endTime - startTime)
- + "ms (LinAlgPhase took "
- + (endTime - linAlgStartTime)
- + "ms)");
- }
- return fe.getFactor();
- }
- }
- return I_1; // fail, too few Knuth-Schroeppel multipliers
- }
- protected void test(long Q_ip1) throws FactorException {
- // initialization for first iteration step
- long i = 0;
- BigInteger A_im2 = null;
- BigInteger A_im1 = I_1;
- BigInteger A_i = BigInteger.valueOf(floor_sqrt_kN);
- long P_im1 = 1;
- long P_i = floor_sqrt_kN;
- long Q_i = 1;
- // first iteration step
- long two_floor_sqrt_kN = floor_sqrt_kN << 1;
- while (true) {
- if (DEBUG) verifyCongruence(i, A_i, BigInteger.valueOf(Q_ip1));
- // [McMath 2004] points out (on SquFoF) that we have to look for square Q_i at some even i.
- // Here I test Q_i+1, so I have to look for square Q_i+1 at odd i.
- // In CFRAC, square congruences are also tested in the CongruenceCollector,
- // but doing it here before trial division is good for performance.
- boolean isSquare = false;
- if (i % 2 == 1) {
- long Q_ip1_sqrt = (long) Math.sqrt(Q_ip1);
- if (Q_ip1_sqrt * Q_ip1_sqrt == Q_ip1) {
- // Q_i+1 is square -> test gcd
- BigInteger gcd = N.gcd(A_i.subtract(BigInteger.valueOf(Q_ip1_sqrt)));
- if (gcd.compareTo(I_1) > 0 && gcd.compareTo(N) < 0) throw new FactorException(gcd);
- isSquare = true;
- }
- }
- if (isSquare == false && (use_all_i || i % 2 == 1)) {
- // Q_i+1 is not square and the i is right, too -> check for smooth relations.
- // Here a constraint on the size of Q would mean a severe performance penalty!
- long Q_test = i % 2 == 1 ? Q_ip1 : -Q_ip1; // make Q congruent A^2
- AQPair aqPair = auxFactorizer.test(A_i, Q_test);
- if (DEBUG) LOG.debug("N = " + N + ": Q_test = " + Q_test + " -> aqPair = " + aqPair);
- if (aqPair != null) {
- // the Q was sufficiently smooth
- if (ANALYZE)
- linAlgStartTime = System.currentTimeMillis(); // just in case it really starts
- boolean addedSmooth = congruenceCollector.add(aqPair);
- if (addedSmooth) {
- int smoothCongruenceCount = congruenceCollector.getSmoothCongruenceCount();
- if (smoothCongruenceCount >= requiredSmoothCongruenceCount) {
- // try to solve equation system
- // LOG.debug("#smooths = " + smoothCongruenceCount + ", #requiredSmooths = " +
- // requiredSmoothCongruenceCount);
- matrixSolver.solve(
- congruenceCollector.getSmoothCongruences()); // throws FactorException
- // no factor exception -> extend equation system and continue searching smooth
- // congruences
- requiredSmoothCongruenceCount += extraCongruences;
- }
- }
- }
- }
- // exit loop ?
- if (++i == maxI) return;
- // keep values from last round
- A_im2 = A_im1;
- A_im1 = A_i;
- P_im1 = P_i;
- long Q_im1 = Q_i;
- Q_i = Q_ip1;
- // compute next values
- long b_i = (floor_sqrt_kN + P_im1) / Q_i; // floor(rational result)
- P_i = b_i * Q_i - P_im1;
- Q_ip1 = Q_im1 + b_i * (P_im1 - P_i);
- // carry along A_i % N from continuant recurrence
- A_i = addModN(mulModN(b_i, A_im1), A_im2);
- // stop when continuant period is complete
- if (b_i == two_floor_sqrt_kN) return;
- }
- }
- /**
- * Verify congruence A_i^2 == (-1)^(i+1)*Q_i+1 (mod N)
- *
- * @param i
- * @param A_i
- * @param Q_ip1
- */
- private void verifyCongruence(long i, BigInteger A_i, BigInteger Q_ip1) {
- // assertTrue(Q_ip1.signum()>=0);
- // verify congruence A^2 == Q (mod N)
- BigInteger Q_test = i % 2 == 1 ? Q_ip1 : Q_ip1.negate().mod(N);
- BigInteger div[] = A_i.pow(2).subtract(Q_test).divideAndRemainder(N);
- // assertEquals(I_0, div[1]); // works
- LOG.debug("A^2-Q = " + div[0] + " * N");
- LOG.debug("A^2 % N = " + A_i.pow(2).mod(N) + ", Q = " + Q_test);
- // assertEquals(Q_test, A_i.pow(2).mod(N)); // works
- }
- /**
- * Addition modulo N, with <code>a, b < N</code>.
- *
- * @param a
- * @param b
- * @return (a+b) mod N
- */
- private BigInteger addModN(BigInteger a, BigInteger b) {
- BigInteger sum = a.add(b);
- return sum.compareTo(N) < 0 ? sum : sum.subtract(N);
- }
- /**
- * Multiplication (m*a) modulo N, with m often small and <code>a < N</code>.
- *
- * @param m
- * @param a
- * @return (m*a) mod N
- */
- private BigInteger mulModN(long m, BigInteger a) {
- if (m < 4) { // 0, 1, 10, 11
- int m_int = (int) m;
- switch (m_int) {
- case 0:
- return I_0;
- case 1:
- return a;
- case 2:
- {
- BigInteger two_a = a.shiftLeft(1); // faster than 2*a or a+a
- return two_a.compareTo(N) < 0 ? two_a : two_a.subtract(N);
- }
- case 3:
- {
- BigInteger ma_modN = a.shiftLeft(1).add(a); // < 3*N
- if (ma_modN.compareTo(N) < 0) return ma_modN;
- ma_modN = ma_modN.subtract(N); // < 2*N
- return ma_modN.compareTo(N) < 0 ? ma_modN : ma_modN.subtract(N);
- }
- }
- // adding case 4 does not help because then bitLength() does not fit exactly
- }
- BigInteger product = a.multiply(BigInteger.valueOf(m));
- return product.compareTo(N) < 0 ? product : product.mod(N);
- }
- }