Build a smaller hash table for ``n_is_prime``

[fredrik-johansson]

Here is the code I used to generate the list of bases for the hash table in n_is_prime: https://gist.github.com/fredrik-johansson/2c568b0b20675a22d77f89b9abdf3215

One can get the necessary list of strong base-2 pseudoprimes up to $2^{64}$ from https://flintlib.org/data/psp2.txt.xz (warning: 186 MB compressed)

Generating the hash table for n = 131072 buckets was reasonably quick, a few hours on my 8-core machine IIRC.

With some patience, one could build a table with smaller n. For example, n = 98304 = 3 * 2^15 looks feasible. Running this for 8 hours on 8 cores completes about 10% of the buckets:

...

max: 399,  average: 324
bucket 5 / 98304 : 21450  (max 21450 for 294 histogram [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1])
bucket 6 / 98304 : 308318  (max 308318 for 295 histogram [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1])
bucket 14 / 98304 : 82434  (max 308318 for 295 histogram [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1])
...
bucket 8052 / 98304 : 12720971  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 61, 122, 206, 454, 760, 1033, 1239, 1290, 1133, 867, 483, 239, 91, 36, 11, 4, 1])
bucket 9522 / 98304 : 370031  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 61, 122, 206, 454, 760, 1034, 1239, 1290, 1133, 867, 483, 239, 91, 36, 11, 4, 1])
bucket 9530 / 98304 : 14197  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 760, 1034, 1239, 1290, 1133, 867, 483, 239, 91, 36, 11, 4, 1])
bucket 9538 / 98304 : 182283  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 761, 1034, 1239, 1290, 1133, 867, 483, 239, 91, 36, 11, 4, 1])
bucket 9546 / 98304 : 593103  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 761, 1034, 1240, 1290, 1133, 867, 483, 239, 91, 36, 11, 4, 1])
bucket 7889 / 98304 : 27939486  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 761, 1034, 1240, 1290, 1133, 867, 483, 240, 91, 36, 11, 4, 1])
bucket 7897 / 98304 : 136819  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 762, 1034, 1240, 1290, 1133, 867, 483, 240, 91, 36, 11, 4, 1])
bucket 7912 / 98304 : 2645123  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 762, 1034, 1240, 1290, 1134, 867, 483, 240, 91, 36, 11, 4, 1])
bucket 7905 / 98304 : 548314  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 762, 1034, 1241, 1290, 1134, 867, 483, 240, 91, 36, 11, 4, 1])
bucket 7913 / 98304 : 507402  (max 736820777 for 372 histogram [0, 0, 0, 0, 0, 0, 0, 1, 2, 4, 5, 11, 18, 42, 62, 122, 206, 454, 762, 1035, 1241, 1290, 1134, 867, 483, 240, 91, 36, 11, 4, 1])

From the histogram one sees that about 95% of the bases fit in 24 bits, so this table could presumably be stored in roughly 98304 * 24 bits ~= 288 KB. If anyone is interested in finishing this computation, it should take a few days on a multicore machine.

One drawback of choosing a non-power-of-two number of buckets is that the hash function is slightly slower to compute, but probably this tradeoff is acceptable.

The table could presumably be optimized a tiny bit by tweaking the hash function to reduce the number of outliers.

Could one reach even smaller n? 2^16 would be a nice target. After 20 core-hours, four buckets complete on my machine:

max: 579,  average: 486
bucket 3 / 65536 : 402547707  (max 402547707 for 478 histogram [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1])
bucket 2 / 65536 : 442735750  (max 442735750 for 435 histogram [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2])
bucket 10 / 65536 : 392401781  (max 442735750 for 435 histogram [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3])
bucket 18 / 65536 : 111066745  (max 442735750 for 435 histogram [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 3])

Ideally most bases will fit in 32 bits in which case this table would fit nicely in 256 KB, though I'm not yet sure if this will actually be the case. It looks like this will take several core-years to complete, which could the target of a distributed computing project.

Question: is there some way to speed up this search substantially? One idea is to run the probable prime test for a SIMD vector of several bases simultaneously. There should also be no problem doing this search on a GPU (you need 64-bit integers, but even slow 64-bit integers would be fine with massive parallelism). Are there some pruning, multipoint evaluation or meet-in-the-middle tricks one could use here?

[math-gout]

I just ran your code with 98304 buckets (took 48 hours on a Xeon E5-2698v4 with 40 threads). I uploaded the results here:
https://gist.github.com/math-gout/3ef9b109e770e5ccc66777edb9fd412d
Note I had to remove the additional information in the parenthesis as the file was too big and I couldn't create a gist. Let me know if you need the full file, I'll try to find another way.

[fredrik-johansson]

Thanks, works perfectly! Updated code in #2507.

I've done some more experiments with n = 2^16 and while generating such a table probably is feasible, it looks unfortunately like the average entry will be a bit larger than 32 bits. The advantage of a fast hash function for a power of two size will then probably be negated by the extra logic to bit-unpack such a table. It might be best to determine an n > 65536 for which most entries fit as snugly as possible in 32 bits.

Another observation: only 1/3 of the bases in the 98304 table are even. It might be possible to speed up the search and gain a fraction of a bit of storage per entry by restricting to odd bases and storing these with an implicit LSB.

[JASory]

You can make this problem easier by eliminating pseudoprimes of certain forms. Checking if an integer is a semiprime of the form (2x+1)(4x+1) can be done considerably faster than a full fermat test. You can also check other similar forms that have many liars like (2x+1)(6x+1). You should be able to produce a 16384-entry table with 16-bit witnesses relatively easily using a total of 32768 bytes. I previously used this approach, but abandoned it in favor of an 262144 element table of 16-bit witnesses calibrated to work by itself against particularly small integers.

I would personally use a reduced operation Lucas test in addition to the base-2 pseudoprime, it's only slightly slower and avoids the extensive tables.

[math-gout]

Hello again!

I did a little more investigating, but haven't implemented anything in FLINT and do not have timings. Here are my results so far:

• Starting with @JASory comment, we can actually generalize to finding liars of the form n = p*q with q = k*(p-1)+1. To find good k values I factored all the liars and checked for this pattern (also allowing for fractions, just in case). These can be caught with something like:

ulong sqrtn = n_sqrt(n);
if (n % ((ulong)(sqrtn/sqrt(k)) + 1) == 0) return 0;

(I double checked that this works by comparing to what I got from the factorization)
(Also, is n_sqrt faster than just the regular sqrt function?)
I ended up choosing k = 2,3,4,5, which reduced the number of liars from 31894014to 24745862 (adding more k values doesn't look promissing as the number of liars we catch that way drops quite quickly). From that I was able to generate a 100kB hash table with 49152 16 bits entries (532 entries do not fit), and a 98kB one with 32768 24 bits (22 entries do not fit).

• Those two tables are essentially the same size, so if we were to keep one of them we need to check wether we should do the 49152 mod, or the 24 bit packing.

• I also tried to force the witnesses to be odd to save the LSB, but I didn't seem to be worth it as the average size increases a bit too much and loosing the bit we just saved.

• Finally, I tried a few different values for the multiplication in the hash function. All values I tried gave similar results, so unless we do a lot of compute and get lucky enough to find a really good one, we might as well remove the multiplication (which is what I did in the two tables from the first point).

Let me know if you are interested in those tables or in me adding them to FLINT in a branch for testing!

[fredrik-johansson]

Interesting. Feel free to implement this in a branch.

I'm concerned that computing a square root and doing several divisions is going to be much slower than the current test, but would be interesting to see just how much.

(Also, is n_sqrt faster than just the regular sqrt function?)

n_sqrt just does sqrt and then corrects the rounding, so n_sqrt is certainly slower.