Skip to content

Commit a44067f

Browse files
TilmanNeumannTilmanNeumann
committed
faster prime test for small N using Lemire division
1 parent a2d92ee commit a44067f

8 files changed

Lines changed: 381 additions & 41 deletions

File tree

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

Lines changed: 9 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -37,18 +37,17 @@ public class LemireIntTrialDivision extends FactorAlgorithm {
3737
limits = new int[primes.length];
3838
for (int i = 0; i < primes.length; i++) {
3939
int prime = primes[i];
40-
// Modulare Inverse für 32-bit (Newton-Verfahren)
40+
// compute modular inverses of p (mod 2^32) using Newton's method
4141
int inv = modularInverseInt(prime);
4242
modularInverse[i] = inv;
43-
// Limit = (2^32 - 1) / prime (Unsigned)
44-
// In Java: Integer.divideUnsigned(-1, prime)
43+
// limit = (2^32 - 1) / prime (unsigned)
4544
limits[i] = Integer.divideUnsigned(-1, prime);
4645
}
4746
}
4847

4948
private static int modularInverseInt(int n) {
50-
int inverse = n;
51-
for (int i = 0; i < 4; i++) { // 4 Iterationen reichen für 32-bit
49+
int inverse = n; // initial estimate
50+
for (int i = 0; i < 4; i++) { // 4 iterations are sufficient for 32 bit numbers
5251
inverse *= 2 - n * inverse;
5352
}
5453
return inverse;
@@ -65,21 +64,21 @@ public BigInteger findSingleFactor(BigInteger N) {
6564
return BigInteger.valueOf(findSingleFactor(N.longValue()));
6665
}
6766

68-
public int findSingleFactor(long numberToFactorize) {
67+
public int findSingleFactor(long N) {
6968
// Lemire can not handle even numbers
70-
if ((numberToFactorize & 1) == 0) return 2;
69+
if ((N & 1) == 0) return 2;
7170

7271
for (int i = 1; i < primes.length; i++) {
7372
// for hard numbers like big semiprimes finding a factor (early) is unlikely and JIT predicts that
7473
// the return branch is unlikely -> always the same data processing; preloading the arrays
75-
if (factorFound (numberToFactorize, i)) return primes[i];
74+
if (factorFound (N, i)) return primes[i];
7675
}
7776

7877
return -1;
7978
}
8079

81-
private boolean factorFound(long numberToFactorize, int i) {
82-
int nInt = (int) numberToFactorize;
80+
private boolean factorFound(long N, int i) {
81+
int nInt = (int) N;
8382
int product = nInt * modularInverse[i];
8483
return Integer.compareUnsigned (product, limits[i]) <= 0;
8584
}

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

Lines changed: 12 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -38,17 +38,17 @@ public class LemireTrialDivision extends FactorAlgorithm {
3838
limits = new long[primes.length];
3939
for (int i = 0; i < primes.length; i++) {
4040
long p = primes[i];
41-
// Berechne die modulare Inverse für ungerade p (Newton-Verfahren)
41+
// compute modular inverses of p (mod 2^64) using Newton's method
4242
long inverse = modularInverse(p);
4343
modularInverse[i] = inverse;
44-
// Limit = (2^64 - 1) / p (Unsigned)
44+
// limit = (2^64 - 1) / prime (unsigned)
4545
limits[i] = Long.divideUnsigned(-1L, p);
4646
}
4747
}
4848

4949
private static long modularInverse(long n) {
50-
long inverse = n; // Initialer Schätzwert
51-
for (int i = 0; i < 5; i++) { // 5 Iterationen reichen für 64-bit
50+
long inverse = n; // initial estimate
51+
for (int i = 0; i < 5; i++) { // 5 iterations are sufficient for 64 bit numbers
5252
inverse *= 2 - n * inverse;
5353
}
5454
return inverse;
@@ -65,30 +65,28 @@ public BigInteger findSingleFactor(BigInteger N) {
6565
return BigInteger.valueOf(findSingleFactor(N.longValue()));
6666
}
6767

68-
public int findSingleFactor(long numberToFactorize) {
68+
public int findSingleFactor(long N) {
6969
// Lemire can not handle even numbers
70-
if ((numberToFactorize & 1) == 0) return 2;
70+
if ((N & 1) == 0) return 2;
7171

7272
for (int i = 1; i < primes.length; i++) {
7373
// for hard numbers like big semiprimes finding a factor (early) is unlikely and JIT predicts that
7474
// the return branch is unlikely -> always the same data processing; preloading the arrays
75-
if (factorFound (numberToFactorize, i)) return primes[i];
75+
if (factorFound (N, i)) return primes[i];
7676
}
7777

7878
return -1;
7979
}
8080

81-
private boolean factorFound(long numberToFactorize, int i) {
82-
// 1. Hole vorberechnete Inverse und Limit
81+
private boolean factorFound(long N, int i) {
82+
// 1. get pre-computed inverse and limit
8383
long inv = modularInverse[i];
8484
long limit = limits[i];
8585

86-
// 2. Multipliziere number * inverse (Überlauf ist beabsichtigt!)
87-
long product = numberToFactorize * inv;
86+
// 2. multiply number * inverse (overflow is intended!)
87+
long product = N * inv;
8888

89-
// 3. Wenn das Produkt (unsigned) kleiner oder gleich dem Limit ist,
90-
// dann ist numberToFactorize restlos durch primes[i] teilbar.
91-
// the calculation is done completely in long -> might be the main speedup 50%
89+
// 3. if the (unsigned) product is less than or equal to the limit, then primes[i] divides N without rest.
9290
return Long.compareUnsigned(product, limit) <= 0;
9391
}
9492
}

src/main/java/de/tilman_neumann/jml/primes/probable/TDivPrimeTest.java renamed to src/main/java/de/tilman_neumann/jml/primes/exact/BarrettPrimeTest.java

Lines changed: 49 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*
22
* 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.
3-
* Copyright (C) 2018-2024 Tilman Neumann - tilman.neumann@web.de
3+
* Copyright (C) 2018-2026 Tilman Neumann - tilman.neumann@web.de
44
*
55
* This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License
66
* as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
@@ -11,41 +11,47 @@
1111
* You should have received a copy of the GNU General Public License along with this program;
1212
* if not, see <http://www.gnu.org/licenses/>.
1313
*/
14-
package de.tilman_neumann.jml.primes.probable;
14+
package de.tilman_neumann.jml.primes.exact;
1515

1616
import org.apache.logging.log4j.Logger;
17-
import org.apache.logging.log4j.LogManager;
1817

19-
import de.tilman_neumann.jml.primes.exact.AutoExpandingPrimesArray;
18+
import de.tilman_neumann.jml.BinarySearch;
19+
20+
import org.apache.logging.log4j.LogManager;
2021

2122
/**
22-
* A deterministic prime test for N < 32 bit using fast trial division.
23+
* A deterministic prime test for N < 32 bit using long-valued Barrett trial division.
24+
*
25+
* Implements the singleton pattern so that the resources will not be allocated twice no
26+
* matter how often the class is used.
2327
*
2428
* @author Tilman Neumann
2529
*/
26-
public class TDivPrimeTest {
30+
public class BarrettPrimeTest {
2731
@SuppressWarnings("unused")
28-
private static final Logger LOG = LogManager.getLogger(TDivPrimeTest.class);
32+
private static final Logger LOG = LogManager.getLogger(BarrettPrimeTest.class);
2933

3034
private static final int NUM_PRIMES_FOR_31_BIT_TDIV = 4793;
3135

36+
private static final BinarySearch binarySeach = new BinarySearch();
37+
3238
private int[] primes;
3339
private long[] pinv;
3440

35-
// lazy-initialized singleton
36-
private static TDivPrimeTest the_instance = null;
41+
/** lazy-initialized singleton */
42+
private static BarrettPrimeTest the_instance = null;
3743

3844
/**
3945
* @return the only TDivPrimeTest instance (singleton)
4046
*/
41-
public static synchronized final TDivPrimeTest getInstance() {
47+
public static synchronized final BarrettPrimeTest getInstance() {
4248
if (the_instance == null) {
43-
the_instance = new TDivPrimeTest();
49+
the_instance = new BarrettPrimeTest();
4450
}
4551
return the_instance;
4652
}
4753

48-
private TDivPrimeTest() {
54+
private BarrettPrimeTest() {
4955
AutoExpandingPrimesArray smallPrimesProvider = AutoExpandingPrimesArray.get();
5056
primes = new int[NUM_PRIMES_FOR_31_BIT_TDIV];
5157
pinv = new long[NUM_PRIMES_FOR_31_BIT_TDIV];
@@ -56,7 +62,37 @@ private TDivPrimeTest() {
5662
}
5763
}
5864

59-
public boolean isPrime(int N) {
65+
// not public because Lemire is faster
66+
boolean isPrime_v1(int N) {
67+
if (N==1) return false;
68+
if ((N&1)==0) return N==2;
69+
70+
// if N is odd and composite then the loop runs maximally up to prime = floor(sqrt(N))
71+
int i=1;
72+
for ( ; i<NUM_PRIMES_FOR_31_BIT_TDIV; i++) {
73+
if ((1 + (int) ((N*pinv[i])>>32)) * primes[i] == N) return primes[i]==N;
74+
}
75+
// N is prime
76+
return true;
77+
}
78+
79+
// not public because Lemire is faster
80+
boolean isPrime_v2(int N) {
81+
if (N==1) return false;
82+
if ((N&1)==0) return N==2;
83+
84+
int pmax = (int) Math.sqrt(N);
85+
int imax = binarySeach.getInsertPosition(primes, NUM_PRIMES_FOR_31_BIT_TDIV-1, pmax);
86+
87+
for (int i = 1; i<=imax; i++) {
88+
if ((1 + (int) ((N*pinv[i])>>32)) * primes[i] == N) return primes[i]==N;
89+
}
90+
// N is prime
91+
return true;
92+
}
93+
94+
// not public because Lemire is faster
95+
boolean isPrimeUnrolled(int N) {
6096
if (N==1) return false;
6197
if ((N&1)==0) return N==2;
6298

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

Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
/*
2+
* 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.
3+
* Copyright (C) 2018-2026 Tilman Neumann - tilman.neumann@web.de
4+
*
5+
* This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License
6+
* as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
7+
*
8+
* This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
9+
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
10+
*
11+
* You should have received a copy of the GNU General Public License along with this program;
12+
* if not, see <http://www.gnu.org/licenses/>.
13+
*/
14+
package de.tilman_neumann.jml.primes.exact;
15+
16+
import org.apache.logging.log4j.Logger;
17+
import org.apache.logging.log4j.LogManager;
18+
19+
import de.tilman_neumann.jml.BinarySearch;
20+
21+
/**
22+
* A deterministic prime test for N < 32 bit using Lemire division.
23+
*
24+
* Implements the singleton pattern so that the resources will not be allocated twice no
25+
* matter how often the class is used.
26+
*
27+
* @author Tilman Neumann
28+
*/
29+
public class LemirePrimeTest {
30+
@SuppressWarnings("unused")
31+
private static final Logger LOG = LogManager.getLogger(LemirePrimeTest.class);
32+
33+
private static final int MAX_N_BITS = 36;
34+
private static final int MAX_PRIME_FACTOR = 1 << ((MAX_N_BITS+1)/2);
35+
36+
private static final BinarySearch binarySeach = new BinarySearch();
37+
38+
private static int MAX_INDEX; // for N<32 bit this would be 4792
39+
40+
private int[] primes;
41+
private long[] modularInverse;
42+
private long[] limits;
43+
44+
// lazy-initialized singleton
45+
private static LemirePrimeTest the_instance = null;
46+
47+
/**
48+
* @return the only TDivPrimeTest instance (singleton)
49+
*/
50+
public static synchronized final LemirePrimeTest getInstance() {
51+
if (the_instance == null) {
52+
the_instance = new LemirePrimeTest();
53+
}
54+
return the_instance;
55+
}
56+
57+
private LemirePrimeTest() {
58+
primes = SmallPrimes.generatePrimes(MAX_PRIME_FACTOR);
59+
MAX_INDEX = primes.length - 1;
60+
modularInverse = new long[primes.length];
61+
limits = new long[primes.length];
62+
for (int i = 0; i < primes.length; i++) {
63+
long p = primes[i];
64+
// compute modular inverses of p (mod 2^64) using Newton's method
65+
long inverse = modularInverse(p);
66+
modularInverse[i] = inverse;
67+
// limit = (2^64 - 1) / p (unsigned)
68+
limits[i] = Long.divideUnsigned(-1L, p);
69+
}
70+
}
71+
72+
private static long modularInverse(long n) {
73+
long inverse = n; // initial estimate
74+
for (int i = 0; i < 5; i++) { // 5 iterations are sufficient for 64 bit numbers
75+
inverse *= 2 - n * inverse;
76+
}
77+
return inverse;
78+
}
79+
80+
// not public because isPrimeUnrolled() is faster
81+
boolean isPrime_v1(int N) {
82+
if (N==1) return false;
83+
if ((N&1)==0) return N==2;
84+
85+
int pmax = (int) Math.sqrt(N);
86+
87+
for (int i = 1; primes[i]<=pmax; i++) {
88+
// for hard numbers like big semiprimes finding a factor (early) is unlikely and JIT predicts that
89+
// the return branch is unlikely -> always the same data processing; preloading the arrays
90+
if (factorFound (N, i)) return false;
91+
}
92+
93+
return true;
94+
}
95+
96+
// not public because isPrimeUnrolled() is faster
97+
boolean isPrime_v2(int N) {
98+
if (N==1) return false;
99+
if ((N&1)==0) return N==2;
100+
101+
int pmax = (int) Math.sqrt(N);
102+
int imax = binarySeach.getInsertPosition(primes, MAX_INDEX, pmax);
103+
104+
for (int i = 1; i<=imax; i++) {
105+
// for hard numbers like big semiprimes finding a factor (early) is unlikely and JIT predicts that
106+
// the return branch is unlikely -> always the same data processing; preloading the arrays
107+
if (factorFound (N, i)) return false;
108+
}
109+
110+
return true;
111+
}
112+
113+
public boolean isPrime/*Unrolled*/(int N) {
114+
if (N==1) return false;
115+
if ((N&1)==0) return N==2;
116+
117+
int pmax = (int) Math.sqrt(N);
118+
int imax = binarySeach.getInsertPosition(primes, MAX_INDEX, pmax);
119+
120+
int i=1;
121+
int unrolledLimit = imax-8;
122+
for ( ; i<unrolledLimit; i++) {
123+
if (factorFound (N, i)) return false;
124+
if (factorFound (N, ++i)) return false;
125+
if (factorFound (N, ++i)) return false;
126+
if (factorFound (N, ++i)) return false;
127+
if (factorFound (N, ++i)) return false;
128+
if (factorFound (N, ++i)) return false;
129+
if (factorFound (N, ++i)) return false;
130+
if (factorFound (N, ++i)) return false;
131+
}
132+
for ( ; i<imax; i++) {
133+
if (factorFound (N, i)) return false;
134+
}
135+
136+
return true;
137+
}
138+
139+
private boolean factorFound(long N, int i) {
140+
// 1. get pre-computed inverse and limit
141+
long inv = modularInverse[i];
142+
long limit = limits[i];
143+
144+
// 2. multiply number * inverse (overflow is intended!)
145+
long product = N * inv;
146+
147+
// 3. if the (unsigned) product is less than or equal to the limit, then primes[i] divides N without rest.
148+
return Long.compareUnsigned(product, limit) <= 0;
149+
}
150+
}

0 commit comments

Comments
 (0)