查找从 2 到 1000 的所有素数的算法不起作用答案

【问题标题】：Algorithm to find all primes from 2 to 1000 not working查找从 2 到 1000 的所有素数的算法不起作用
【发布时间】：2016-02-02 16:41:18
【问题描述】：

这是一段代码，使用以下语句计算从 2 到 1000 的所有素数，即数字 n 是素数 iff：

在第一个版本中，我认为我正确地实现了算法：

public class Giuga {

    public static void main(String[] args){
        int n = 2;

        while(n<=1000){
            int k = 1;
            long sum = 0;
            while(k<=n-1){
                sum = sum+(long)Math.pow((double)k,(double)n-1);
                k++;
            }
            if(sum%n==n-1){
                System.out.println(n + " is a prime.");
            }
            n++;
        }
    }
}

但是，由于变量sum 快速增长，会发生溢出，在素数 17 之后将不再有输出。

为了防止我不得不使用这个：

嗯，我做到了，这是我的 2. 版本：

public class Giuga {

    public static void main(String[] args){
        int n = 2;

        while(n<=1000){
            int k = 1;
            long sum = 0;
            while(k<=n-1){
                sum = sum+((long)Math.pow((double)k%n,(double)n-1))%n; //Here are the changes
                k++;
            }
            if(sum%n==n-1){
                System.out.println(n + " is a prime.");
            }
            n++;
        }
    }
}

我认为我做对了，但现在输出在素数 13 之后停止。

我试图找出我的错误已经有一段时间了。我究竟做错了什么？从 2 到 1000 必须有 168 个素数。

【问题讨论】：

嗯...这个方程看起来就像适用于所有可能基的费马小定理。但如果是这样，那么它只是一个必要条件，而不是一个足够找到所有素数的条件。 Charmicael 数也应满足该等式。

标签： java algorithm primes

【解决方案1】：

正如已经指出的那样，doubles 的精度只有大约 16 位，不够精确，无法在足够高的数字上保持正确的余数。

您可以切换到longs 并执行您自己的模幂运算。

int k = 1;
long sum = 0;
while(k<=n-1){
    long pow = 1;
    for (int i = 0; i < n - 1; i++)
        pow = (pow * k) % n;
    sum = (sum + pow)%n;
    k++;
}

可以通过将这种简单的模幂运算更改为通过重复平方使用模幂运算来改进该算法，它不是最有效的素数查找算法，但它现在是正确的。

2 is a prime.
3 is a prime.
5 is a prime.
7 is a prime.
11 is a prime.
13 is a prime.
17 is a prime.
19 is a prime.
23 is a prime.
29 is a prime.
31 is a prime.

（截图）

977 is a prime.
983 is a prime.
991 is a prime.
997 is a prime.

要通过重复平方使其模幂，替换

long pow = 1;
for (int i = 0; i < n - 1; i++)
    pow = (pow * k) % n;

与

long pow = 1;
long square = k;
int exp = n - 1;
while (exp > 0)
{
    if ((exp & 1) == 1)
    {
        pow = (pow * square) % n;
    }
    square = (square * square) % n;
    exp >>= 1;
}

它连续测试指数的每一位，如果设置了，则将当前平方乘以pow。

【讨论】：

【解决方案2】：

Java doubles（即pow 的输出）不能精确地表示足够大的数字以产生正确的余数。你应该切换到modular exponentiation。

【讨论】：

但是除了 Java 中的 BigInteger 之外，还有其他使用模幂的方法吗？
@de_dust 自己实现？平方的方法很快，大概10行代码。
@de_dust 大约在“从右到左二进制方法”下的 wiki 文章中包含该方法的伪代码。把它翻译成Java，你就会得到答案！刚刚自己测试过，它就像一个魅力。

【解决方案3】：

您始终可以使用 BigInteger 类进行大量计算

private static boolean isPrime(int n) {
    BigInteger N = new BigInteger(String.valueOf(n));
    BigInteger N_MINUS_1 = N.subtract(BigInteger.ONE); 

    BigInteger sum = BigInteger.ZERO;
    for (int k = 1;  k < n; k++)
        sum = sum.add(new BigInteger(String.valueOf(k)).modPow(N_MINUS_1,N)).mod(N);
    return sum.equals(N_MINUS_1);
}

有趣的是这是Fermat's little theorem的变体，对于和Σ中的每个k

k^(n-1)%n 应该是 1，否则这个数不是素数！因此，如果我们找到 k^(n-1)%n != 1 的 k，我们就可以停止计算。上述算法可改写为：

private static boolean isPrimeFermat(int n) {
    BigInteger N = new BigInteger(String.valueOf(n));
    BigInteger N_MINUS_1 = N.subtract(BigInteger.ONE); 

    for (int k = 1;  k < n; k++){
        if (new BigInteger(String.valueOf(k)).modPow(N_MINUS_1, N).equals(BigInteger.ONE) == false)
            return false;
    }       
    return true;
}

瞧！

【讨论】：

它有效，但我不太明白，为什么每个 summand 都必须为 1（mod n）。你能解释一下吗？我明白，如果每个加数都是 1（mod n），它是一个素数，但为什么这是唯一的情况？
旁注：这当然比我的解决方案要快得多（因为它是由一群聪明的软件工程师而不是我写的^^）。
这实际上是费马小定理的总和变分
感谢您将其添加到答案中。所以基本上你看到加数对应于一个可以检查费马小定理的表达式。我虽然您从 OP 初始定理 (sum mod n==n-1 <=> n is prime) 推导出它 (a^k mod n == 1 forall a <=> n is prime)，但不知道您是如何做到的。 :)

【解决方案4】：

你的力量也太大了（猜猜 999^999 是什么）。因此，您还必须使用逐步的(ab) mod n = ((a mod n)(b mod n)) mod n（模幂运算）来计算功率：

public class Giuga {

    public static void main(String[] args) {
        int count = 0;
        for (int i = 2; i < 1000; i++) {
            if (isPrime(i)) {
                count++;
                System.out.printf("%d is prime\n", i);
            }
        }
        System.out.printf("Found %d primes.\n", count);
    }

    // checks if number is a prime
    private static boolean isPrime(int number) {
        int bigSum = 0;
        for (int k = 1; k <= number - 1; k++) {
            bigSum = (bigSum + calcPowMod(k, number - 1, number)) % number;
        }
        return bigSum % number == number - 1;
    }

    // calculates (a^b)%mod, making sure that intermediate results are always < max(a^2,mod)
    private static int calcPowMod(int a, int b, int mod) {
        int pow = 1;
        for (int k = 1; k <= b; k++) {
            pow = (pow * a) % mod;
        }
        return pow;
    }
}

【讨论】：