费马-欧拉定理和快速幂算法

费马小定理（Fermat’s little theorem）：当p是质数，且a与p互质时，a^(p-1)≡1(mod p)。这其实是欧拉（Euler）定理的一个特例。所以，费马小定理又被称为费马-欧拉定理。

欧拉定理：若a,n是正整数，且a,n互质，则：a^φ(n)≡1(mod n)，其中φ(n)称为欧拉函数，指的是小于n且与n互质的整数的个数。欧拉定理的证明并不复杂。

设整数x₁、x₂、……、x_i、……、x_φ(n)是所有小于n且与n互质的整数，一共φ(n)个。考察ax₁、ax₂、……、ax_i、……、ax_φ(n)，因为a和x_i都与n互质，所以ax_i也与n互质。并且，对于任意的0<i,j<φ(n)，有ax_i(mod n)≠ax_j(mod n)，因为假如ax_i≡ax_j(mod n)，那么就有a(x_i-x_j)≡0(mod n)，由于a与n是互质的，所以只能是x_i-x_j≡0(mod n)，但是x_i、x_j都比n小，又不相等，所以x_i-x_j≡0(mod n)不可能成立，所以假设不成立。也就是说ax₁、ax₂、……、ax_i、……、ax_φ(n)也是φ(n)个与n互质的整数，虽然x_i与ax_i不一定模n同余，但是集合{x_i}与集合{ax_i(mod n)}是相同的。即：x₁•x₂•……•x_i•……•x_φ(n)≡ax₁•ax₂•……•ax_i•……•ax_φ(n)(mod n)，所以x₁•x₂•……•x_i•……•x_φ(n)≡a^φ(n)x₁•x₂•……•x_i•……•x_φ(n)(mod n)，所以(a^φ(n)-1)x₁•x₂•……•x_i•……•x_φ(n)≡0(mod n)，所以a^φ(n)-1≡0(mod n)，所以a^φ(n)≡1(mod n)，证毕。

对于欧拉函数，有如下一些性质：

1.如果p是质数，由p数的定义和欧拉函数的定义，显然有φ(p)=p-1。代入欧拉定理不难得到a^p≡a(mod p)也就是费马小定理。

2.如果p是质数，φ(p^k)=p^k-p^k-1，这是因为小于p^k的整数一共有p^k-1个，其中，p和p的倍数p、2p、……、np、……、p²、……、p^k-1、p^k-1+p、……、(p^k-1-1)p一共有p^k-1-1个，所以，φ(p^k)=(p^k-1)-(p^k-1-1)=p^k-p^k-1=p^k(1-1/p)

3.如果m、n互质，那么φ(mn)=φ(m)φ(n)

证明：考察所有小于mn的整数u，设u=am+r(其中r=0,1,……,m-1，a=0,1,……,n-1)，若r为φ(m)个与m互质的数之间的一个，那么u与m互质。也就是说，r_i、m+r_i、2m+r_i、……、(n-1)m+r_i这nφ(m)个数与m互质。在i固定的情况下，这n个数并不一定也与n互质，不过模n的余数互不相同（反证法易证），正好组成模n的完全剩余类，所以其中恰好有φ(n)个与n互质。所以在所有小于mn的整数中，既与m互质，又与n互质的一共有φ(m)φ(n)个，亦即：φ(mn)=φ(m)φ(n)，证毕。

上述证明，构造一个矩阵看上去清楚一点。

1、2、……、r、……、m

m+1、m+2、……、m+r、……、2m

……、……、……、……、……、……

(n-1)m+1、(n-1)m+2、……、(n-1)m+r、……、mn

上述证明的意思是说，这个n行、m列的行列式，其中只有φ(m)列与m互质，而每一个与m互质的列中，又只能有φ(n)列与n互质。

所以，如果n=p₁^m₁•p₂^m₂•……•p_i^m_i，那么φ(n)=φ(p₁^m₁)•φ(p₂^m₂)•……•φ(p_i^m_i)=p₁^m₁(1-1/p₁)•p₂^m₂(1-1/p₂)•……•p_i^m_i(1-1/p_i)=p₁^m₁•p₂^m₂•……•p_i^m_i•(1-1/p₁)(1-1/p₂)……(1-1/p_i)=n(1-1/p₁)(1-1/p₂)……(1-1/p_i)

例如：φ(10)=φ(2×5)=10(1-1/2)(1-1/5)=10×1/2×4/5=4，小于10且与10互质的数是1、3、7、9一共4个。

回到费马小定理，费马小定理又可写作a^p≡a(mod p)当p为质数时成立，所以就有人猜想，满足a^p≡a(mod p)的p就一定是质数（中国猜想），但显然这个猜想不成立，如2³⁴¹≡2(mod 341)但是341=11×31，可以写一个程序来验证，当a和p比较大时，采用传统的连乘算法时，循环几次就变成了两个大数的乘积，很容易溢出，并且关键是最后又要求模。所以这里最好采用快速幂算法。快速幂的原理就是把p转换成b_n2ⁿ+b_n-12^n-1+……+b_i2ⁱ+……+b₀2⁰(其中，b_i=0或1)也就是把求a^p转换成求a^{b_n2n+b_n-12n-1+……+b_i2i+……+b₀20}=a^b_n2n•a^b_n-12n-1•……•a^b_i2i•……•a^b₀20，所以，本来需要算p次乘法，变成了算n次乘法，n≒log₂p，并且，由于b_i=0或1，这n个连乘当中有一些项直接为1，是可以不用计算的。当然a^b_i2i也是一个幂运算，不过，可以从a²⁰开始，把结果保存起来，每一次计算a²ⁱ的时候直接利用a^2i-1的结果，也就是说a²ⁱ=a^2×2i-1=(a^2i-1)²。而把p转换成b_n2ⁿ+b_n-12^n-1+……+b_i2ⁱ+……+b₀2⁰(其中，b_i=0或1)其实就是p的二进制形式，用C语言的话根本不用特殊处理，直接从p的最低位开始，计算a²⁰，判断最低位是否为1，为1就把a²⁰乘到结果里去，为0就跳过，然后右移一位，继续循环，直到p为0为止，程序很好写，代码如下：

#include <stdio.h>
unsigned long long pow(unsigned long long x, long y,long mod) {
    unsigned long long p = 1;
    while (y) {
        if (y & 1)p = x * p%mod;
        x = x * x%mod;
        y >>= 1;
    }
    return p;
}
int prime(long long a) {
    int i;
    if (a == 2)
        return 1;
    for (i = 2; i*i <= a; i++)
        if (a%i == 0)
            return 0;
    return 1;
}
void main() {
    long p,a;
    while (1) {
        scanf("%ld %ld", &p, &a);
        if (a == 0 || p == 0)break;
        if (!prime(p)&&pow(a,p,p)==a)
            printf("yes\n");
        else
            printf("no\n");
    }
}

Sample Input
3 2
10 3
341 2
341 3
1105 2
1105 3
0 0
Sample Output
no
no
yes
no
yes
yes

这个程序是验证中国猜想的，最终是要求模a的余数，所以在快速幂里每一次算乘法都求了模，这样避免中间结果溢出或者太大，但是算法本身并没有要求非要每一步都求模不可。只是在这种模运算的题目当中，这样做能够使可以计算的a、p的范围宽许多。