使用 Pari 的卡迈克尔数

Question

尝试编写 Pari 代码来解决上述问题。

Answer 1

我没有使用 Pari 的经验，但这里有一些有用的建议：

n 是 Carmichael 当且仅当它是复合的，并且对于所有与n 互质的a 和1 < a < n，同余a^(n-1) = 1 (mod n) 成立。要直接使用此定义，您需要：

1) 一种测试a 和n 是否相对质数的有效方法

2) 一种计算a^(n-1) (mod n)的有效方法

对于第一个 - 使用 Euclidean algorithm 作为最大公约数。它在循环中最有效地计算，但也可以通过简单的递归gcd(a,b) = gcd(b,a%b) 和基础gcd(a,0) = a 来定义。在 C 中，这只是：

unsigned int gcd(unsigned int a, unsigned int b){
    return b == 0? a : gcd(b, a%b);
}

对于第二点——在计算a^k (mod n) 时，您可以做的几乎最糟糕的事情是首先通过重复乘法计算a^k，然后通过n 修改结果。相反——使用exponentiation by squaring，在中间阶段取余数（mod n）。它是一种基于观察的分而治之的算法，例如a^10 = (a^5)^2 和 a^11 = (a^5)^2 * a。一个简单的 C 实现是：

unsigned int modexp(unsigned int a, unsigned int p, unsigned int n){
    unsigned long long b;
    switch(p){
        case 0:
            return 1;
        case 1:
            return a%n;
        default:
            b = modexp(a,p/2,n);
            b = (b*b) % n;
            if(p%2 == 1) b = (b*a) % n;
            return b;
        }
} 

注意unsigned long long在计算b*b时使用，防止溢出。

要测试n 是否为 Carmichael，您不妨先测试n 是否为偶数，在这种情况下返回0。否则，在2 到n-1 的范围内逐个遍历数字a。首先检查gcd(a,n) == 1 是否是复合的，如果n 是复合的，那么在达到n 和gcd(a,n) > 1 的平方根之前，您必须至少有一个a。保留一个布尔标志，用于跟踪是否遇到过这样的a，如果超过平方根而没有找到这样的a，则返回0。对于a 和gcd(a,n) == 1，计算模幂a^(n-1) (mod n)。如果这与 1 不同，则返回 0。如果您的循环完成检查所有n 下方a 而不返回0，那么数字是Carmichael，所以返回1。一个实现是：

int is_carmichael(unsigned int n){
    int a,s;
    int factor_found = 0;
    if (n%2 == 0) return 0;
    //else:
    s = sqrt(n);
    a = 2;
    while(a < n){
        if(a > s && !factor_found){
            return 0;
        }
        if(gcd(a,n) > 1){
            factor_found = 1;
        }
        else{
            if(modexp(a,n-1,n) != 1){
                return 0;
            }
        }
        a++;
    }
    return 1; //anything that survives to here is a carmichael
}

一个简单的驱动程序：

int main(void){
    unsigned int n;
    for(n = 2; n < 100000; n ++){
    if(is_carmichael(n)) printf("%u\n",n);
    }
    return 0; 
}    

输出：

C:\Programs>gcc carmichael.c

C:\Programs>a
561
1105
1729
2465
2821
6601
8911
10585
15841
29341
41041
46657
52633
62745
63973
75361

这仅需大约 2 秒即可运行，并且与 this 列表的初始部分匹配。

这可能是一种比较实用的方法来检查高达一百万左右的数字是否是卡迈克尔数。对于较大的数字，您可能应该让自己获得一个好的因式分解算法并使用 Korseldt 的标准，如在 Carmichael 数字上的 Wikipedia entry 中所述。