仅用C基础语法实现Lucas Lehmer测试至M31的高效方法求助

我刚学C语言一天，目前只掌握了基本数据类型（char、int、float）、while循环、运算符和printf输出。想仅用这些知识实现Lucas Lehmer测试，检测所有p≤31的梅森素数候选数（即Mₚ=2ᵖ-1，其中M₃₁是有符号int的最大值）。

最初的实现代码如下：

#include <stdio.h>
#include <stdbool.h>

int main() {
    int p = 3;
    while (p < 32) {
        int e = 0;
        int Mp = 1;
        while (e < p) {
            Mp *= 2;
            ++e;
        }
        Mp -= 1;
        printf("Testing p=%d: %d\n", p, Mp);
        int i = 0;
        int Si = 4;
        while (i < (p - 2))
        {
            Si = (Si * Si - 2) % Mp;
            ++i;
        }
        if (Si == 0) {
            printf("2^%d is prime: %d\n", p, Mp);
        }
        ++p;
        ++p;
    }
    return 0;
}

这段代码在p=13时出现溢出问题，因为Si的平方超出了有符号int的范围。

为避免溢出，我把乘法改成了重复加法，但效率大幅下降，且只能检测到2¹⁹-1为素数：

int main() {
    int p = 3;
    while (p < 32) {
        int e = 0;
        int Mp = 1;
        while (e < p) {
            Mp *= 2;
            ++e;
        }
        Mp -= 1;
        int i = 0;
        int Si = 4;
        while (i < (p - 2))
        {
            int newSi;
            newSi = 0;
            int multiplications = 0;
            while (multiplications < Si) {
                newSi = (newSi + Si) % Mp;
                ++multiplications;
            }
            Si = newSi - 2;
            ++i;
        }
        if (Si == 0) {
            printf("2^%d is prime: %d\n", p, Mp);
        }
        ++p;
        ++p;
    }
    return 0;
}

我不想用unsigned类型，希望仅用有符号int实现，问是否有高效的方法？

可以实现一个二进制分解的模乘方法，替代直接乘法或循环加法，这样既不会溢出，又能把时间复杂度降到O(log Si)，完全符合你已学的语法知识。

核心思路是：计算(a * b) % m时，把b拆成二进制位，每次只处理一位，同时不断对m取模，保证数值始终在有符号int的范围内。比如：

• 初始化结果为0
• 当b大于0时：
  • 如果b是奇数，就把a加到结果里，然后对m取模
  • 把a左移一位（相当于乘2），然后对m取模（避免a溢出）
  • 把b右移一位（相当于除以2）

把这个逻辑用while循环实现，替换掉原来的乘法或循环加法部分，就能高效计算Si*Si % Mp而不溢出。

修改后的完整代码：

#include <stdio.h>

int main() {
    int p = 3;
    while (p < 32) {
        // 计算Mp = 2^p -1
        int e = 0;
        int Mp = 1;
        while (e < p) {
            Mp *= 2;
            ++e;
        }
        Mp -= 1;
        printf("Testing p=%d: %d\n", p, Mp);

        int i = 0;
        int Si = 4;
        while (i < (p - 2)) {
            // 计算 (Si * Si) % Mp，用二进制分解法避免溢出
            int a = Si;
            int b = Si;
            int mod_result = 0;
            // 二进制分解计算a*b mod Mp
            while (b > 0) {
                if (b % 2 == 1) {
                    mod_result = (mod_result + a) % Mp;
                }
                // a乘2后取模，防止溢出
                a = (a * 2) % Mp;
                // b除以2
                b = b / 2;
            }
            // 更新Si为 (Si² -2) mod Mp
            Si = (mod_result - 2) % Mp;
            // 处理负数情况：如果Si是负的，加上Mp转成正数
            if (Si < 0) {
                Si += Mp;
            }
            ++i;
        }

        if (Si == 0) {
            printf("2^%d -1 is prime: %d\n", p, Mp);
        }
        // 只测试奇素数p（梅森素数的p必须是素数，这里直接跳奇数）
        ++p;
        ++p;
    }
    return 0;
}

代码说明：

1. 模乘部分：通过二进制分解b，每次只处理一位，每次操作后都对Mp取模，确保mod_result和a始终小于Mp，而Mp是2^{p-1，当p≤31时Mp≤2}31-1（有符号int的最大值），所以不会溢出。
2. 负数处理：因为mod_result -2可能得到负数，所以加上Mp转成正数，保证后续计算正确。
3. 效率提升：二进制分解的循环次数是log₂(Si)，比原来的循环Si次快得多，能轻松处理到p=31的情况。

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