OEIS link. This is the sequence of primes p such that

    \[\gcd(\text{ord}_p(2), \text{ord}_p(3)) = 1.\]

I first thought of this sequence after reading a mathoverflow post about the Diophantine equation

    \[(2^x-1)(3^y-1)=2z^2.\]

The sequence arises when you consider solving the equation for x,y coprime. In this case you might also expect 2^x-1 and 3^y-1 to be coprime, in which case you would have

    \[2^x-1=p^2\quad 3^y-1=2q^2\]

for integers p,q 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 10^{13}. 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 10^{13}, 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 \text{ord}_p(2) for each prime and then checking if \text{ord}_p(3) 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 (\mathbb{Z}/p\mathbb{Z})^\times generated by 2 is also a subgroup, thus by Lagrange’s theorem its order divides p-1. So we don’t have to check whether all the powers of two are congruent to 1, just the ones which are divisors of p-1.

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 ab the integer multiplication mod 2^{64}. 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 (\alpha,\beta) (not just (2,3) as we’ve considered so far) are all finite some real coincidences must happen. So I suspect they’re all infinite.