I am trying to implement exponentiation by squaring of a power of two modulo 2^n − 1. I do not store the power of two, just its exponent (e.g, 4 = 2^2, the exponent is 2).
For example, I want co compute (2^2)^30 (mod 2^5-1)
. The correct answer is 1. However, my method below gives wrong result which is 16. The GMP code:
/* input:
mp_bitcnt_t pos ... the exponent of a base, base = 1 << pos, pos is initially 2 (thus 1 << pos is 4)
const mpz_t exp ... the exponent for squaring (30)
mp_bitcnt_t m_exp ... Mersenne number exponent (5) (thus the modulus is 2^5-1)
*/
while (mpz_cmp_ui(exp, 0UL) > 0) { /* exp > 0 */
if (mpz_odd_p(exp)) { /* exp is odd */
pos += 1;
}
pos = 2 * pos;
mpz_fdiv_q_2exp(exp, exp, 1); /* exp >>= 1 */
pos = pos % m_exp;
}
/* the result is 1 << pos */
Where is the error in my code?
I can compute this exponentiation using the mpz_powm
method and the answer is correct. However, this is EXTREMELY slow. I want to compute this manually.
It turns out that I cannot to modify pos
in an odd step, and I have to introduce a separate variable result
. The correct code is:
mp_bitcnt_t result = 0;
while (mpz_cmp_ui(exp, 0UL) > 0) {
if (mpz_odd_p(exp)) {
result += pos;
result %= m_exp;
}
pos += pos;
pos %= m_exp;
mpz_fdiv_q_2exp(exp, exp, 1);
}
/* the result is 1 << result */