涉及威尔逊定理的素性检验没有按计划进行

Question

我刚刚开始编程，我遇到了一个我似乎无法解决的问题。我已经编写了这个函数，isPrime似乎总是能通过相等性测试。我可以确认factorial函数可以工作，因为我已经单独测试过它了。

威尔逊定理指出，如果(p - 1)！+1是p的倍数，则p是素数。

#include <stdio.h>
#include <math.h>

void isPrime(double p);
double factorial(double n);

int main(void) {
    double userInput;
    while(1) {
        scanf("%lf", &userInput);
        isPrime(userInput);
    }
    return 0;
}

//

double factorial(double n) {
    if(n <= 1)
        return n;
    else
        return n * factorial(n - 1);
}

void isPrime(double p) {
    if(modf(factorial(p - 1) + 1, &p) == 0)
        printf("Prime!\n");
    else
        printf("Not prime!\n");
}

Answer 1

不要使用doubles来保存整数。舍入误差将使相等测试完全是假的。参见"What Every Computer Scientist Should Know about Floating Point“。

Answer 2

如果n大于23左右，double就不能准确地存储阶乘。从那时起，你的结果将完全是假的。如果要使用威尔逊定理进行素性测试，请使用模算术，

uint64_t modFactorial(uint64_t n, uint64_t m) {
    uint64_t f = 1;
    for(; n > 1; --n) {
        f = (n*f) % m;
    }
    return f;
}

int isPrime(uint64_t p) {
    if (p < 2) return 0;
    return modFactorial(p-1,p) + 1 == p;
}

它可以正确处理小于2^32的输入。威尔逊定理在实践中并不是检验素性的好方法，它的价值在于它在数论中的应用。

如果您需要处理大整数，请使用GMP (它附带-概率-质数测试)。

Answer 3

modf没有像您期望的那样执行模数运算。引述如下：

 The modf() break value into integral and fractional parts, each of which
 has the same sign as the argument.  They return the fractional part, and
 store the integral part (as a floating-point number) in the object
 pointed to by iptr

你想要fmod和朋友，而不是modf。