next PollardRhoBrentMontgomery64 version working with 63-bit inputs

Author: TilmanNeumann

src/main/java/de/tilman_neumann/jml/base/Uint128.java

@@ -245,6 +245,29 @@ public static Uint128 spMul64_MH(long a, long b) {
 
 		return new Uint128(r_hi, r_lo);
 	}
+	
+	/**
+	 * Multiplication of two signed 64-bit integers using Math.multiplyHigh().
+	 * Pretty fast if supported by intrinsics, which needs newer hardware and Java 10+.<br><br>
+	 * 
+	 * @param a
+	 * @param b
+	 * @return
+	 */
+	// TODO the method doesn't belong to this class anymore (because it's signed); several method names and coments need a workover
+	public static Uint128 mul64_MH_signed(long a, long b) {
+		final long r_lo = a*b;
+		long r_hi = Math.multiplyHigh(a, b);
+		
+		if (DEBUG) {
+			// compare to pure Java implementation
+			Uint128 testResult = mul64(a, b); // TODO compare with signed implementation
+			Ensure.ensureEquals(testResult.high, r_hi);
+			Ensure.ensureEquals(testResult.low, r_lo);
+		}
+
+		return new Uint128(r_hi, r_lo);
+	}
 
 	/**
 	 * The square of an unsigned 64 bit integer.

src/main/java/de/tilman_neumann/jml/factor/pollardRho/PollardRhoBrentMontgomery64MH.java

@@ -68,7 +68,7 @@ public String getName() {
 
 	@Override
 	public BigInteger findSingleFactor(BigInteger N) {
-		// there is a complication with this check: The algorithm works for some 63-bit numbers and there are tests for it, but there may be 63-bit numbers where it fails
+		// this version works for all 63 bit numbers!
 		if (N.bitLength() > 63) { // this check should be negligible in terms of performance
 			throw new IllegalArgumentException("N = " + N + " has " + N.bitLength() + " bit, but " + getName() + " only supports arguments <= 63 bit");
 		}
@@ -108,8 +108,7 @@ public long findSingleFactor(long nOriginal) {
 	    	        final int iMax = Math.min(m, r-k);
 	    	        for (int i=iMax; i>0; i--) {
 	    	            y = montMul64(y, y+1, n, minusNInvModR);
-	    	            final long diff = x<y ? y-x : x-y; // XXX would be nice if we could get rid of this like in PollardRhoBrentMontgomery32
-	    	            q = montMul64(diff, q, n, minusNInvModR);
+	    	            q = montMul64(y-x, q, n, minusNInvModR);
 	    	        }
 	    	        G = gcd.gcd(q, n);
 	    	        // if q==0 then G==n -> the loop will be left and restarted with new y
@@ -122,8 +121,7 @@ public long findSingleFactor(long nOriginal) {
 	    	if (G==n) {
 	    	    do {
 	    	        ys = montMul64(ys, ys+1, n, minusNInvModR);
-    	            final long diff = x<ys ? ys-x : x-ys;
-	    	        G = gcd.gcd(diff, n);
+	    	        G = gcd.gcd(ys-x, n);
 	    	    } while (G==1);
 	    	    if (DEBUG) LOG.debug("G = " + G);
 	    	}
@@ -170,7 +168,7 @@ private void setUpMontgomeryMult() {
 	 */
 	public static long montMul64(long a, long b, long N, long Nhat) {
 		// Step 1: Compute a*b
-		Uint128 ab = Uint128.spMul64_MH(a, b);
+		Uint128 ab = Uint128.mul64_MH_signed(a, b);
 
 		// Step 2: Compute t = ab * (-1/N) mod R
 		// Since R=2^64, "x mod R" just means to get the low part of x.
@@ -180,7 +178,7 @@ public static long montMul64(long a, long b, long N, long Nhat) {
 
 		// Step 3: Compute r = (a*b + t*N) / R
 		// Since R=2^64, "x / R" just means to get the high part of x.
-		long r = ab.add_getHigh(Uint128.spMul64_MH(t, N));
+		long r = ab.add_getHigh(Uint128.mul64_MH_signed(t, N));
 		// If the correct result is c, then now r==c or r==c+N.
 		// This is fine for this factoring algorithm, because r will 
 		// * either be subjected to another Montgomery multiplication mod N,