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.

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.

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.

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.

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?

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).

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
sr1sieve has a significant performance advantage with no I/O overhead
to communicate between the 2 parts.

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.

21181*2^37598180+1 is divisible by 3799647223

21181*2^37598180+1 is divisible by 3799647223

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

How is sr1sieve determining the group order?

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.

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.

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

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

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.

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

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.

It seems to me that I'm doing it backwards, trying to find a prime b * 2^k + 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 * 2^k + 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 * 2^k + 1 the recurrence relation is

S₀ = b + 1
S_k+1 = 2 * S_k - 1

And this recurrence is performed (mod p).

S_0,p = b + 1 (mod p)
S_(k+1),p = 2 * S_k,p - 1 (mod p)

For each p, the inner loop searches in t = 0 .. (order of 2) - 1 for a value S_t,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 * 2^k + 1 share the divisor p with b * 2^t + 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 S_k 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 S₀ = 21182 are 2, 7, 17, and 89
Since S_k 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 S_k.
Every value S_k has a divisor at offset k in the cycles corresponding to the factors
factors of S_k. These are the only factors which have this offset. If S_k 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 S_k.
For example 7 is a factor of S₀ = 21182 and also of S₃ = 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.

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

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

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

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?

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.

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

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.

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

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.

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.

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:

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