While implementing Paul Underwood's algorithm[1] for a small machine with only 16 bit native wordsize I found that $y$ ($y$ in OEIS, $a$ in the paper) cannot exceed 177 without using big integers and I need to find out, how many primes would slip through with that limit.
As computing time is cheap these days I took the liberty to just use brute-force to find the $y$ for the primes up to $2^{39}$. I also kept trace of the amount of primes per $y$ (starting at zero).
[10575966557, 5287990242, 0, 2643997054, 0, 1321982935, 661002427, 0, 0,
330494926, 0, 165259666, 0, 0, 0, 82615577, 0, 41313770, 0, 0, 0,
20655767, 0, 0, 0, 0, 0, 10331891, 0, 5163031, 0, 0, 0, 0, 0, 2581046, 0,
0, 0, 1289199, 0, 642747, 0, 0, 0, 318400, 0, 0, 0, 0, 0, 156966, 0, 0,
0, 0, 0, 76550, 0, 37065, 0, 0, 0, 0, 0, 17491, 0, 0, 0, 8071, 0, 3687,
0, 0, 0, 0, 0, 1655, 0, 0, 0, 709, 0, 0, 0, 0, 0, 314, 0, 0, 0, 0, 0, 0,
0, 131, 0, 0, 0, 42, 0, 24, 0, 0, 0, 13, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0]
(The overall sum is too large by 8 but that was my error, I ran it in chunks and started the new chunks at $p$ instead of $p+1$).
Computing the distribution of the numbers above gave a large hint that it might be $1/2^n$ with $n$ being the count of $y$. Example:
$y = 0$ in $1/2^1$
$y = 1$ in $1/2^2$
$y = 3$ in $1/2^3$
$y = 5$ in $1/2^4$
...
$y = 21$ in $1/2^{10}$ of all cases
and so on.
Is that actually the case and if so, why?
Also based on that but maybe different enough for a new question: is there a way to compute the probability of $y > 177$, that is, is there a simple way to count all $y$ up to and including 177 for primes?
[1] Underwood, Paul. "Quadratic Frobenius probable prime tests costing two selfridges." arXiv preprint arXiv:1706.01265 (2017).
I would say the probability of y (or a) exceeding 177 is about
1/2^(21/32*177) = 1/2^(112.875) = 1/10^34
based on the first 21 values for minimal y (or a) given in the paper, although there might be some log(log) adjustments.
The maximum for odd n<2^50 was 81, and 50*32/21. = 76.19.