# Finding number of primes less than 10 million using a multithreaded program

I recently decided to learn how to multithread Java programs, so I made a small program to compare the performance of serial and multithreaded programs that perform the same task.

I created a serial program that calculates the number of primes from 1 to 10 million, and timed it 50 times using a test program. Here's the code for the serial program:

import java.util.Locale;
import java.text.NumberFormat;
/**
* A serial program that calculates the number of primes less than a
* certain number (which is hard coded). Used as a basis for
* benchmarking the multi-threaded programs that do the same thing.
*
* @author Tirthankar Mazumder
* @version 1.2
* @date 2nd May, 2021
*/
public class PrimalityTestSerial {
public static final long max = 10_000_000;
public static void main(String[] args) {
final long startTime = System.currentTimeMillis();

long num_primes = primeCalculator();

final long endTime = System.currentTimeMillis();

NumberFormat nf = NumberFormat.getInstance(Locale.US);

System.out.println("Number of primes less than " + nf.format(max) + ": " + num_primes);
System.out.println("Took " + (endTime - startTime) + " ms.");
System.out.println();
}

private static boolean isPrime(long l) {
long upper_bound = (long) Math.floor(Math.sqrt(l));

for (long i = 2; i <= upper_bound; i++) {
if (l % i == 0)
return false;
}

return true;
}

public static long primeCalculator() {
long num_primes = 0;

for (long i = 2; i <= max; i++) {
if (isPrime(i))
num_primes++;
}
return num_primes;
}
}


Here's the code for the worker class used in the multithreaded version:

/**
* A worker class for calculating the number of primes from start to end,
* which are private member variables. Instances of this class are used
* in the multithreaded version of PrimalityTestSerial.
*
* @author Tirthankar Mazumder
* @version 1.2
* @date 3rd May, 2021
*/
public class PrimalityTestWorker1 implements Runnable {
//Member variables
public static long totalPrimeCount = 0;
private final long start;
private final long end;

public PrimalityTestWorker1(long start, long end) {
this.start = start;
this.end = end;
}

private synchronized void increment(long num) {
totalPrimeCount += num;
}

private static boolean isPrime(long l) {
long upper_bound = (long) Math.floor(Math.sqrt(l));

for (long i = 2; i <= upper_bound; i++) {
if (l % i == 0)
return false;
}

return true;
}

private void numPrimes() {
long primeCount = 0;
for (long i = start; i <= end; i++) {
if (isPrime(i))
primeCount++;
}
increment(primeCount);
}

public void run() {
numPrimes();
}
}


Here's the main program which uses instances of PrimalityTest1Worker as threads:

import java.util.Locale;
import java.text.NumberFormat;
/**
* The master program for the multithreaded primality test that creates
* objects of the PrimalityTestWorker1 to make threads, and then collates
* the results and prints them to stdout.
*
* @author Tirthankar Mazumder
* @version 1.0=2
* @date 3rd May, 2021
*/
public class PrimalityTestParallel1Runner {
public static final int cores = Runtime.getRuntime().availableProcessors();
//We will spawn as many threads as there are cores on the system, and not
//more than that because we are not I/O bound here.

public static final long max = PrimalityTestSerial.max;
//For consistency.

public static void main(String[] args) {
long startTime = System.currentTimeMillis();

primeCalculator();

long endTime = System.currentTimeMillis();

NumberFormat nf = NumberFormat.getInstance(Locale.US);

System.out.println("Number of primes less than " + nf.format(max) + ": " +
PrimalityTestWorker1.totalPrimeCount);

System.out.println("Took " + (endTime - startTime) + " ms.");
System.out.println();
}

public static void primeCalculator() {

long chunk = max / cores;

for (int i = 0; i < cores; i++) {
t.start();

}

for (int i = 0; i < cores; i++) {
try {
} catch (InterruptedException e) {
System.out.println("Was interrupted.");
return;
}
}
}
}


Finally, here's the code for the testing program, which runs each program 50 times and then calculates the average runtimes:

import java.util.Arrays;
/**
* A wrapper class that handles benchmarking the performances of
* PrimalityTestSerial and PrimalityTestParallel1Runner and then
* prints information about the results to stdout.
*
* @author Tirthankar Mazumder
* @version 1.0
* @date 8th May, 2021
*/
public class PrimalityTestSuite {
public static final int n = 50;
//Number of test runs to perform

public static void main(String[] args) {
long totalSerialTime = 0;
long totalParallelTime = 0;

long serialTimes[] = new long[n];

double avgSerialTime = 0;
double avgParallelTime = 0;

System.out.println("Starting Serial runs...");
long startTime = System.currentTimeMillis();
for (int i = 0; i < n; i++) {
PrimalityTestSerial.primeCalculator();
serialTimes[i] = System.currentTimeMillis();
}

for (int i = 0; i < n; i++) {
serialTimes[i] -= startTime;
for (int j = 0; j < i; j++) {
serialTimes[i] -= serialTimes[j];
//to get rid of the time taken by the previous runs
}
avgSerialTime += serialTimes[i];
}

avgSerialTime /= n;

long parallelTimes[] = new long[n];

System.out.println("Starting parallel runs...");
startTime = System.currentTimeMillis();
for (int i = 0; i < n; i++) {
PrimalityTestParallel1Runner.primeCalculator();
parallelTimes[i] = System.currentTimeMillis();
}

for (int i = 0; i < n; i++) {
parallelTimes[i] -= startTime;
for (int j = 0; j < i; j++) {
parallelTimes[i] -= parallelTimes[j];
//to get rid of the time taken by the previous runs
}
avgParallelTime += parallelTimes[i];
}

avgParallelTime /= n;

Arrays.sort(serialTimes);
Arrays.sort(parallelTimes);

double bestThreeSerialAvg = (serialTimes[0] + serialTimes[1]
+ serialTimes[2]) / 3;
double bestThreeParallelAvg = (parallelTimes[0] + parallelTimes[1]
+ parallelTimes[2]) / 3;

System.out.println();
System.out.println("Results:");

System.out.println("Average of " + n + " Serial Runs: " + avgSerialTime + " ms.");
System.out.println("Average of " + n + " Parallel Runs: " + avgParallelTime + " ms.");

System.out.println();
System.out.println("Average speed-up: " + avgSerialTime / avgParallelTime + "x");
System.out.println();

System.out.println("Average of best 3 Serial Runs: " + bestThreeSerialAvg + " ms.");
System.out.println("Average of best 3 Parallel Runs: " + bestThreeParallelAvg + " ms.");

System.out.println();
System.out.println("Average speed-up (w.r.t. best run times): "
+ bestThreeSerialAvg / bestThreeParallelAvg + "x");
System.out.println();
}
}


Here are the results from the test program:

Starting Serial runs...
Starting parallel runs...

Results:
Average of 50 Serial Runs: 4378.92 ms.
Average of 50 Parallel Runs: 1529.2 ms.

Average speed-up: 2.8635364896678x

Average of best 3 Serial Runs: 4328.0 ms.
Average of best 3 Parallel Runs: 1297.0 ms.

Average speed-up (w.r.t. best run times): 3.3369313801079414x


From here, it is obvious that the average speed-up is simply around 3x. However, this is surprising because I expect the multi-threaded program to run 7 to 8 times faster (because the serial program uses just 1 core, whereas the multi-threaded program should use all 8 cores on my system.)

So my question is, why is the multi-threaded program not as fast as I expect it to be?

• Every time I see a prime-counting function program I feel obligated to tell you there are much faster algorithms that run in sub-linear time using en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm My pure python implementation runs in 0.01 sec. – qwr May 9 at 6:10
• @qwr Thanks, I didn't know about this algorithm. But even if I did, I wouldn't have used it here, because my objective was to make the threads do some nontrivial work so that I can measure my performance gains. But for tasks where prime-counting is necessary and performance is required, I will most definitely use that algorithm from now on :D – wermos May 9 at 9:15
• @qwr: That's still way too much work. I use 1e7 / math.log(1e7) as an approximation and get a result in 7ns. ;) – Eric Duminil May 9 at 19:23
• Yes, the PNT provides a very good approximation – qwr May 9 at 19:47
• The dynamic programming algorithm for finding primes, known as the Sieve of Eratosthenes, is a very good exercise in principles of HPC if you want to up-level your skills. Parallel tasks must synchronize access to a common array of bits. You'll need a queue of tasks (each task is to strike every Nth number), then dispatch them to threads, decide when it's cheaper to lump tasks (every 2nd runs a long time, every 1E6th is fast). And how to optimize it: find first 1E4 primes then schedule tasks only for prime N<1E4 up to 1E6, rinse&repeat. Try 1E11, it needs only a 5G array, as primes>2 are odd. – kkm May 10 at 5:14

threadEnd = (threadEnd + chunk > max) ? max : threadEnd + chunk;


Your primeCalculator's equally sized chunks run in unequal time. The calls to primeCalculator for, say, 1-20 are going to be significantly faster than the calls to primeCalculator for, say, 10000-10020. When 1-20 is done 10000-10020 will still be running, hence performance isn't cut to an 1/8.

Following on from Peilonrayz' answer, to counteract that problem as much as possible, you can do multiple "chunks" per thread, making each such chunk smaller, and spreading them out across the values.

In the extreme, the chunks are of size 1. If you know how many threads there will be (we see you have cores threads), then this is fairly simple to implement: tell each thread t the total number of threads n and an initial value x to check for primality, then add n after each check. That is, instead of providing a threadStart and threadEnd value, provide a threadStart value between 1 and n, and a threadIncrement value (that is just n; alternatively, you could make that n a global variable that all threads can see so that you don't need to pass it into the worker constructor). In each thread, you then do something like this:

x = threadStart;

while (x < max) {
if (isPrime(x)) {
// ...
}
x += n;
}

• Note that every even number is composite (except 2), and will leave the division loop right away, on the first divisor. Dividing up the odd numbers between threads is a very good idea, though, with x += nthreads, but if you do this naively then half your threads will exit much sooner than the others. (Which is maybe ok on a CPU with 2 logical cores per physical core, otherwise you're wasting half your available cores.) – Peter Cordes May 9 at 19:17
• @PeterCordes Of course there are more efficient ways to do this, even much more efficient that what you have stated. – Jivan Pal May 9 at 19:43
• My point wasn't primarily about making it more efficient (that's a side benefit), but making the work done in each thread actually about equal. So you can use static scheduling with 8 equal chunks, and not have some threads still working long after others exit, assuming they all had a CPU core to run on the whole time. But that problem would happen if some threads are only checking even numbers, and some are only checking odd. (e.g. starting at 2 vs. starting at 3, and incrementing by 8. Every even number has n % i == 0 for i=2 so they exit on the first iteration.) – Peter Cordes May 9 at 19:55
• @PeterCordes, ah, yes, this is actually a good point. If the number of threads n is not prime, then those threads whose index is a factor of n will finish very quickly. – Jivan Pal May 9 at 20:01
• Oh yes, I was assuming an even number of threads, so you could add nthreads and always get odd numbers. Or even if you had a setup with 6 cores for example (e.g. Coffee Lake without hyperthreading), a thread that started with n=3 would get all multiples of 3. So yeah, there's a possible problem for any non-power-of-2 number of threads, I think. (It's not nthreads being non-prime that's the issue, e.g. nt=7 would make the thread starting with n=7 a problem). – Peter Cordes May 9 at 20:13

Many modern CPUs can run at a higher clock speed when only one core is active, vs. when all cores are active. This ratio between one-core vs. all-core max turbo / boost-clock depends entirely on the CPU model. (It can be higher in laptops, where thermal constraints are a bigger deal, for the similar reasons to max turbo vs. sustained frequency being a higher ratio in laptops.)

Also:

If your CPU has SMT (e.g. Intel Hyper Threading), one thread per physical core is probably mostly maxing out the integer-division throughput. (Composite numbers with small prime factors will mean you break out of the loop relatively often, letting SMT do its job by letting the other thread get more done while this thread is recovering from a branch misprediction.)

But still, with SMT, best-case speedup may only scale with number of physical cores, not the number you detect with availableProcessors. E.g. 4x on a 4c8t (4 core 8 thread) CPU like i7-6700k which the OS will detect as having 8 logical cores, but has 4 physical.

(Division is a very slow operation, and not fully pipelined even in modern high-end x86 CPUs. I.e. it can't start a new division every clock cycle, the way it can for multiply, let alone do multiple per clock cycle like with add/sub. For example, Intel Skylake has throughput of one idiv r32 per 6 cycles (https://uops.info/ / https://agner.org/optimize/). It's microcoded as 10 uops, but in 6 cycles the 4-wide front-end can issue 24 uops, so there's plenty of room for loop overhead of a test/jcc and inc / cmp/jcc in that simple loop that's just running idiv repeatedly to do signed integer division.)

This may be why your unequal work distribution identified by @Peilonrayz doesn't hurt as much as it could: once the "easy" thread finishes, whichever thread has a physical core all to itself can make faster progress.

This dependency on thread-scheduling may also be why your best-3 average (3.33x) is significantly better than your mean of 2.86x speedup. The best case only happens when the last (most work-intensive) thread shares a core with the first thread (finishes first). I'm guessing you do have a CPU with SMT. (But I'm not at all confident in this prediction.)

BTW, uneven chunk costs are what OpenMP schedule(dynamic) is for, when using C with OpenMP to auto-parallelize a loop. (https://stackoverflow.com/questions/10850155/whats-the-difference-between-static-and-dynamic-schedule-in-openmp). Break the problem up into more smaller chunks, and have a worker thread start a new one every time you finish one.

@Jivan Pal's suggestion to interleave the numbers you check across threads is an alternative way to make the work uniform, as long as you skip even numbers! Every even number above 2 is composite. (OTOH, if you insist on brute-force checking every one, then ideally the OS would schedule threads checking even numbers to the same cores as the real work so they don't compete as much. Or at least those 4 threads would exit much sooner than the rest, leaving 4 odd-number checking threads.)

BTW, this is fine as an experiment / benchmark for multi-threading, but if you actually wanted to count primes fast, you'd certainly want to improve the algorithm each thread is running. There's a trivial factor of 2 here, by only checking the odd numbers as divisors: for(int i=3 ; i < rootn ; i+=2). And by only considering odd numbers as candidate primes in the outer loop. On the plus side, you do only go up to sqrt(n) which means you're doing vastly less work than even more naive i < n brute-force loops.

Or the much better option, use a , and yes CodeReview.SE has a whole tag for it. (But that's harder to parallelize.) Of course it doesn't matter how inefficient your algorithm is when you're just comparing relative speedups.

In general, throwing more threads / cores at a problem is sometimes easy to do (for embarrassingly parallel problems like brute-force checking), but making each thread more efficient is also very valuable. (And costs less total CPU cycles / energy.) For this problem specifically, sieving on one core should be faster than brute-forcing on all cores, even with a big Xeon or maybe even a cluster.

Or as @qwr comments, even more specific math for primes (https://en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm) allows counting them without finding them. Prime numbers are kind of a special case, compared to other real-world problems, given how much clever math has been done on them - the algorithmic optimizations are even huger than what's often possible.

For primes, the amount of work scales with $$\O(\sqrt{n})\$$, the number of trial divisions. OTOH primes are less frequent the larger the numbers are: most numbers have some smaller factors, so you aren't frequently going all the way to $$\\sqrt{n}\$$ iterations.

• you should try codegolf.stackexchange.com/questions/74269/… but for even bigger inputs. the best published algorithms are almost O(n^2/3) time – qwr May 9 at 19:44
• @PeterCordes Thank you for this extremely informative answer, I learned a lot from this. And I'm amazed that you managed to actually infer that I have a SMT processor from just looking at the mean and best 3 speed-ups! I do have an Intel SMT processor. I'm on a laptop which has an Intel i5-1135G7. It has 4 physical cores and 8 logical cores. Moreover, as you mentioned, the OS does indeed recognize it as having 8 cores. (I've checked by printing the result of Runtime.getRuntime().availableProcessors() to stdout, and by running cat /proc/cpuinfo on the terminal in Ubuntu.) – wermos May 11 at 6:31
• @wermos: Part of my guesswork was that 4c8t CPUs are very common, but 8c8t CPUs aren't, and I knew you had 8 threads. (Even rarer than I thought; maybe just Intel Coffee Lake-refresh i7 desktop - I think all AMDs ship with SMT enabled. Or I guess ancient dual-socket Core2Quad Xeon systems.) But even without that, the fact that the best was significantly better than the average hinted that there was some inter-thread competition, and SMT is the only plausible mechanism if your system was mostly idle for benchmarking. Not mem or cache. – Peter Cordes May 11 at 17:03
• So the choice between 4c8t or 8c8t was fairly solid once I consider the rarity of 8c8t; most people don't disable HT / SMT in the BIOS. But if I'd had to consider geometries like 5c8t (e.g. a VM exposing 8 cores across 5 physical cores), I wouldn't have been able to estimate the right number of physical cores, just that it was less than #threads. – Peter Cordes May 11 at 17:08
• @wermos: Fun fact: Your Ice Lake / Tiger Lake introduced a new improved integer divide unit that no longer competes with FP division, but it's still 6c throughput. Fewer uops for the front-end (just 4, no longer indirecting to the microcode sequencer). uops.info/…. So naive brute-force trial division still sucks just as much :/ – Peter Cordes May 11 at 17:12