Skip to content

Commit 66b0685

Browse files
TilmanNeumannTilmanNeumann
committed
add Lemire division implementations thanks to Thilo Harich
1 parent d46f815 commit 66b0685

4 files changed

Lines changed: 204 additions & 1 deletion

File tree

‎src/main/java/de/tilman_neumann/jml/factor/tdiv/LemireIntTrialDivision.java‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
package de.tilman_neumann.jml.factor.tdiv;
2+
3+
import java.math.BigInteger;
4+
5+
import de.tilman_neumann.jml.factor.FactorAlgorithm;
6+
import de.tilman_neumann.jml.primes.exact.SmallPrimes;
7+
8+
/**
9+
* Lemire division for int arguments.
10+
*
11+
* @author Thilo Harich
12+
*/
13+
public class LemireIntTrialDivision extends FactorAlgorithm {
14+
15+
private static final int MAX_PRIME_FACTOR = 1<<16; // sufficient for numbers to factor with 32 bits
16+
17+
private static int[] primes;
18+
private static int[] modularInverse;
19+
private static int[] limits;
20+
21+
static {
22+
primes = SmallPrimes.generatePrimes(MAX_PRIME_FACTOR);
23+
modularInverse = new int[primes.length];
24+
limits = new int[primes.length];
25+
for (int i = 0; i < primes.length; i++) {
26+
int prime = primes[i];
27+
// Modulare Inverse für 32-bit (Newton-Verfahren)
28+
int inv = modularInverseInt(prime);
29+
modularInverse[i] = inv;
30+
// Limit = (2^32 - 1) / prime (Unsigned)
31+
// In Java: Integer.divideUnsigned(-1, prime)
32+
limits[i] = Integer.divideUnsigned(-1, prime);
33+
}
34+
}
35+
36+
private static int modularInverseInt(int n) {
37+
int inverse = n;
38+
for (int i = 0; i < 4; i++) { // 4 Iterationen reichen für 32-bit
39+
inverse *= 2 - n * inverse;
40+
}
41+
return inverse;
42+
}
43+
44+
@Override
45+
public String getName() {
46+
return "LemireIntTrialDivision";
47+
}
48+
49+
@Override
50+
public BigInteger findSingleFactor(BigInteger N) {
51+
if (N.bitLength() > 32) throw new IllegalArgumentException("LemireIntTrialDivision.findSingleFactor() does not work for N>32 bit, but N=" + N + " has " + N.bitLength() + " bits.");
52+
return BigInteger.valueOf(findSingleFactor(N.longValue()));
53+
}
54+
55+
public int findSingleFactor(long numberToFactorize) {
56+
// Lemire can not handle even numbers
57+
if ((numberToFactorize & 1) == 0) return 2;
58+
59+
for (int i = 1; i < primes.length; i++) {
60+
// for hard numbers like big semiprimes finding a factor (early) is unlikely and JIT predicts that
61+
// the return branch is unlikely -> always the same data processing; preloading the arrays
62+
if (factorFound (numberToFactorize, i)) return primes[i];
63+
}
64+
65+
return -1;
66+
}
67+
68+
private boolean factorFound(long numberToFactorize, int i) {
69+
int nInt = (int) numberToFactorize;
70+
int product = nInt * modularInverse[i];
71+
return Integer.compareUnsigned (product, limits[i]) <= 0;
72+
}
73+
}

‎src/main/java/de/tilman_neumann/jml/factor/tdiv/LemireTrialDivision.java‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
package de.tilman_neumann.jml.factor.tdiv;
2+
3+
import java.math.BigInteger;
4+
5+
import de.tilman_neumann.jml.factor.FactorAlgorithm;
6+
import de.tilman_neumann.jml.primes.exact.SmallPrimes;
7+
8+
/**
9+
* Lemire division for long arguments.
10+
*
11+
* @author Thilo Harich
12+
*/
13+
public class LemireTrialDivision extends FactorAlgorithm {
14+
15+
private static final int MAX_N_BITS = 36;
16+
private static final int MAX_PRIME_FACTOR = 1 << ((MAX_N_BITS+1)/2);
17+
18+
private static int[] primes;
19+
private static long[] modularInverse;
20+
private static long[] limits;
21+
22+
static {
23+
primes = SmallPrimes.generatePrimes(MAX_PRIME_FACTOR);
24+
modularInverse = new long[primes.length];
25+
limits = new long[primes.length];
26+
for (int i = 0; i < primes.length; i++) {
27+
long p = primes[i];
28+
// Berechne die modulare Inverse für ungerade p (Newton-Verfahren)
29+
long inverse = modularInverse(p);
30+
modularInverse[i] = inverse;
31+
// Limit = (2^64 - 1) / p (Unsigned)
32+
limits[i] = Long.divideUnsigned(-1L, p);
33+
}
34+
}
35+
36+
private static long modularInverse(long n) {
37+
long inverse = n; // Initialer Schätzwert
38+
for (int i = 0; i < 5; i++) { // 5 Iterationen reichen für 64-bit
39+
inverse *= 2 - n * inverse;
40+
}
41+
return inverse;
42+
}
43+
44+
@Override
45+
public String getName() {
46+
return "LemireTrialDivision";
47+
}
48+
49+
@Override
50+
public BigInteger findSingleFactor(BigInteger N) {
51+
if (N.bitLength() > MAX_N_BITS) throw new IllegalArgumentException("LemireTrialDivision.findSingleFactor() does not work for N>" + MAX_N_BITS + " bit, but N=" + N + " has " + N.bitLength() + " bits.");
52+
return BigInteger.valueOf(findSingleFactor(N.longValue()));
53+
}
54+
55+
public int findSingleFactor(long numberToFactorize) {
56+
// Lemire can not handle even numbers
57+
if ((numberToFactorize & 1) == 0) return 2;
58+
59+
for (int i = 1; i < primes.length; i++) {
60+
// for hard numbers like big semiprimes finding a factor (early) is unlikely and JIT predicts that
61+
// the return branch is unlikely -> always the same data processing; preloading the arrays
62+
if (factorFound (numberToFactorize, i)) return primes[i];
63+
}
64+
65+
return -1;
66+
}
67+
68+
private boolean factorFound(long numberToFactorize, int i) {
69+
// 1. Hole vorberechnete Inverse und Limit
70+
long inv = modularInverse[i];
71+
long limit = limits[i];
72+
73+
// 2. Multipliziere number * inverse (Überlauf ist beabsichtigt!)
74+
long product = numberToFactorize * inv;
75+
76+
// 3. Wenn das Produkt (unsigned) kleiner oder gleich dem Limit ist,
77+
// dann ist numberToFactorize restlos durch primes[i] teilbar.
78+
// the calculation is done completely in long -> might be the main speedup 50%
79+
return Long.compareUnsigned(product, limit) <= 0;
80+
}
81+
}
82+

‎src/main/java/de/tilman_neumann/jml/primes/exact/SmallPrimes.java‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
package de.tilman_neumann.jml.primes.exact;
2+
3+
/**
4+
* Generation of primes for Lemire division.
5+
*
6+
* @author Thilo Harich
7+
*/
8+
public class SmallPrimes {
9+
10+
public static int[] generatePrimes(int limit) {
11+
boolean[] isComposite = getIsComposite(limit);
12+
13+
int count = getCount(limit, isComposite);
14+
15+
// Array mit Primzahlen füllen
16+
int[] primes = new int[count];
17+
int idx = 0;
18+
for (int i = 2; i <= limit; i++) {
19+
if (!isComposite[i]) primes[idx++] = i;
20+
}
21+
22+
return primes;
23+
}
24+
25+
private static int getCount(int limit, boolean[] isComposite) {
26+
int count = 0;
27+
for (int i = 2; i <= limit; i++) {
28+
if (!isComposite[i]) count++;
29+
}
30+
return count;
31+
}
32+
33+
private static boolean[] getIsComposite(int limit) {
34+
boolean[] isComposite = new boolean[limit + 1];
35+
for (int i = 2; i * i <= limit; i++) {
36+
if (!isComposite[i]) {
37+
for (int j = i * i; j <= limit; j += i) {
38+
isComposite[j] = true;
39+
}
40+
}
41+
}
42+
return isComposite;
43+
}
44+
}

‎src/test/java/de/tilman_neumann/jml/factor/FactorizerTest.java‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -111,7 +111,9 @@ public FactorizerTest() {
111111
//new TDiv63(),
112112
// new TDiv63Inverse(1<<21),
113113
//new TDiv().setTestLimit(1<<21),
114-
114+
// new LemireIntTrialDivision(),
115+
// new LemireTrialDivision(),
116+
115117
// Hart's one line factorizer
116118
//new HartSimple(),
117119
// new HartFast(true),
@@ -243,6 +245,8 @@ private void testRange(int bits) {
243245
if (bits>31 && algName.startsWith("PollardRhoBrent31")) continue; // int implementation
244246
if (bits>31 && algName.startsWith("PollardRhoTwoLoops31")) continue; // int implementation
245247
if (bits>31 && algName.startsWith("PollardRhoBrentMontgomery32")) continue; // int implementation
248+
if (bits>32 && algName.startsWith("LemireIntTrialDivision")) continue; // int implementation
249+
if (bits>36 && algName.startsWith("LemireTrialDivision")) continue; // not enough primes stored
246250
if (bits>42 && algName.startsWith("TDiv63Inverse")) continue; // not enough primes stored
247251
if (bits>52 && algName.startsWith("SquFoF31")) continue; // int implementation
248252
if (bits>59 && algName.startsWith("Lehman")) continue;

0 commit comments

Comments
 (0)