v1.4.0 tilman_neumann.jml java-math-library from github repo

- see:
https://github.com/TilmanNeumann/java-math-library/releases/tag/v1.4.0

symja_android_library/matheclipse-external/pom.xml

@@ -44,6 +44,12 @@
 			<groupId>org.apache.logging.log4j</groupId>
 			<artifactId>log4j-api</artifactId>
 		</dependency>
+
+		<dependency>
+			<groupId>org.hipparchus</groupId>
+			<artifactId>hipparchus-core</artifactId>
+		    <version>4.0-SNAPSHOT</version>
+		</dependency>
 		<!-- test dependencies -->
 		<dependency>
 			<groupId>${project.groupId}</groupId>

symja_android_library/matheclipse-gpl/src/main/java/de/tilman_neumann/jml/base/Uint128.java

@@ -350,28 +350,6 @@ public static Uint128 square64(long a) {
 		return new Uint128(r_hi, r_lo);
 	}
 
-	/**
-	 * Computes the low part of the product of two unsigned 64 bit integers.
-	 * 
-	 * Overflows of the "middle term" are not interesting here because they'ld only
-	 * affect the high part of the multiplication result.
-	 * 
-	 * @param a
-	 * @param b
-	 * @return (a*b) & 0xFFFFFFFFL
-	 */
-	// XXX a*b should give the same result !?
-	public static long mul64_getLow(long a, long b) {
-		final long a_hi = a >>> 32;
-		final long b_hi = b >>> 32;
-		final long a_lo = a & 0xFFFFFFFFL;
-		final long b_lo = b & 0xFFFFFFFFL;
-		final long lo_prod = a_lo * b_lo;
-		final long med_term = a_hi * b_lo + a_lo * b_hi;
-		final long r_lo = ((med_term & 0xFFFFFFFFL) << 32) + lo_prod;
-		return r_lo;
-	}
-
 	/**
 	 * Multiplication of two unsigned 128 bit integers.
 	 * 

symja_android_library/matheclipse-gpl/src/main/java/de/tilman_neumann/jml/factor/CombinedFactorAlgorithm.java

@@ -1,6 +1,6 @@
 /*
  * 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-2024 Tilman Neumann - [email protected]
+ * Copyright (C) 2018-2025 Tilman Neumann - [email protected]
  *
  * 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.
@@ -61,7 +61,7 @@ public class CombinedFactorAlgorithm extends FactorAlgorithm {
 
 	private TDiv31Barrett tDiv31 = new TDiv31Barrett();
 	private HartFast2Mult hart = new HartFast2Mult(true); // for general factor arguments, trial division is needed
-	private TinyEcm64MHInlined tinyEcm = new TinyEcm64MHInlined();
+	private TinyEcm64MHInlined tinyEcm = new TinyEcm64MHInlined(true); // for general factor arguments, trial division is needed
 	private PollardRhoBrentMontgomery64MH pollardRhoBrentMontgomery64MH = new PollardRhoBrentMontgomery64MH();
 	private TDiv tdiv = new TDiv();
 	private EllipticCurveMethod ecm = new EllipticCurveMethod(0);
@@ -140,6 +140,7 @@ public void searchFactors(FactorArguments args, FactorResult result) {
 		}
 		else if (NBits<46) hart.searchFactors(args, result);
 		else if (NBits<63) tinyEcm.searchFactors(args, result);
+		else if (NBits<64) pollardRhoBrentMontgomery64MH.searchFactors(args, result);
 		else {
 			if (SEARCH_SMALL_FACTORS) {
 				int actualTdivLimit;

symja_android_library/matheclipse-gpl/src/main/java/de/tilman_neumann/jml/factor/FactorAlgorithm.java

@@ -13,13 +13,13 @@
  */
 package de.tilman_neumann.jml.factor;
 
-import static de.tilman_neumann.jml.base.BigIntConstants.I_0;
-import static de.tilman_neumann.jml.base.BigIntConstants.I_1;
-import static de.tilman_neumann.jml.base.BigIntConstants.I_2;
-import static de.tilman_neumann.jml.base.BigIntConstants.I_MINUS_1;
+import static de.tilman_neumann.jml.base.BigIntConstants.*;
+
 import java.math.BigInteger;
-import org.apache.logging.log4j.LogManager;
+
 import org.apache.logging.log4j.Logger;
+import org.apache.logging.log4j.LogManager;
+
 import de.tilman_neumann.jml.factor.base.FactorArguments;
 import de.tilman_neumann.jml.factor.base.FactorResult;
 import de.tilman_neumann.jml.primes.probable.BPSWTest;
@@ -87,7 +87,7 @@ public SortedMultiset<BigInteger> factor(BigInteger N) {
 	 * @param N Number to factor.
 	 * @param primeFactors a map to which found factors are added
 	 */
-    public void factor(BigInteger N, SortedMultiset<BigInteger> primeFactors) {
+	public void factor(BigInteger N, SortedMultiset<BigInteger> primeFactors) {
 		// Make N positive
 		if (N.signum()<0) {
 			primeFactors.add(I_MINUS_1);
@@ -114,7 +114,6 @@ public void factor(BigInteger N, SortedMultiset<BigInteger> primeFactors) {
 		}
 
 		// N contains larger factors...
-		FactorArguments args = new FactorArguments(N, 1);
 		FactorResult factorResult = new FactorResult(primeFactors, new SortedMultiset_BottomUp<BigInteger>(), new SortedMultiset_BottomUp<BigInteger>(), 3);
 		SortedMultiset<BigInteger> untestedFactors = factorResult.untestedFactors; // ArrayList would be faster
 		untestedFactors.add(N);
@@ -149,9 +148,7 @@ public void factor(BigInteger N, SortedMultiset<BigInteger> primeFactors) {
 
 				BigInteger compositeFactor = factorResult.compositeFactors.firstKey();
 				int exp = factorResult.compositeFactors.removeAll(compositeFactor);
-				args.N = compositeFactor;
-				args.NBits = compositeFactor.bitLength();
-				args.exp = exp;
+				FactorArguments args = new FactorArguments(compositeFactor, exp);
 				searchFactors(args, factorResult);
 				if (DEBUG) LOG.debug("3: factorResult: " + factorResult);
 			}

symja_android_library/matheclipse-gpl/src/main/java/de/tilman_neumann/jml/factor/base/matrixSolver/BlockLanczos.java

@@ -586,7 +586,6 @@ private void colexchange(int[] XmY, int[] V, int[] V1, int[] V2, int col1, int c
 	    for (row = V.length - 1; row >= 0; row--) {             
 	    	// If both bits are different toggle them.
 	    	if (((matr1[row] & mask1) == 0) != ((matr2[row] & mask2) == 0)) {
-	    		// If both bits are different toggle them.
 	    		matr1[row] ^= mask1;
 	    		matr2[row] ^= mask2;
 	    	}

symja_android_library/matheclipse-gpl/src/main/java/de/tilman_neumann/jml/factor/ecm/EllipticCurveMethod.java

@@ -980,8 +980,15 @@ void LongToBigNbr(long Nbr, int Out[]) {
 	void BigNbrToBigInt(BigInteger N, int[] TestNbr, int NumberLength) {
 		// Temp needs 1 extra int because of TestNbr[j] = p; in BigNbrToBigInt(long[])
 	    long[] Temp = new long[NumberLength+1]; // zero-initialized
+
+	    // Convert32To31Bits does not work for negative N, so we switch the sign before and after
+	    boolean chSign = N.signum()<0;
+	    if (chSign) N = N.negate();
+
 		BigNbrToBigInt(N, Temp, NumberLength);
 	    Convert32To31Bits(Temp, TestNbr, NumberLength);
+
+		if (chSign) ChSignBigNbr(TestNbr);
 	}
 
 	/**
@@ -1233,7 +1240,6 @@ private void ChSignBigNbr(int Nbr[]) {
 		}
 	}
 
-	// TODO This may be the problem in in-out-conversion31 of negative inputs
 	void Convert31To32Bits(int[] nbr31, long[] nbr32) {
 		int i, j, k;
 		i = 0;
@@ -1257,6 +1263,15 @@ void Convert31To32Bits(int[] nbr31, long[] nbr32) {
 		}
 	}
 
+	/**
+	 * convert nbr32 into nbr31.<br/><br/>
+	 * 
+	 * <strong>Warning: This implementation does not work for negative arguments!</strong>
+	 * 
+	 * @param nbr32
+	 * @param nbr31
+	 * @param NumberLength
+	 */
 	private void Convert32To31Bits(long[] nbr32, int[] nbr31, int NumberLength) {
 		int i, j, k;
 		j = 0;
@@ -1683,7 +1698,14 @@ void ModInvBigNbr(int[] a, int[] b, int[] inv) {
 	    if (i < 0 || B[i] < CalcAuxModInvMu[i]) { // If B < Mu
 	    	SubtractBigNbr32(CalcAuxModInvMu, B, CalcAuxModInvMu); // Mu <- Mu - B
 	    }
+	    // It doesn't matter here that Convert32To31Bits() does not work for negative inputs, because CalcAuxModInvMu is always positive
 	    Convert32To31Bits(CalcAuxModInvMu, inv, NumberLength);
+	    if (DEBUG) {
+		    BigInteger inv32 = BigIntToBigNbr(CalcAuxModInvMu);
+		    Ensure.ensureGreater(inv32.signum(), 0);
+		    BigInteger inv31 = BigIntToBigNbr(inv);
+		    Ensure.ensureEquals(inv32, inv31);
+	    }
 	}
 
 	BigInteger BigIntToBigNbr(int[] nbr) {

symja_android_library/matheclipse-gpl/src/main/java/de/tilman_neumann/jml/factor/ecm/TinyEcm64.java

@@ -43,15 +43,18 @@
 import de.tilman_neumann.jml.factor.FactorAlgorithm;
 import de.tilman_neumann.jml.factor.base.FactorArguments;
 import de.tilman_neumann.jml.factor.base.FactorResult;
-import de.tilman_neumann.jml.factor.tdiv.TDiv;
+import de.tilman_neumann.jml.factor.tdiv.TDiv63Inverse;
 import de.tilman_neumann.jml.gcd.Gcd63;
 import de.tilman_neumann.jml.primes.probable.BPSWTest;
 import de.tilman_neumann.jml.random.SpRand32;
 import de.tilman_neumann.util.Ensure;
 
 /**
  * A port of Ben Buhrow's tinyecm.c (https://www.mersenneforum.org/showpost.php?p=521028&postcount=84),
- * an ECM implementation for unsigned 64 bit integers.
+ * an ECM implementation for unsigned 64 bit integers.<br/><br/>
+ * 
+ * <strong>WARNING: Without trial division, tinyEcm might fail (i.e. run forever) for some N having >=2 small factors, like 735 = 3*5*7^2.</strong>
+ * If you are not sure about the nature of your test numbers, use doTDivFirst==true.
  * 
  * @author Tilman Neumann
  */
@@ -99,7 +102,7 @@ public ecm_work() {
 
 	private static final boolean DEBUG = false;
 
-	private static final int MAX_BITS_SUPPORTED = 62;
+	private static final int MAX_BITS_SUPPORTED = 62; // seems to work for 63 bit numbers now, but very slow - so not completely fixed for that
 
 	// The reducer R is 2^64, but the only constant still required is the half of it.
 	private static final long R_HALF = 1L << 63;
@@ -155,19 +158,30 @@ public ecm_work() {
 		0, 12, 0, 13, 0, 0, 0, 14, 0, 15,
 		0, 0, 0, 16, 0, 0, 0, 0, 0, 17,
 		18, 0 }; // last entry 0 or 1 makes no performance difference
-
-	private long LCGSTATE;
 
-	private TDiv tdiv = new TDiv();
+	private TDiv63Inverse tdiv = new TDiv63Inverse(1<<21);
 
 	private BPSWTest bpsw = new BPSWTest();
 
 	private Gcd63 gcd63 = new Gcd63();
 
 	private SpRand32 spRand32 = new SpRand32();
-
+
+	private boolean doTDivFirst;
+
+	/**
+	 * Standard constructor.<br/><br/>
+	 * <strong>WARNING: Without trial division, tinyEcm might fail (i.e. run forever) for some N having >=2 small factors, like 735 = 3*5*7^2.</strong>
+     *
+	 * @param doTDivFirst
+	 */
+	public TinyEcm64(boolean doTDivFirst) {
+		this.doTDivFirst = doTDivFirst;
+	}
+
 	public String getName() {
-		return "TinyEcm64";
+		String tdivStr = doTDivFirst ? "with tdiv" : "without tdiv";
+		return "TinyEcm64(" + tdivStr + ")";
 	}
 
 	/**
@@ -1086,28 +1100,25 @@ public static long montMul64(long a, long b, long N, long Nhat) {
      * @return factor or 0 if no factor
 	 */
 	long check_factor(long Z, long n) {
-		long f = gcd63.gcd(Z,  n);
+		long f = gcd63.gcd(Z, n);
 		if (DEBUG) LOG.debug("check_factor: gcd(" + Z + ", " + n + ") = " + f);
         return (f>1 && f<n) ? f : 0;
 	}
 
-	/**
-	 * {@inheritDoc}
-	 * 
-	 * <br><br>
-	 * WARNING: tinyEcm's findSingleFactor() implementation is unstable when called with N having only small factors, like 735=3*5*7^2.
-	 * In such cases, an suitable amount of trial division should be carried out before.
-	 */
 	@Override
 	public BigInteger findSingleFactor(BigInteger N) {
-		// original rng, not comparable with C version
-		//Random rng = new Random();
-		//rng.setSeed(42);
-		//LCGSTATE = 65537 * rng.nextInt();
-		// rng comparable with C version
-		LCGSTATE = 4295098403L;
-		if (DEBUG) LOG.debug("LCGSTATE = " + LCGSTATE);
-
+		return findSingleFactor(N, this.doTDivFirst);
+	}
+
+	private BigInteger findSingleFactor(BigInteger N, boolean doTDivFirst) {
+		if (doTDivFirst) {
+			// Do trial division before ecm.
+			// The required amount of trial division has been derived experimentally.
+			// With a tdiv limit of 2^8 it still fails sometimes; 2^9 is stable; 2^10 is best in terms of performance.
+			final long factor = tdiv.setTestLimit(1<<10).findSingleFactor(N.longValue());
+			if (factor > 1) return BigInteger.valueOf(factor);
+		}
+
 		int NBits = N.bitLength();
 		if (NBits > MAX_BITS_SUPPORTED) throw new IllegalArgumentException("N=" + N + " has " + NBits + " bit, but tinyEcm only supports arguments <= " + MAX_BITS_SUPPORTED + " bit.");
 		// TODO Try to make it work for 63, 64 bit numbers
@@ -1142,27 +1153,35 @@ public BigInteger findSingleFactor(BigInteger N) {
 
 	@Override
 	public void searchFactors(FactorArguments args, FactorResult result) {
-		// tinyEcm is unstable when it comes to find factors of small composites like 735 = 3*5*7^2
-		// so we need some trial division first
-		tdiv.setTestLimit(1<<16).searchFactors(args, result);
-
-		if (result.untestedFactors.isEmpty()) return; // N was "easy"
-		// Otherwise we continue
-		BigInteger N = result.untestedFactors.firstKey();
-		int exp = result.untestedFactors.removeAll(N);
-		if (DEBUG) Ensure.ensureEquals(1, exp); // looks safe, otherwise we'ld have to consider exp below
-
-		if (bpsw.isProbablePrime(N)) {
-			result.primeFactors.add(N);
-			return;
+		BigInteger N;
+		int exp;
+		if (this.doTDivFirst) {
+			// Do trial division before ecm.
+			// The required amount of trial division has been derived experimentally.
+			// With a tdiv limit of 2^8 it still fails sometimes; 2^9 is stable; 2^10 is best in terms of performance.
+			tdiv.setTestLimit(1<<10).searchFactors(args, result);
+			if (DEBUG) LOG.debug("result after tdiv = " + result);
+
+			if (result.untestedFactors.isEmpty()) return; // N was "easy"
+			// Otherwise we continue
+			N = result.untestedFactors.firstKey();
+			exp = result.untestedFactors.removeAll(N); // can be > 1
+
+			if (bpsw.isProbablePrime(N)) {
+				result.primeFactors.add(N, exp);
+				return;
+			}
+		} else {
+			N = args.N;
+			exp = args.exp;
 		}
-
-		BigInteger factor1 = findSingleFactor(N);
-		//LOG.debug("After findSingleFactor()...");
+		
+		BigInteger factor1 = findSingleFactor(N, false);
+		if (DEBUG) LOG.debug("result after findSingleFactor() = " + result);
 		if (factor1.compareTo(I_1) > 0 && factor1.compareTo(N) < 0) {
 			// We found a factor, but here we cannot know if it is prime or composite
-			result.untestedFactors.add(factor1, args.exp);
-			result.untestedFactors.add(N.divide(factor1), args.exp);
+			result.untestedFactors.add(factor1, exp);
+			result.untestedFactors.add(N.divide(factor1), exp);
 		}
 
 		// nothing found