OEIS link. This is the sequence of primes such that I first thought of this sequence after reading a mathoverflow post about the Diophantine equation The sequence arises when you consider solving the equation for coprime. In this case you might also expect and to be coprime, in which case you would have for integers because a product of coprime integers which equals twice a square must one be a square and the other twice a square. But you would be wrong to think that. However, the primes that can make that conjecture fail are somewhat rare: 683, 599479, 108390409, 149817457, 666591179, 2000634731, 4562284561, 14764460089, 24040333283, 2506025630791, 5988931115977, and none more up to . OEIS generates these nice graphs: So it looks like the sequence is more or less linearly distributed in the logarithmic graph, so actually exponentially distributed.

I, and someone called sheerluck/Andrew, wrote a program to calculate these numbers. It’s somewhat interesting in its own right. It uses the Sieve of Eratosthenes to generate batches of primes and then spawns a thread to check all the primes in a batch. The version on github is configured to use 4 threads, and it took about a month to check all the way to , but you can edit it to use as many as you want. Each thread then checks each prime in its designated batch one by one. It does this by calculating for each prime and then checking if shares a factor with it. This is more efficient than calculating both orders because we can return early when checking the order of 3 nearly always. The function that calculates the order of 2 is as follows:

std::vector<uint64_t> order_two(std::map<uint64_t, uint64_t> factors, uint64_t p) {
namespace view = std::ranges::views;
uint64_t group_order = p - 1;
std::vector<uint64_t> order;
for (const auto& [P, e] : factors)
{
uint64_t exponent = group_order;
for (const auto f: view::iota(0ull, e + 1))
{
if (modpow_two(exponent, p) != 1)
{
uint64_t alpha = P;
order.push_back(alpha);

break;
}
exponent /= P;
}

}
return order;
}

It uses the fact that the subset of generated by 2 is also a subgroup, thus by Lagrange’s theorem its order divides . So we don’t have to check whether all the powers of two are congruent to 1, just the ones which are divisors of .

Also, the reason why I calculate the order of 2 and then check if it has a common factor with the order of 3 and not the other way around is because calculating the order of 2 is slightly faster using this trick. And I profiled several ways of multiplying 64 bit unsigned ints without overflowing: peasant multiplication, using 128 bit ints (__uint128_t from gcc), using some assembly, and using a trick based on how a * b, the c++ modular multiplication is equal to the integer multiplication mod . The one using assembly turned out to be fastest.

I conjecture that the sequence is infinite. If the analogous sequences for all pairs of naturals (not just as we’ve considered so far) are all finite some real coincidences must happen. So I suspect they’re all infinite.