I'm the author of that essay.

I mentioned in passing in a couple of locations that the code I describe is only an approximation to the code which I am actually optimizing. The fingerprint objects are bit patterns of (in this case) 2048 bits, and the similarity between two fingerprints is defined as the number of bits in the intersection divided by the number of bits in the union. (This is the "Tanimoto" or "Jaccard" similarity.) Simple math says that if the one fingerprint has B bits, and the threshold is t, then fingerprints to search must have between Bt and B/t bits.

What I do then to speed up the search is to precompute the popcounts for each fingerprint, sort them by popcount, and make an auxiliary index so I know where the start and end of each popcount count region is. As a consequence, this reduces the O(NN) term by a factor of two. With that in place, my sort algorithm is closer to:

for (row=0; row<n; row++) {
query_popcount = popcount(fingerprint[row]);
start_popcount = num_bits*threshold;
end_popcount = num_bits/threshold+1;
for (target_popcount = start_popocount; target_popcount<end_popcount; target_popcount++) {
for (col = start_of_popcount[target_popcount]; col<start_of_popcount[target_popcount+1]; col++) {
num_bits_in_commmon = intersection_popcount(fingerprint[row], fingerprint[col]);
if (num_bits_in_common / double(query_popcount + target_popcount - num_bits_in_common > threshold) {
count++;
}
}
}
}

This is less comprehensible than the version I posted in the essay, but still has the essential problem I wanted to resolve.

With these optimizations in place, it's more complicated to determine the number of calculations which will be done in a row. For one, the popcount distribution is highly data dependent. I don't know that number until I'm in an inner loop, and so it isn't worthwhile to try and attempt to assign a "constant amount of effort per Unit work."

In the usual case of high thresholds, and if the fingerprint popcount distribution is broad, then I am mostly working along the near-diagonal of the matrix, which improves my cache coherency. However, if I have enough fingerprints then this simple linear scan will blow the cache, and the above needs to be replaced by some cache oblivious search pattern. If that's what you meant by "2d shapes" in your comment then perhaps you missed my mention of Morton Z-ordering?

My goal, however, was to explain the thinking needed to minimize lock contention and the difficulties in developing a parallel algorithm, so I used a simplified version of what I actually did to get that point across. Your comments deal with other aspects of the algorithm, and don't contribute to the point.

BTW, what do you mean by "insane number of threads"? The target is a 8 core machine, and my timings numbers show a scaleup of about 6x for that case. I'm used to working on parallel molecular dynamics codes which use 10-100 processors, so this seems rather small.

Also, "he can just keep a separate count per worker and then add it to the total when that worker finishes" isn't correct, or rather, I can't see how to apply it to the symmetric form. The code I came up with has one single-threaded code, of O(N) time, and after the final barricade. I don't see how any other solution has less contention.