排查基于平方求幂的2的幂模2^n−1计算GMP代码错误

问题背景 #

尝试用**平方求幂（exponentiation by squaring）**实现2^k mod (2^n−1)的快速计算，仅存储2的幂的指数（比如4=2^2就存指数2）。但现有代码计算(2^2)^30 mod (2^5−1)时，正确结果应为1，代码却输出16，需要排查错误并优化（替代速度较慢的mpz_powm）。

错误代码示例 #

/* input:
    mp_bitcnt_t pos ... 底数的指数，底数 = 1 << pos，初始值为2（即底数是4=2^2）
    const mpz_t exp ... 平方操作的指数（30）
    mp_bitcnt_t m_exp ... 梅森数的指数（5）（即模数是2^5-1=31）
*/

while (mpz_cmp_ui(exp, 0UL) > 0) { /* exp > 0 */
    if (mpz_odd_p(exp)) { /* exp 是奇数 */
        pos += 1;
    }

    pos = 2 * pos;

    mpz_fdiv_q_2exp(exp, exp, 1); /* exp >>= 1 */

    pos = pos % m_exp;
}

/* 结果为 1 << pos */

核心错误解析 #

你的代码存在三个逻辑错误：

1. 奇数指数的累积操作完全错误 #

当exp为奇数时，需要把当前底数（对应pos指数）乘到结果中，也就是结果指数 += 当前pos（模m_exp），而不是给pos加1。这一步混淆了平方求幂中「结果累积」的核心逻辑，导致奇数位的贡献完全计算错误。

2. 未区分结果指数与当前底数指数 #

你只用了一个pos变量同时承担「结果存储」和「当前底数指数」的角色，导致结果累积和底数平方的操作互相干扰。正确的平方求幂需要分开维护两个变量：一个存储最终结果的指数，另一个存储当前待平方的底数指数。

3. 循环内操作顺序错误 #

你先执行底数平方（pos=2*pos）再取模，虽然取模顺序不影响结果，但结合前两个错误，整个循环的逻辑链完全混乱，最终得到的pos已经和正确结果毫无关联。

修正后的代码 #

/* input:
    mp_bitcnt_t pos ... 初始底数的指数（例如2，对应2^2=4）
    const mpz_t exp ... 幂次（例如30）
    mp_bitcnt_t m_exp ... 梅森数指数（例如5，对应模数2^5-1=31）
*/

mp_bitcnt_t res_pos = 0;       // 存储最终结果的指数
mp_bitcnt_t curr_pos = pos;    // 存储当前待平方的底数指数
mpz_t exp_copy;
mpz_init_set(exp_copy, exp);   // 复制输入的exp，避免修改原变量

while (mpz_cmp_ui(exp_copy, 0UL) > 0) {
    if (mpz_odd_p(exp_copy)) {
        // 奇数位：将当前底数乘到结果中，指数相加后取模
        res_pos = (res_pos + curr_pos) % m_exp;
    }
    // 底数平方：指数翻倍后取模，避免数值过大
    curr_pos = (curr_pos * 2) % m_exp;
    // exp右移一位，处理下一个二进制位
    mpz_fdiv_q_2exp(exp_copy, exp_copy, 1);
}

mpz_clear(exp_copy);

/* 最终结果：
   - res_pos=0时，2^0=1，符合模2^m_exp-1的结果
   - 否则直接返回1 << res_pos
*/

示例验证 #

针对输入pos=2、exp=30、m_exp=5：

• 30的二进制为11110，循环步骤如下：
  1. exp=30（偶）：curr_pos=(2*2)%5=4，exp=15
  2. exp=15（奇）：res_pos=(0+4)%5=4；curr_pos=(4*2)%5=3；exp=7
  3. exp=7（奇）：res_pos=(4+3)%5=2；curr_pos=(3*2)%5=1；exp=3
  4. exp=3（奇）：res_pos=(2+1)%5=3；curr_pos=(1*2)%5=2；exp=1
  5. exp=1（奇）：res_pos=(3+2)%5=0；curr_pos=(2*2)%5=4；exp=0
• 最终res_pos=0，对应2^0=1，和正确结果一致。

额外优化提示 #

利用梅森数2^n-1的特性：2^k mod (2^n-1)等价于2^(k mod n)（当k mod n≠0时），若k mod n=0则结果为1。你的需求本质是计算(pos * exp) mod m_exp，修正后的代码其实就是用平方求幂的方式快速计算大整数乘法模m_exp，避免了直接计算大整数乘法的溢出问题，效率远高于mpz_powm。

内容的提问来源于stack exchange，提问作者DaBler