## Other

drummers-lowrise

Message boards : Number crunching : hey! a more efficient sieve for b*2^k+1

 Subscribe SortOldest firstNewest firstHighest rated posts first
Author Message
composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155166 - Posted: 22 Apr 2022 | 11:41:38 UTC

I've developed a two-part sieve for b*2^k+1 that uses the order of the multiplicative group (mod p) for test primes p that sieves k-values rather than sieving b*2^k+1. It's exponentially faster to sieve k-values so there will be significant speedup for large k. However the trade-off is that determining the group order for a given p is not a trivial amount of work. Fortunately the group order of each p can be computed independently in parallel and fed at any time to a single k-sieve, which runs very fast. I thought about keeping this a trade secret and sieving the heck out of SoB and other conjectures to get a significant LLR testing advantage, but that's not really in keeping with the spirit of community scientific research, and I'm not about to rent a gazillion DO droplets to find prime divisors and compute group orders. So here's how it works.

Program A receives "divisor" primes p (aka potential factors, by whatever method) and outputs the group order and k-value whose b*2^k+1 == 0 (mod p) "the offset", or just the order if p is not a divisor for any k, or "no order" if p divides b. Only a fraction of these will be divisors for a given b. The real beauty of this is that the group order (mod p) is independent of b, so the output of program A can be reused by another instance of program A with a different b value, and this short-circuits the looping because it doesn't need to find the group order. In fact because it's independent of b, the residue logic is identical to finding the group order for 2^k+1. The core logic without too much detail is r0 = b + 1; ... some initial processing for small p and set k = 1 or 2; ... r = r0; h=(p-1)/2; while (k <= h && r != 0 && r != r0) { r = 2 * r - 1; if (r >= p) r -= p; k++;} THAT'S ALL FOLKS. There's no testing for r < 0 since it happens only once in the loop when r == -1 immediately following r == 0, so that test is brought outside the loop for a 10% speedup (which I saw experimentally), and a second loop resumes the search for the order (if not already supplied) sans testing for r == 0 (just until r == r0 or k == h). If the group order hasn't been found found by the time (k == h) then it's automatically (p-1) and the second loop can exit as soon as it finds a 0 residue. When (r == r0) then the group order is (k - 1) (because of k++ loop mechanics), and if the 0 residue hasn't been found then p is not a divisor for any k-value. All the group orders should be included in the output whether p is divisor or not so that it can be reused for another b.

So the group order work is done once for all b values, whether Proth or PSP or SoB or 321. The group order for p is the step size used by the k-sieve with initial k offset determined by program A (only the offset differs for each b value for a give p).

Program A is the one that takes a "long time" to run, finding the group order and offset for around 1,600 candidate primes per minute (single thread) at p ~ 31M. I realize this is very much smaller than the ~ 1P leading edge of sr5sieve, and also much slower than direct sieves on the values b*2^k+1, but the domain being sieved is way smaller (k rather than 2^k) so there's a benefit to be expected for larger k. And also not to berate the benefit of computing the group order of each p just once which is then reused for all other b-values.

Program A processed p up to 51M (~ 3.1 million primes) in 880 minutes producing 1.7 million actual divisors for 321. The output fed to program B sieves the k-values in about 3.4 seconds for k up to 50M. Using divisors up to 51M isn't a very deep sieve, leaving 3.9 million b*2^k+1 values left not factored for k up to 50M.

For comparison, this machine produces primes in 8 GB memory with the Sieve of Eratosthenes up to (2^37) in 977.5 seconds, getting significant speed by storing only odd numbers in the bit map. It might have gone faster but all cores were already churning out group orders and offsets.

I haven't run sr5sieve, so I'd like some feedback on runtimes from experienced users for comparison.

The values of p I processed so far are quite small - they fit in a 64-bit long integer - so I would convert these to using GMP math functions and it would slow down a lot. It's simple to convert since programs A and B total to less than 500 lines of code, compared to the 15k+ lines of code in sr5sieve.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155168 - Posted: 22 Apr 2022 | 18:25:57 UTC

How many candidates are left in PG's sieved files for 321 from the current leading edge of n=19,278,606 up to 50M?

My k-sieve is now at p=67,342,087 and reports it has 2,393,276 candidates left up to 50M.
It's currently eliminating around 50 candidates per minute between that leading edge and 50M.

Edit: It reports 6,287,537 candidates from 19,278,606 up to 100M.
I just recompile program B with a different sieve size and it gives me the result for that in 5 seconds.
It's eliminating around 77 candidates per minute between 50M and 100M.

stream
Volunteer moderator
Volunteer developer
Volunteer tester

Joined: 1 Mar 14
Posts: 1051
ID: 301928
Credit: 563,881,725
RAC: 1,288

Message 155170 - Posted: 22 Apr 2022 | 20:42:02 UTC - in response to Message 155168.

How many candidates are left in PG's sieved files for 321 from the current leading edge of n=19,278,606 up to 50M?

Proth part only (+1) - 1,105,981.

321 was sieved to 95P.

My k-sieve is now at p=67,342,087 and reports it has 2,393,276 candidates left up to 50M.

Is 67,342,087 a maximum prime factor you've tested? If yes, then I have to disappoint you. Right now I'm running sr1sieve on 321 (4 cores) and it advances p by 40,000,000 per second. We're usually throwing away factors below p=1e9 after processing because it's faster to resieve them if we ever will need them again.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155176 - Posted: 23 Apr 2022 | 0:02:10 UTC - in response to Message 155170.

My k-sieve is now at p=67,342,087 and reports it has 2,393,276 candidates left up to 50M.

Is 67,342,087 a maximum prime factor you've tested? If yes, then I have to disappoint you. Right now I'm running sr1sieve on 321 (4 cores) and it advances p by 40,000,000 per second. We're usually throwing away factors below p=1e9 after processing because it's faster to resieve them if we ever will need them again.

When you put it that way, it seems pretty disappointing. But...
(Yes, all primes p up to 67M were tested for the given b, separating them into
a set of proven factors and a set of proven non-factors for a given b.)

TLDR; Program_B output in real time, with web page updates in 15 to 20 second intervals.
live sieve report for b=21181

However I'm doing something quite different than sr1sieve.
Program_A classifies the p's into factors and non-factors
(for any k's; we just don't yet know which k's the proven factors will divide until trying with the k-sieve).
It distills these down to a couple of million good ones which could be used for sr1sieve (for a given b).
When I say I've checked 67M p, it means I know for certain which of these will never produce results
for a particular b, for any k up to infinity. The proven non-factor p's would fail to produce any sieve hits
even with sr1sieve for a particular b, ever. So if you hand me a large prime, then Program_A tests in a
finite amount of time whether it would be good to use with sr1sieve on a given b.

sr1sieve can only prove that p is a factor when it stumbles upon some value b*2^k+1 which it divides.
sr1sieve can't prove the converse so you must retain (or regenerate) all values of p as possible divisors.

The speed of my sieve goes across all k values, for as many 0 <= k <= N as you can hold in a bitmap.
It's very fast to sieve with a proven factor for all 0 <= k <= billions. But this is not so useful for PrimeGrid.
We are both trying to eliminate candidate k values with factors, so that as few as possible are left over
to test for primality with LLR. Your approach is much faster and targets directly the range which will be
submitted for LLR testing, which makes a lot of sense for PrimeGrid because we want to produce
candidates that have the shortest LLR test times, first. It doesn't matter that Program_B can eliminate
huge values of k because there are so many smaller k that are left to test with LLR.

To summarize, we're kind of talking apples and oranges, even though they can both be eaten.

sr1sieve checks a range of flat space between MIN and MAX which contains some pre-initialized
candidates b*2^k+1 with multiple b values and multiple k values which fall in that range, right?

Program_A determines whether p divides *any* value of b*2^k+1, and Program_B tries to find one.
Many of these proven factors don't have a k in the specified sieve range, but we are absolutely certain
some k's will be found given a large enough sieve range.

Program_A reports primes which are proven factors because 0 <= offset < prime for such primes.
Currently I use -1 offset to report a non-factor, but i'll change that to be able to use unsigned integers.

I could run more simultaneous instances of Program_A to increase the rate of checking p & b combinations.
The speed increase for multiple b values occurs in Program_A comes from having the group order supplied to it already.
One instance has to be a leader on this so it can compute the order for other instances.

I've run it this way.

primes 1000000000 > p1G
cat p1G | ./program_A 21181 | tee out.21181 | ./program_A 22699 | tee out.22699 | ...

I don't like the way this works because tee uses buffered output even though program_A flushes
the output for every p record it writes. So the next instance has to wait until the O/S flushes multiple lines of output.

I like it better with

nohup cat p1G | ./program_A 21181 >out.21181 &
nohup tail -f out.21181 | cat out.21181 - | ./program_B 37596884 50000001 >report.21181 &

and then I can casually browse the current status with tail -n 3 out.21181 and with tail -n 3 report.21181
This is pretty close to what's happening to create the web page.

Program_B produces a line only when fil= changes. This is counted as the number of "effective factors".

stream
Volunteer moderator
Volunteer developer
Volunteer tester

Joined: 1 Mar 14
Posts: 1051
ID: 301928
Credit: 563,881,725
RAC: 1,288

Message 155178 - Posted: 23 Apr 2022 | 13:42:19 UTC - in response to Message 155176.

sr1sieve checks a range of flat space between MIN and MAX which contains some pre-initialized
candidates b*2^k+1 with multiple b values and multiple k values which fall in that range, right?

Assume you need to solve b*2^k+1 = 0 (mod p), where b is constant, k is a set of candidates inside sieving range k_min <= k < K, and p is a prime to test.

Of course we're not testing all existing k lying inside sieving range or stored in our candidate file. Instead, rewrite this equation as:

2^k = -1/b (mod p) [eq. 1]

or more algorithmically efficient form

1/(2^k) = -b/c (mod p)

since c=1 for Proth primes, it simplifies further to

1/(2^k) = -b (mod p) [eq. 2]

Equations [1] and [2] are discrete logarithm problem which is solved using BSGS algorithm.

Applying the algorithm and solving the equation, we'll quickly have all possible k values which must be removed from the sieve.

This how sr1sieve / sr2sieve works. sr1sieve is optimized for single sequence, sr2sieve have some tricks to work faster with many sequences (sieving few different b with same k-range).

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155182 - Posted: 23 Apr 2022 | 17:36:50 UTC - in response to Message 155178.

Applying the algorithm and solving the equation, we'll quickly have all possible k values which must be removed from the sieve.

This how sr1sieve / sr2sieve works. sr1sieve is optimized for single sequence, sr2sieve have some tricks to work faster with many sequences (sieving few different b with same k-range).

OK, I didn't realize sr1sieve isn't testing with all primes.
So we're doing the same thing with k values, except I have O(n) search on each prime and BSGS is O(n^.5).
That's why sr1sieve selects primes much faster than Program_A.
How is sr1sieve determining the group order?

My list of possible factors from Program_A for a given b should be identical to yours,
and the sieve step length and starting point should be identical unless we're doing something different.
Can you check what sr1sieve does with the prime 89667313?
The group order I have for this one is 14944552
and the first k which is sieved out is 4605044
so the next k sieved out is 4605044+14944552=19549596
If we don't match then one of us is missing k values, or sieving out too many k's,
or sieving out the wrong k's entirely.

I think the sieve part (Program_B) would be about the same speed except that
to communicate between the 2 parts.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155185 - Posted: 23 Apr 2022 | 22:52:48 UTC - in response to Message 155182.

Applying the algorithm and solving the equation, we'll quickly have all possible k values which must be removed from the sieve.

This how sr1sieve / sr2sieve works. sr1sieve is optimized for single sequence, sr2sieve have some tricks to work faster with many sequences (sieving few different b with same k-range).

OK, I didn't realize sr1sieve isn't testing with all primes.
So we're doing the same thing with k values, except I have O(n) search on each prime and BSGS is O(n^.5).
That's why sr1sieve selects primes much faster than Program_A.
How is sr1sieve determining the group order?

My list of possible factors from Program_A for a given b should be identical to yours,
and the sieve step length and starting point should be identical unless we're doing something different.
Can you check what sr1sieve does with the prime 89667313?
The group order I have for this one is 14944552
and the first k which is sieved out is 4605044
so the next k sieved out is 4605044+14944552=19549596
If we don't match then one of us is missing k values, or sieving out too many k's,
or sieving out the wrong k's entirely.

I think the sieve part (Program_B) would be about the same speed except that
to communicate between the 2 parts.

I neglected to say that the above check is for b=21181

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155186 - Posted: 24 Apr 2022 | 1:27:03 UTC

Wow, k=37598180 has been very resistant to being knocked out of my sieve as a candidate,
as the lowest candidate above the cutoff ever since prime=167.

Is it still in the official "n" candidates list for SoB 21181?
I suppose so, it must be just the next one on the list, as it is likely to have an astronomically large factor.

JimB
Honorary cruncher

Joined: 4 Aug 11
Posts: 920
ID: 107307
Credit: 989,555,179
RAC: 23,710

Message 155198 - Posted: 24 Apr 2022 | 18:38:33 UTC - in response to Message 155186.

21181*2^37598180+1 is divisible by 3799647223

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155199 - Posted: 24 Apr 2022 | 19:22:10 UTC - in response to Message 155198.

21181*2^37598180+1 is divisible by 3799647223

Ah, only 10 digits, not such a big factor after all.
What method did you use to get this factor?

JimB
Honorary cruncher

Joined: 4 Aug 11
Posts: 920
ID: 107307
Credit: 989,555,179
RAC: 23,710

Message 155200 - Posted: 24 Apr 2022 | 19:36:54 UTC - in response to Message 155199.

It was either sr1sieve or sr2sieve. I can't find the spiral notebooks I used at the time.

stream
Volunteer moderator
Volunteer developer
Volunteer tester

Joined: 1 Mar 14
Posts: 1051
ID: 301928
Credit: 563,881,725
RAC: 1,288

Message 155217 - Posted: 25 Apr 2022 | 17:24:45 UTC - in response to Message 155182.

How is sr1sieve determining the group order?

Sorry, this question is above my math skills. As far as I know output of prime generator is fed directly into BSGS, all magic happens there. srsieve /sr1sieve/sr2sieve source is available, version with latest updates was posted somewhere on mersenneforum.org.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155219 - Posted: 25 Apr 2022 | 18:10:08 UTC - in response to Message 155217.

How is sr1sieve determining the group order?

Sorry, this question is above my math skills. As far as I know output of prime generator is fed directly into BSGS, all magic happens there. srsieve /sr1sieve/sr2sieve source is available, version with latest updates was posted somewhere on mersenneforum.org.

It's OK, there's more than one way to do it. I'll look at the code.

You said p advanced 40M per second in your sr1sieve with 4 cores at 95P.
That's not only due to a faster algorithm and implementation.
It's mainly because the density of primes is less at larger primes,
so there are bigger jumps between successive primes.
I'm currently seeing a difference of around 330 between successive primes at p ~ 119M.

With 19 instances of program_A running in parallel my throughput is about 125 primes per second,
so p is advancing at 125 * 330 ~ 41K per second now.

Ravi Fernando
Volunteer tester
Project scientist

Joined: 21 Mar 19
Posts: 211
ID: 1108183
Credit: 14,457,450
RAC: 3,173

Message 155223 - Posted: 26 Apr 2022 | 4:30:51 UTC - in response to Message 155219.

I assume by "group order" you mean the order of 2 in the multiplicative group mod p. This will always be a divisor of p-1, and usually it will be (p-1)/d where d is small. So the first algorithm that comes to my mind is:

1. Factor p-1.
2. Initialize a "result" variable to equal p-1.
3. Loop through all distinct prime divisors q of p-1:
a) Check whether 2^(result/q) is 1 mod p. If so, set result = result/q and repeat until either it fails or result is no longer divisible by q.

4. Return result.

Disclaimer: I haven't read the source code of sr1sieve or any other sieving program, and I am virtually certain that it's smarter than anything I can come up with in 30 seconds.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155224 - Posted: 26 Apr 2022 | 13:28:26 UTC - in response to Message 155223.

I assume by "group order" you mean the order of 2 in the multiplicative group mod p.

Yes. Thanks for that correction.

What I notice is that when the order of 2 in the multiplicative group formed by 2^k+1 (mod p) is p-1,
then 2 is a generator with the same order in the multiplicative group formed by b*2^k+1 (mod p),
and the cycle of residues b*2^k+1 (mod p) for k = 0 .. p-1 is simply a rotation by "e" positions
of the same cycle of residues of 2^k+1 (mod p). I haven't proven this, it's just a pattern I see.

Since p divides 2^{(p-1)/2}+1 the pattern dictates that p also divides b*2^{e+(p-1)/2}+1
so then it's sufficient to have the order of 2 and the rotation "e" to find the offset = e+(p-1)/2 for the sieve.

As a concrete example with p = 29, the order of 2 in multiplicative group 2^k+1 (mod 29) is 28.
With k = 0, 2^k+1 = 2 and with b = 21181, b*2^k+1 = 2 when k = 3, so e = 3.
The k for which p divides 2^k+1 is (p-1)/2 = 14, so with knowledge of the rotation e
we easily calculate that p divides b*2^k+1 with k = (p-1)/2 + e = 17.

The ideal situation would be to have a formula which relates b, p, and e.

rogue
Volunteer developer

Joined: 8 Sep 07
Posts: 1259
ID: 12001
Credit: 18,565,548
RAC: 0

Message 155225 - Posted: 26 Apr 2022 | 15:04:28 UTC

And if you think sr1sieve is fast, I suggest that you give srsieve2cl a try.

Ravi Fernando
Volunteer tester
Project scientist

Joined: 21 Mar 19
Posts: 211
ID: 1108183
Credit: 14,457,450
RAC: 3,173

Message 155226 - Posted: 26 Apr 2022 | 16:33:11 UTC - in response to Message 155224.

What I notice is that when the order of 2 in the multiplicative group formed by 2^k+1 (mod p) is p-1,
then 2 is a generator with the same order in the multiplicative group formed by b*2^k+1 (mod p),
and the cycle of residues b*2^k+1 (mod p) for k = 0 .. p-1 is simply a rotation by "e" positions
of the same cycle of residues of 2^k+1 (mod p). I haven't proven this, it's just a pattern I see.

How I would phrase this: the subgroup H generated by 2 in (Z/pZ)* has the same size as its coset bH. (That is, powers of 2 mod p cycle with the same period as b times powers of 2.) Moreover, if 2 is a generator of (Z/pZ)* (a primitive root), then the list b*2^k is the same as the list 2^k shifted cyclically by e steps, where e (a discrete logarithm) is chosen so that 2^e = 1/b mod p. This is all true, provided of course that b is not itself a multiple of p.
The ideal situation would be to have a formula which relates b, p, and e.

It sounds like you're looking for an efficient solution to the discrete log problem. If you find one, I'm sure the NSA would love to have a word with you.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155227 - Posted: 26 Apr 2022 | 16:53:02 UTC - in response to Message 155226.

What I notice is that when the order of 2 in the multiplicative group formed by 2^k+1 (mod p) is p-1,
then 2 is a generator with the same order in the multiplicative group formed by b*2^k+1 (mod p),
and the cycle of residues b*2^k+1 (mod p) for k = 0 .. p-1 is simply a rotation by "e" positions
of the same cycle of residues of 2^k+1 (mod p). I haven't proven this, it's just a pattern I see.

How I would phrase this: the subgroup H generated by 2 in (Z/pZ)* has the same size as its coset bH. (That is, powers of 2 mod p cycle with the same period as b times powers of 2.) Moreover, if 2 is a generator of (Z/pZ)* (a primitive root), then the list b*2^k is the same as the list 2^k shifted cyclically by e steps, where e (a discrete logarithm) is chosen so that 2^e = 1/b mod p. This is all true, provided of course that b is not itself a multiple of p.
The ideal situation would be to have a formula which relates b, p, and e.

It sounds like you're looking for an efficient solution to the discrete log problem. If you find one, I'm sure the NSA would love to have a word with you.

I believe that they would. I took a couple of undergraduate courses instructed by H. C. Williams at around the time there was a rumour among the students that one of his publications was censored.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155238 - Posted: 27 Apr 2022 | 19:39:13 UTC

It seems to me that I'm doing it backwards, trying to find a prime b * 2k + 1
by using primes p to sieve out all the k's in the test range (except the first such k for each p)
where p divides b * 2k + 1 so that the leftovers can be tested for primality.
The primality test is needed because a leftover k can still produce a composite
which hasn't been identified as such by a higher p than checked so far,
and more than 50% of primes tested fail as they are not factors with that b.

For numbers b * 2k + 1 the recurrence relation is

S0 = b + 1
Sk+1 = 2 * Sk - 1

And this recurrence is performed (mod p).

S0,p = b + 1 (mod p)
S(k+1),p = 2 * Sk,p - 1 (mod p)

For each p, the inner loop searches in t = 0 .. (order of 2) - 1 for a value St,p = 0
and upon finding one, using (2 * t) as the starting offset for an iteration of the k-sieve
with step length = (order of 2) to knock out k's where b * 2k + 1 share the divisor p with b * 2t + 1
or "do nothing" in the k-sieve if this p isn't a divisor for any t (no residue of 0 is found in the cycle).

That's how I do it now, and it's pretty slow. Here's another idea.

Instead of looping (mod p), I could leave Sk as-is,
and attempt to factor this huge number with feasibly small divisors,
using these factors as primes for determining cycle-lengths
to use in the k-sieve with offset k. The unfactorable part is used as if it
is also prime, with some loss of sieve coverage as a tradeoff to speed.

For example the factors of S0 = 21182 are 2, 7, 17, and 89
Since Sk is always odd for k > 0, we disregard 2 and we have
primes 7, 17, and 89 forming cycles with lengths 3, 8, and 11 (respectively)
all having cycle offset 0 to pick off k's in the k-sieve.

The 3 cycle lengths are coprime so there is a long repeating supercycle in the k-sieve
which knocks out (3 * 8 * 11 = 264) k values each supercycle. A simple implementation
just iterates through the whole sieve range for each small cycle length independently.
A space-time tradeoff is realized by implementing the supercycle hits as a list of differences
between offsets within the supercycle, so that the supercycle is used with one pass of
the sieving range, resulting no repeated CPU cache refills of the sieve range, and
sometimes knocking multiple k's in the same cache line.

Unlike the first method above, we never have to search a prime's cycle for a 0 divisor at
some k-offset, as we are given the offset in each iteration of the k-loop.
On the other had, we have to use small trial divisors to produce factors of Sk.
Every value Sk has a divisor at offset k in the cycles corresponding to the factors
factors of Sk. These are the only factors which have this offset. If Sk is
prime, then there's only one factor and it's not being tested for here for primality.

Small factors often reappear as factors of another Sk.
For example 7 is a factor of S0 = 21182 and also of S3 = 169449 = 3 * 7 * 8069.
The k-sieve disregards factors whose cycle lengths (in this case 2, 3, 8068) are <= k and as a result we
would sieve for offset k=3 only with cycle length 8068.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155242 - Posted: 28 Apr 2022 | 5:57:34 UTC

I wondered if that last idea would even be useful, so I took a random small prime, 163

\$ factor 163 163: 163

computed the k-sieve parameters with program_A
\$ echo 163 | ./cycle_filter 21181 pri=163 len=162 off=122

finding that it produces a starting value of 122 for the k-sieve
then computed S122 = 21181 * 2122 + 1 and factored it
\$ echo "2 122 ^ 21181 * 1 + p" | dc | factor 112617512714881212415902149375191913857025: 5 5 163 21673 353869 807707 198857917 22434688079929

Then I ran the unique factors of S122 through program_A
obtaining the unexpected result that all the factors > 163 of S122
have the same starting offset as 163.
\$ echo "2 122 ^ 21181 * 1 + p" | dc | factor | cut -d : -f 2 | tr ' ' '\n' | uniq | ./cycle_filter 21181 pri=5 len=4 off=2 pri=163 len=162 off=122 pri=21673 len=2709 off=122 pri=353869 len=353868 off=122 pri=807707 len=807706 off=122 pri=198857917 len=198857916 off=122 pri=22434688079929 len=934778669997 off=122

THAT'S A FANTASTIC RELATIONSHIP BETWEEN THE FACTORS.
A few more examples have convinced me that this result is general.

Does anyone have a proof for this?
EDIT: I stated this in the previous post but it didn't sink in so deeply due to the density of my brain.

\$ for ((k=12; k<23; k++)); do s=\$(echo "2 \$k ^ 21181 * 1 + p" | dc) echo -e "\nS[\$k] = \$s ==========\n" echo \$s | factor | cut -d : -f 2 | tr ' ' '\n' | uniq | ./cycle_filter 21181 done S[12] = 86757377 ========== pri=7 len=3 off=0 pri=941 len=940 off=12 pri=13171 len=2634 off=12 S[13] = 173514753 ========== pri=3 len=2 off=1 pri=19279417 len=4819854 off=13 S[14] = 347029505 ========== pri=5 len=4 off=2 pri=6469 len=6468 off=14 pri=10729 len=5364 off=14 S[15] = 694059009 ========== pri=3 len=2 off=1 pri=7 len=3 off=0 pri=53 len=52 off=15 pri=71 len=35 off=15 pri=8783 len=4391 off=15 S[16] = 1388118017 ========== pri=11 len=10 off=6 pri=13 len=12 off=4 pri=17 len=8 off=0 pri=19 len=18 off=16 pri=41 len=20 off=16 pri=733 len=244 off=16 S[17] = 2776236033 ========== pri=3 len=2 off=1 pri=29 len=28 off=17 pri=739 len=246 off=17 pri=1489 len=744 off=17 S[18] = 5552472065 ========== pri=5 len=4 off=2 pri=7 len=3 off=0 pri=158642059 len=158642058 off=18 S[19] = 11104944129 ========== pri=3 len=2 off=1 pri=47 len=23 off=19 pri=8750941 len=583396 off=19 S[20] = 22209888257 ========== pri=83077 len=1932 off=20 pri=267341 len=53468 off=20 S[21] = 44419776513 ========== pri=3 len=2 off=1 pri=7 len=3 off=0 pri=23 len=11 off=10 pri=67 len=66 off=21 pri=1372633 len=228772 off=21 S[22] = 88839553025 ========== pri=5 len=4 off=2 pri=89 len=11 off=0 pri=139 len=138 off=22 pri=287251 len=57450 off=22

Ravi Fernando
Volunteer tester
Project scientist

Joined: 21 Mar 19
Posts: 211
ID: 1108183
Credit: 14,457,450
RAC: 3,173

Message 155243 - Posted: 28 Apr 2022 | 7:04:35 UTC - in response to Message 155242.

Then I ran the unique factors of S122 through program_A
obtaining the unexpected result that all the factors > 163 of S122
have the same starting offset as 163.

That seems very expected to me. The offset of p is by definition the smallest k such that p divides 21181 * 2^k + 1. So you found a bunch of prime factors of 21181 * 2^122 + 1, calculated their offsets, and discovered that... all of them divide 21181 * 2^122 + 1, and most of them (the big ones, at least) don't divide any smaller number of the form 21181 * 2^k + 1. What's surprising about that?

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155245 - Posted: 28 Apr 2022 | 13:25:59 UTC - in response to Message 155243.

The offset of p is by definition the smallest k such that p divides 21181 * 2^k + 1.

That's true but not helpful to explaining why the larger factors don't also divide some value with smaller k.

I understand the reason, now, why the larger factors all have offset == k.

It is because the smaller factors have (order of 2) < k so they cover k with k = n * (order of 2) + offset with n > 0,
while the larger factors have (order of 2) > k, and the only way to equal k when n is 0 is to have their offset == k.

I think composite factors would follow this behaviour as well.
Indeed:
\$ echo 25 | ./cycle_filter 21181 pri=25 len=24 off=2 \$ echo "21673 5 * p" | dc | ./cycle_filter 21181 pri=108365 len=10836 off=122

NB the program just prints "pri=", it doesn't expect the divisors to be prime.

Ravi Fernando
Volunteer tester
Project scientist

Joined: 21 Mar 19
Posts: 211
ID: 1108183
Credit: 14,457,450
RAC: 3,173

Message 155249 - Posted: 28 Apr 2022 | 18:10:16 UTC - in response to Message 155245.

It is because the smaller factors have (order of 2) < k so they cover k with k = n * (order of 2) + offset with n > 0,
while the larger factors have (order of 2) > k, and the only way to equal k when n is 0 is to have their offset == k.

I think composite factors would follow this behaviour as well.
Indeed:
\$ echo 25 | ./cycle_filter 21181 pri=25 len=24 off=2 \$ echo "21673 5 * p" | dc | ./cycle_filter 21181 pri=108365 len=10836 off=122

NB the program just prints "pri=", it doesn't expect the divisors to be prime.

Looks correct to me, except that the order of 2 mod 25 is 20, not 24. Maybe a slight bug in the case of divisors with repeated prime factors?

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155251 - Posted: 28 Apr 2022 | 18:49:48 UTC - in response to Message 155249.

Looks correct to me, except that the order of 2 mod 25 is 20, not 24. Maybe a slight bug in the case of divisors with repeated prime factors?

I know why. I assumed that the input is prime, so then assumed that the order is (p-1)
if it isn't less than or equal to (p-1)/2. Obviously I violated that input assumption.
It's a case of garbage in, garbage out. I'll have to redo the order computation to handle composite input.

JeppeSN

Joined: 5 Apr 14
Posts: 1850
ID: 306875
Credit: 52,448,841
RAC: 36,751

Message 155264 - Posted: 30 Apr 2022 | 10:04:36 UTC - in response to Message 155251.

Of course Euler's totient function Ï†(p) has to be used instead of p-1 if p may be nonprime. /JeppeSN

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155270 - Posted: 1 May 2022 | 3:26:35 UTC

I see some interesting symmetries.

There's the obvious rotational symmetry between 2^k+1 (mod p) and b * 2^k + 1 (mod p).
Call the 2^k+1 (mod p) the principal cycle, and the other one the rotated cycle,
which is the same as the principal cycle rotated by "r" elements.

Looking at p = 29

The principal cycle and rotated cycle repeat with the (order of 2) = 28 = (p-1).

The first difference between successive elements of the principal cycle is a cycle of length (p-1)
that contains a pair of shorter cycles each with (p-1)/2 elements, the second one identical to the first,
but with a change of sign of each element, provided that 0 in the principal cycle is replaced with p for
taking differences.

The first difference of the rotated cycle is the same thing as the first difference of the principal cycle,
but rotated by "r" elements.

The "inverse cycle" is constructed by taking the inverse of each element of the principal cycle,
using (p+1)/2 as the inverse of the 0 element to fix symmetries for first differences.
Like the principal cycle, inverse cycle has length (p-1).

The "rotated inverse" cycle is the inverse of the rotated cycle.
It is the same as the inverse of the principal cycle, rotated by "r" elements.

The first difference of the inverse cycle has total length (p-1) and contains a reflection symmetry
at (p-1)/2, providing that the inverse of the 0 in the principal cycle taken as (p+1)/2.

The first difference of the rotated inverse cycle also has total length (p-1) but
splits into 2 shorter reflection symmetries with one rotation of the total cycle.
The first symmetry contains (2r) elements (reflected about r).
This is seen in the "tid" column, as {7 -4 -5 -5 4 7} in the 6 elements after the rotated -1.
The remaining (p-1-2r) elements are in the second reflection below the first, reflecting around (p-1)/2 - r.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155274 - Posted: 1 May 2022 | 13:35:45 UTC - in response to Message 155270.

This is seen in the "tid" column, as {7 -4 -5 -5 4 7} in the 6 elements after the rotated -1.

The rotation by 1 element is superfluous.
It is eliminated by defining the first difference @ k as data[k+1] - data[k].

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155483 - Posted: 19 May 2022 | 7:13:24 UTC - in response to Message 155170.

How many candidates are left in PG's sieved files for 321 from the current leading edge of n=19,278,606 up to 50M?

Proth part only (+1) - 1,105,981.

321 was sieved to 95P.

My k-sieve is now at p=67,342,087 and reports it has 2,393,276 candidates left up to 50M.

Is 67,342,087 a maximum prime factor you've tested? If yes, then I have to disappoint you. Right now I'm running sr1sieve on 321 (4 cores) and it advances p by 40,000,000 per second. We're usually throwing away factors below p=1e9 after processing because it's faster to resieve them if we ever will need them again.

I decided to focus my energies for now on SOB 21181.
I find that conjectures are more appealing for sieve work.
How many candidates are left up to n=50M in SOB 21181, and how high was it sieved to?

I'm currently at p=1.44 billion and have 61'282 candidates remaining for SOB 21181 for 37'760'709 <= n <= 50'000'000.
At this stage it takes about 12 minutes to eliminate 4 candidates in that range by processing 114k primes.

I think we are comparing apples and oranges when talking about sieving.
The sr1sieve file is multi-gigabytes large, correct?
My sieve file is a bit map of candidate exponents up to 50M, so it's just 6 MB.

stream
Volunteer moderator
Volunteer developer
Volunteer tester

Joined: 1 Mar 14
Posts: 1051
ID: 301928
Credit: 563,881,725
RAC: 1,288

Message 155486 - Posted: 19 May 2022 | 11:18:17 UTC - in response to Message 155483.

How many candidates are left up to n=50M in SOB 21181, and how high was it sieved to?

There are 135375 candidates total between 0 and 50M.
There are 32908 candidates left between 37,7M (current leading edge) and 50M.

SoB sieving was one of PG projects, it was run for many months (may be years) on hundreds of computers. It's sieved up to something insane, our log files have factors around p=10000000P.

The sr1sieve file is multi-gigabytes large, correct?
My sieve file is a bit map of candidate exponents up to 50M, so it's just 6 MB.

There is no need to sieve everything from 0 to 50M using sr1sieve. There is a first-stage program - "srsieve". You need only to specify K and N-range to sieve. It'll take just a few minutes to sieve, for example, up to p=1e9. This will eliminate majority of candidates and resulting output file will be only few megabytes long. Factors found during first stage are usually just thrown away, it's much faster to resieve again if necessary. Candidates survived during first stage should be sieved further using "sr1sieve" (tuned for single sequence) or "sr2sieve" (tuned for multiple sequences).

Please note that format "b*2^k+1" which you've used in first post is a bit unusual. In sieving software and everywhere on PrimeGrid we're using "kbnc" format - "k*b^n+c". In case of Proth primes it can be written as "k*2^n+1". In my previous replies I tried to adapt to your notation, but now I finally switched to standard "k" and "n".

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 155488 - Posted: 19 May 2022 | 16:07:00 UTC

Thanks, I'll switch to the standard notation. My title should say k*2^n+1.

For k=21181, I have about double the number of candidates left in both cases, 0 to 50M and 37'7M to 50M.
Given the exponential decay in number of candidates eliminated per prime p,
I have a very looooong way to go to catch up to PrimeGrid.

Sieving the very first divisor (p==3) eliminates 25M candidates (exactly 1/2) from 0 to 50M.
The next divisor (5) eliminates a further 12'5M candidates (exactly 1/4).
Then (7) eliminates 4'17M more candidates (exactly 1/12), leaving 8'3M candidates from 0 to 50M.

The first odd prime to be ignored is p==31.
It has an order of 2 (a cycle length in the 'n' sieve) equal to 5.
Since 31 divides none of those 21181*2^n+1 for n = (0 .. 4) then 31 never divides 21181*2^n+1 for any n.

Sieving to eliminate candidates is not the slow part.
Program B processes the first 79'2 million primes (up to p=1'44e9) in 50 seconds with one thread.
It gets information from Program A on which of these primes to ignore.
If Program A simply didn't output these primes, Program B would have 23'9 million fewer lines of input.
That should save a few seconds. Technically, if we want to keep the ignored primes, they should be
written to a different file than is used as input for Program B.

Computing the remainders to see whether a prime should be ignored is the job of Program A.
Computing each remainder is very fast.
It's just r[0]=21182 (when p > 21182) or r[0] = remainder(21182/p); then r[n+1]=2*r[n]-1 (mod p);
but how many of these remainders there are to calculate depends on the factors of (p - 1).
Program A gradually slows down with larger p, reaching a steady pace where it is rejecting
almost all primes as having an offset (the least n where p divides 21181*2^n+1) which exceeds
the sieve's upper limit.

There is no test in the loop over n to adjust a negative remainder, since the only case is r[n+1] == -1
which immediately follows r[n] == 0, at which point we've found p to be a divisor of 21181*2^n+1.

The first loop continues until (r[n+1] == r[0]) or (r[n+1] == 0) or n reaches the sieve limit (with tricks to reduce it).
If the loop ends on the first condition, we have found the order of 2 without finding a remainder of 0,
so p is not a divisor and we go onto the next prime. If it ends on the third condition then
we have exceeded the range we are sieving and we go onto the next prime. That n and prime and
remainder can be saved for later reuse if the sieve upper limit is ever extended, and we have the
massive disk space to store these residues.

Otherwise p is a divisor so we set r[n+1] == p-1 then check whether it matches r[0] and if not
continue with another loop over n to find the order of 2, once again limiting the size of n with some tricks.

If the second loop terminates before it finds the order of 2, then at least it has found the
first n where r[n] == 0 which potentially eliminates one candidate in the sieve.

When Program A finds the offset AND the order of 2, then Program B loops through the sieve
with that offset and cycle length to try to eliminating multiple candidates for that p.

I strongly doubt that I've eliminated different candidates than the ones left in PrimeGrid's sieve.

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 156226 - Posted: 8 Jul 2022 | 6:19:30 UTC

I've been fiddling with AVX2 in C using Intel intrinsics.
Not quite there yet, but mostly worked out.
It's probably slower than the serial version because modulo requires
so many computations in the vector version.

Without the nitpicky edge cases, the serial way to see if p divides x = k*2^n+1 is:

x0 = k+1 x = x1 = 2 * x0 - 1 n = 1 result = -1 // assume not a divisor while (x != x0) { n++ x = 2 * x - 1 if (x >= p) x -= p // conditional subtraction performs mod p if (x == 0) { result = n break } }

For AVX2 the SIMD calculation skips ahead by 4 serial calculations.
For AVX-512 the vector length would be 8, and SIMD X = 256 * X - 255
x0 = k+1 x1 = 2 * x0 - 1 X = ( x1, 2*x1-1, 2*(2*x1-1)-1, 2*(2*(2*x1-1)-1)-1 ) // Initialize vector with the first 4 iterations of the serial loop // conditional subtractions for modulo P8 = ( 8*p, 8*p, 8*p, 8*p ) P4 = ( 4*p, 4*p, 4*p, 4*p ) P2 = ( 2*p, 2*p, 2*p, 2*p ) P = ( p, p, p, p ) E = ( x0, x0, x0, x0 ) // loop ending condition can occur in any element result = -1 n = 1 while (true) { // SIMD comparisons are tricky, destination argument is replaced // matching elements set to all one bits, otherwise 0 T = E SIMD compare T == X // updates T break if any(X) != 0 n += 4 SIMD X = 16 * X - 15 // In general the SIMD computation is (2^w)X - (2^w-1) with vector length w T = X SIMD compare T >= P8 // updates T SIMD AND T, P8 // updates T SIMD X - T // updates X T = X SIMD compare T >= P4 SIMD AND T, P4 SIMD X - T T = X SIMD compare T >= P2 SIMD AND T, P2 SIMD X - T T = X SIMD compare T >= P SIMD AND T, P SIMD X - T if any(X) == 0 { result = n + "offset of element which is 0" break; } }

composite
Volunteer tester

Joined: 16 Feb 10
Posts: 1172
ID: 55391
Credit: 1,211,016,878
RAC: 1,196,437

Message 156249 - Posted: 9 Jul 2022 | 9:20:49 UTC

The pseudocode in the previous post wasn't quite right.
The actual code snippets in this post are correct.
This isn't the full code, just a chunk for comparing speed.

As expected, the serial version is faster.
It takes around 24% less time than the AVX2 version.

I've begun to think that Intel's vector processing doesn't add more integer ALU resource to a core,
just more registers, and that vector processing is emulated in microcode with the existing serial ALU.
If Intel really added more ALU resource for SIMD, it would be largely wasted since very little code
uses it effectively.

When I get around to it, I'll try an AVX2 version that works on 4 independent primes in parallel
using the serial algorithm so that the modulo code is kept simple. Each element of the vector
must have independent loading and book-keeping since each prime will cause the loop to exit
while other primes are still in progress.

serial version

#include <stdio.h> static long offset; static long cycle_len; static void serial_version (long k, long prime) { const long r0 = k+1; // n == 0 const long ini = 2*k+1; // n == 1 for (int i = 0; i < 100000; i++) { // repetitions for benchmarking offset = -1L; cycle_len = -1L; long vec = ini; // vector of length 1 is a scalar long n = 2; for (; n <= prime - 1; n++) { vec = vec * 2 - 1; if (vec >= prime) vec -= prime; if (vec == 0) { offset = n; break; } if (vec == r0) { cycle_len = n+1; break; } } } }

AVX2 version
#include <immintrin.h> #include <stdlib.h> #include <stdio.h> static long offset; static long cycle_len; static void avx2_version (long k, long prime) { /* memory for printing a vector */ __m256i *mem = (__m256i *) aligned_alloc(32, 32); // aligned on 32-byte boundary long *ptr = (long *) mem; #define STORE(x) _mm256_store_si256(mem, x) const __m256i c0 = _mm256_set_epi64x(0L, 0L, 0L, 0L); const __m256i c15 = _mm256_set_epi64x(15L, 15L, 15L, 15L); const __m256i pt1 = _mm256_set_epi64x(prime, prime, prime, prime); const __m256i pt2 = _mm256_set_epi64x(prime*2, prime*2, prime*2, prime*2); const __m256i pt4 = _mm256_set_epi64x(prime*4, prime*4, prime*4, prime*4); const __m256i pt8 = _mm256_set_epi64x(prime*8, prime*8, prime*8, prime*8); const __m256i pt1m1 = _mm256_set_epi64x(prime-1, prime-1, prime-1, prime-1); const __m256i pt2m1 = _mm256_set_epi64x(prime*2-1, prime*2-1, prime*2-1, prime*2-1); const __m256i pt4m1 = _mm256_set_epi64x(prime*4-1, prime*4-1, prime*4-1, prime*4-1); const __m256i pt8m1 = _mm256_set_epi64x(prime*8-1, prime*8-1, prime*8-1, prime*8-1); /* initial remainder */ const __m256i r0 = _mm256_set_epi64x(k+1, k+1, k+1, k+1); // n == 0 /* next remainder in 4 parallel streams */ const __m256i ini = _mm256_set_epi64x(2*k+1, 4*k+1, 8*k+1, 16*k+1); // n = (1, 2, 3, 4) __m256i tst; for (int i = 0; i < 100000; i++) { // repetitions for benchmarking offset = -1L; cycle_len = -1L; __m256i vec = ini; // STORE(vec); // printf("vec: %ld %ld %ld %ld\n", ptr[3], ptr[2], ptr[1], ptr[0]); long n = 5; for (; n <= prime - 1; n += 4) { // vec = vec * 16 - 15; vec = _mm256_sub_epi64(_mm256_slli_epi64(vec, 4), c15); // vec = mod(vec, prime); vec = _mm256_sub_epi64(vec, _mm256_and_si256(_mm256_cmpgt_epi64(vec, pt8m1), pt8)); vec = _mm256_sub_epi64(vec, _mm256_and_si256(_mm256_cmpgt_epi64(vec, pt4m1), pt4)); vec = _mm256_sub_epi64(vec, _mm256_and_si256(_mm256_cmpgt_epi64(vec, pt2m1), pt2)); vec = _mm256_sub_epi64(vec, _mm256_and_si256(_mm256_cmpgt_epi64(vec, pt1m1), pt1)); // STORE(vec); // printf("vec: %ld %ld %ld %ld\n", ptr[3], ptr[2], ptr[1], ptr[0]); // if (any(vec) == 0) tst = _mm256_cmpeq_epi64(vec, c0); if (! _mm256_testz_si256(tst, tst)) { STORE(tst); if (ptr[3] != 0) offset = n; if (ptr[2] != 0) offset = n+1; if (ptr[1] != 0) offset = n+2; if (ptr[0] != 0) offset = n+3; break; } // if (any(vec) == k+1) tst = _mm256_cmpeq_epi64(vec, r0); if (! _mm256_testz_si256(tst, tst)) { STORE(tst); if (ptr[3] != 0) cycle_len = n; if (ptr[2] != 0) cycle_len = n+1; if (ptr[1] != 0) cycle_len = n+2; if (ptr[0] != 0) cycle_len = n+3; break; } } } }

Message boards : Number crunching : hey! a more efficient sieve for b*2^k+1