I just came up with this:
// calculates 2^n mod p
uint64_t twopow(uint64_t n, uint64_t p) {
const constexpr uint64_t bitmask {18446744073709551615ull}; // 2^64-1
uint64_t res = 1 << (n % 64);
res *= modpow(1 + (bitmask ^ p), n >> 6, p);
return res % p;
}
What it does is invert the bits of p by XOR-ing them against the integer that's all ones to get 2^64-p-1 which is equivalent to 2^64-1 mod p then add one and you get 2^64 mod p. Then writing n=64*a+b and 2^n=(2^64)^a*2^b.
In general I have p, n between 2^32 and 2^64. In this interval you can't multiply the integers directly or they overflow, so you have to use a workaround such as Montgomery modular arithmetic or Russian peasant multiplication which hurts the performance. The trick with the bitmask reduces n by 6 digits but in general it can still be as large as 2^58 and since the base of the exponential isn't 2 anymore we lose that (possible) advantage. Is there a trick for calculating 2^n mod p which is as fast as repeated squaring would be if p were less than 2^32?
