这篇文章我们将介绍数论当中几个很重要的定理：威尔逊定理、费马小定理以及欧拉定理，并讨论一些基于这些定理的算法。
首先我们给出费马小定理：如果p是素数，并且gcd(a,p) = 1 , 那么有a^(p-1) = 1(mod p)。
而关于这个定理的证明，也是不难理解的。
在证明之前，我们先需要知道这样两个引理。
引理1：若a、b、c为任意3个整数，且gcd(m,c) = 1,则当ac = bc(mod m)时，有a = b (mod m)。
引理2：设m是一个整数,且m>1,b是一个整数且gcd(m,b)=1.如果a1,a2,a3,a4,…am是模m的一个完全剩余系,那么ba[1],ba[2],ba[3],ba[4],…ba[m]也构成模m的一个完全剩余系。
这里需要解释一下的是，所谓完全剩余系，即{a1,a2,a3,a4,…am}中任意一个元素ai % m != 0。针对引理1，通过同余式本身的特点，都十分易证，随后再基于引理1，引理2也不难得证。这里便不再累述。 有了这两个引理做铺垫，下面便可以开始费马小引理的证明了。
构造素数p的完全剩余系P = {1,2,3,……p-1}。
因为gcd(a,p) = 1,由引理2得素数p的新的完全剩余系：A = {a,2a,3a,……(p-1)a}。
我们选出A、P中第i个元素，由引理1得，Ai = Pi(mod p)，因此我们容易得到如下等式。
1 * 2 * 3 *……(p-1) = a * 2a * 3a * ……(p-1)a (mod p)
等式两边同除(p - 1)!，得到a^(p-1) = 1 (mod p)，定理得证。
既然知道了费马小定理的内容和证明了，这里我们就用它进行一下拓展利用。
首先我们知道，在费马小定理的等式a^(p - 1) = 1 (mod p)中，如果p是素数，那么等式是成立的。
即 isprime(p) = 1 => a^(p - 1) = 1 (mod p) 。 那么我们根据逆否命题与原命题的等价性，可以得到，a^(p - 1) != 1 (mod p) => isprime(p) = 0。
那么逆命题和否命题与原命题是否等价呢？如果等价，我们岂不是可以通过判断是否满足费马小定理的等式来判断p是否为素数了？
实践表明，满足费马小定理的a、p组合中，p是可以为合数的，不妨自己举几个例子。 因此对于满足费马小定理等式的a、p，有这样的定义。如果p是合数，则称p是以a为基的伪素数。
我们不妨通过一个题目来理解一下所谓伪素数的概念。(Problem source:pku 3641)
Description
Fermat's theorem states that for any prime number p and for any integer a > 1, ap = a (mod p). That is, if we raise a to the pth power and divide by p, the remainder is a. Some (but not very many) non-prime values of p, known as base-a pseudoprimes, have this property for some a. (And some, known as Carmichael Numbers, are base-a pseudoprimes for all a.)
Given 2 < p ≤ 1000000000 and 1 < a < p, determine whether or not p is a base-a pseudoprime.
Input
Input contains several test cases followed by a line containing "0 0". Each test case consists of a line containing p and a.
Output
For each test case, output "yes" if p is a base-a pseudoprime; otherwise output "no".
题目大意：基于伪素数的概念，给出a、p，请你判断p是否是以a为基的伪素数。
数理分析：有了上文对伪素数概念的较少，我们知道，伪素数首先必须是合数，然后需要满足费马小定理。
编程实现：在判断是否是伪素数的时候，由于是判断单个数字，所以使用朴素的判断法即可。而判断是否满足费马小定理的时候，则用到了我们在之前讨论过的乘法快速幂算法。 参考代码如下。
#include<iostream> #include<math.h> using namespace std; bool ifprime(long long n) { long long i; int iffind = 0; if(n == 2 || n == 1) return 1; else { for(i = 2;i <= sqrt(n + 0.0);i++) if(n%i == 0) { iffind = 1; break; } if(iffind) return 0; else return 1; } } long long quick_mod(long long a , long long b , long long m) //快速幂取模，计算a^b(mod m)的值 { long long ans = 1; while(b) { if(b&1) { ans = (ans * a)%m; b--; } b /= 2; a = a * a %m; } return ans; } int main() { long long n , a, p ,b; while(cin >> p >> a) { if(p == 0 && a == 0) break; if(ifprime(p)) cout<<"no"<<endl; else { long long ans = quick_mod(a,p,p); if(ans == a) cout<<"yes"<<endl; else cout<<"no"<<endl; } } return 0; }
通过上文对费马小定理初步的理解，我们这里介绍在其基础上建立的一种更加高效的素数检测法——Miller-Rabin素数测试法。
通过上面的介绍，我们已经知道，满足费马小定理仅仅是一个数是素数的必要条件，那么对于一个满足费马小定理的数字，我们再通过什么方法来判断它是素数还是合数呢？
可能有人会想到之前学过的筛法，但是还是不够高效，在判断较大的数是否是素数的时候，只要一用筛法，时间复杂度就会激增，而避免这一情况的发生的办法就是想出线性时间复杂度的算法。
这里我们给出数论中一个线性时间复杂度判断素数的定理。
二次探测定理：如果p是一个素数，且0<x<p，则方程x^2 = 1(mod p)的解为x = 1或x = p - 1。
通过这个名字我们也可以看到，好像这个定理就是为了弥补利用费马小定理检测素数的诞生的。那么我们可以简单的概括这个新的检测素数的方法了。
step1：对于整数p，利用费马小定理判断。如果不满足，输出结果，如果满足，进入step2。
step2：利用二次探测定理判断，输出结果。 这就是基于费马小定理的改进版的Miller-Rabin素数检测法。由于两个定理中都包含未知量，所以再编程实现的时候我们都需要随机生成符合要求的参量，因此,Miller-Rabin算法本质上是一种随机算法。
那么下面我们开始讨论具体实现该算法的一些细节的处理。
虽然我们在描述算法步骤的时候，似乎将两个定理的应用分开了，但是通过观察两个定理的形式，会发现其实有着相似点，我们尝试将两个定理结合在一起。
假设我们想要判断n是否为素数，并随机生成了整数a(为了满足应用二次探测定理的条件),基于费马小定理，我们想要判断a^(n-1) = 1 (mod n)是否成立。
为了将此表达式更近得像二次探测定理靠拢，我们要找到含有平方的形式。那么我们就将n - 1化成 d * 2^r的形式(此处d为奇数)，这样费马小定理的左式就可以化成平方的形式了。
现在判断a^(d*2^r) = 1 (mod n)是否成立，如果成立，n是素数。 即(((a^d)^2)^2)^2)... = 1(mod n) ①
此时左式有r + 1层括号，并且这个等式同时包含了两个定理的形式。我们从最内层的a^d分析。
1.如果a^d = 1(mod n)，则①显然成立，n是素数。
2.如果a^d ≠ 1 (mod n)，结合二次探测定理，会有如下两种情况。
1'.a^d = n - 1,则(n-1)^2 = 1(mod n)，这就引导出了情况1.，n为素数。
2'.a^d ≠n - 1，等式不成立， n不是素数。
以上便完成了二次探测定理和费马小定理的巧妙结合。 我们通过一个题目来具体实现以下改进版的Miller-Rabin算法。(Problem source : hdu 2138)
题目大意：给出你一些数，让你判断素数的个数。
数理分析：由于题目给出的素数的范围很大，因此筛法就显得太过耗时，因此这里就可以用Miller-Rabin算法了。
编程实现：结合上文给出的分析，利用右移位运算来实现n - 1 = d*(2^r)的转换，并且还需要我们曾经讨论过的乘法快速幂的算法。
参考代码如下。
#include<stdio.h> #include<stdlib.h> #include<cmath> long long quick_mod(long long a,long long b,long long m) { long long ans = 1; while(b) { if(b&1) { ans = (ans * a) % m; b--; } b /= 2; a = a * a % m; } return ans; } bool witness(long long a,long long n) { long long m = n - 1; int j = 0; while(!(m&1)) //进行n-1 = d*(2^r)的转化 { j++; m >>= 1; } long long x = quick_mod(a, m , n); if(x == 1 || x == n -1) return false; while(j--) { x = x * x % n; if(x == n-1) return false; } return true; } bool miller_rabin(__int64 n) { if(n==2) return true; if(n==1 || ((n&1)==0)) return false; for(int i=0;i<50;i++)//随机生成次数为50次，该次数怎么得来的涉及一些概率方面的实验，在这里不作展开论述，会在以后关于随机算法的深入讨论中分析。 { __int64 a=rand()*(n-2)/RAND_MAX +1; //生成随机数a保证a<n if(witness(a, n)) return false; } return true; } int main() { int n,cnt; __int64 a; while(scanf("%d",&n)!=EOF) { cnt=0; while(n--) { scanf("%I64d",&a); if(miller_rabin(a)) cnt++; } printf("%d ",cnt); } return 0; }
基于上文对费马小定理的引入及其相关周边概念的介绍，我们再来通过一个题目来看看费马小定理在解决数论问题中有着怎样的应用。(Problem source : hdu 4704)
题目大意：定义函数s(k)，表示一个整数n被分成k个数的方案数，现在给定n，求解∑s(i) mod (1e9 + 7)的值，其中i∈[1,n]。
数理分析：这是一道比较综合的数论题目，其中涉及了很多的数论知识点。下面我们依次来击破。
首先我们分析函数s(i)的值，其功能是给出整数n被分成i个数字的方案数，我们不妨将整数n写成1+1+……+1的形式，这n个1中间形成了n-1个空间，我们插入i-1个隔板，那么我们整数n就被分成了i份，显然s(i) = C(n-1,i-1),这里运用到了组合数学中的鸽巢定理。
所以，我们所需要求的变成了∑C(n-1,i-1) , i∈[1,n]，由组合数学中二项式的相关性质，我们容易知道上式等于2^(n-1)。
但是基于n的取值太大，我们要考虑简化方案。
记m = 1e9 + 7。
当前我们所求的是：2^(n-1) = ? (mod m)。
通过费马小定理我们看到，a^(p-1) = 1(mod p)，当gcd(a,p)=1的时候。容易看到，gcd(a , m) = 1，那么有2^(m - 1) = 1(mod m)。
我们假设有n-1 = x(m-1) + q,那么基于同余等式的基本性质，2^(n-1) = 2^(x*(m - 1)) * 2^q (mod m) = 2^(x*(m-1)) (mod m) * 2^q(mod m) = 2^q(mod m)。
这样，一个很大的指数n-1，以上的等价转换，幂次就大大降低了，随后用快速幂取模来实现2^p (mod m)来得到最终结果即可。
总结起来，这道题目综合了组合数学的二项式的性质、鸽巢原理、费马小定理、快速幂取幂等多个算法，是一个很综合的数论题目。
参考代码如下。
#include<stdio.h> using namespace std; const int MAX = 100000 + 10; const int mod = 1000000000 + 7; char s[MAX]; long long MOD(char *a , int Mod) { long long sum = 0; for(int i = 0;a[i]!=' ';++i) { sum=(10*sum + a[i] - '0')%Mod; } return sum; } long long quick_mod(long long a,long long b,long long m)//快速幂取模：求a^b%m { long long ans = 1; while(b) { if(b&1) { ans = (ans * a) % m; b--; } b /= 2; a = a * a % m; } return ans; } int main() { while(scanf("%s",s) != EOF) { long long n = MOD(s , mod - 1) - 1; printf("%I64d ",quick_mod(2,n,mod)); } }
我们你再来看一道利用费马小定理在具体的数论问题中的应用。(Problem source : hdu 1098)
题目大意：给出了f(x)的形式，给出k，让你求解最小的a，使得对于任意的整数x，都有65|f(x)。
数理分析：首先我们对f(x)做简单的因式分解。
f(x)=5*x^13+13*x^5+k*a*x
=x(5*x^12 + 13*x^4 + k*a)
对于任意的x，要使得f(x)整除65，那么我们只好使得后面的因子——(5*x^12 + 13*x^4 + k*a)整除65恒成立。
即：5*x^12 + 13*x^4 +k*a = 0(mod 65)，容易看到，65 = 13 * 5，而这个式子中恰好有5和13，那么我们就尝试找到一些联系。对于这个式子能够整除65，我们不妨将其转换既能整除5，又能整除13。即这个问题被分解成两个子问题：
1. 5*x^12 + 13*x^4 + k*a = 0 (mod 13)
2. 5*x^12 + 13*x^4 + k*a = 0(mod 5)
基于这种形式，我们显然可以继续转化。
1'. 5*x^12 + k*a = 0 (mod 13)
2'. 13*x^4 + k*a = 0 (mod 5)
针对1'. , 有gcd(13,12) = 0 , 因此5*x^12 = 5(mod 13)，因此当且仅当k*a = 8(mod 13)，这种情况是成立的，此多项式是13的整数倍。
针对2'. ,同理分析，当且仅当k*a = 2(mod 5)，此多项式是5的整数倍。
这两种综合起来，两个条件同时成立，即这个多项式是65的整数倍，对于任意的x，该函数值f(x)都是能够整除65的。
通过这个问题我们不难看出，费马小定理不但能够在特定的情况下做降幂的恒等变化，而且在解决整除性问题在简化问题也有着用武之地。
参考代码如下。
#include<stdio.h> using namespace std; int main() { int k , i , flag; while(scanf("%d",&k)!=EOF) { for(i=1;i<66;i++) { if(i*k%13==8 && i*k%5==2) { flag = i; break; } else flag = 0; } if(flag) printf("%d ",flag); else printf("no "); } }
我们再来看一道利用费马小定理进行优化计算的题目。(Problem soure:hdu 4549)
F[0] = a F[1] = b F[n] = F[n-1] * F[n-2] ( n > 1 )
现在给出a, b, n，你能求出F[n]的值吗？
我们通过写出F[n]的前几项，会发现F[n]其实和斐波那契数列fib[n]联系起来，即F[n] = a^fib[n-1] * b^fib[n]。记mod = 1000000007.
因此F[n] % m = (a^fib[n-1] * b^fib[n]) mod m
= (a^fib[n-1] mod mod) * (b^fib[n] mod mod)
观察到mod的取值和a、b的取值范围，能够看到gcd(mod,a) = gcd(mod ,b) = 1,并且mod是素数，由此我们利用费马小定理可知：
a^(mod - 1) ≡ 1(% mod)
b^(mod - 1)≡ 1(% mod)
基于这两个等式便可很好的简化计算，随后利用矩阵快速幂求fib[n]即可。
参考代码如下。
#include<cstdio> const long long mod = 1000000007; struct Matrix { long long mat[2][2]; }; const Matrix P = //得到斐波那契数列的元矩阵 { 1 , 1, 1 , 0, }; const Matrix I = //单位矩阵 { 1 , 0 , 0 , 1 , }; Matrix matrixmul(Matrix a , Matrix b) { Matrix c; int i , j , k; for(i = 0;i < 2;i++) for(j = 0;j < 2;j++) { c.mat[i][j] = 0; for(k = 0;k < 2;k++) c.mat[i][j] += a.mat[i][k]*b.mat[k][j]; c.mat[i][j] %= (mod - 1);//费马小定理 } return c; } Matrix quickpow(long long n) //矩阵快速幂 { Matrix m = P , ret = I; while(n) { if(n&1) ret = matrixmul(ret , m); n >>= 1; m = matrixmul(m,m); } return ret; } long long quickpow(long long a , long long b) //a^b % mod { long long ret = 1; while(b) { if(b&1) ret = ret*a%mod; b >>= 1; a = a * a % mod; } return ret; } int main() { long long a , b , n; Matrix q; while(~scanf("%I64d%I64d%I64d",&a,&b,&n)) { q = quickpow(n); printf("%I64d ",quickpow(a,q.mat[1][1])*quickpow(b,q.mat[1][0])%mod); }//( a^(f[n-1]%mod) * b^(f[n]%mod) ) % mod return 0; }
威尔逊定理：如果p是素数的充分必要条件是(p-1)! ≡ -1(mod p)。
关于这条定理的详细证明，笔者将在《<具体数学>——数论》中介绍，这里我们直接来看一道应用的题目。(Problem source : hdu 5391)
Problem Description
Tina Town is a friendly place. People there care about each other.
Tina has a ball called zball. Zball is magic. It grows larger every day. On the first day, it becomes
time as large as its original size. On the second day,it will become
times as large as the size on the first day. On the n-th day,it will become
times as large as the size on the (n-1)-th day. Tina want to know its size on the (n-1)-th day modulo n
题目大意：给出变量n，求解(n-1)! mod n的结果。
数理分析：结合威尔逊定理，我们通过判断n是否为素数进行分别讨论。
当n是素数，显然结果是-1，考虑到取余结果一般不考虑负数，这里将其看成n-1.
当n是合数，我们容易看到n(n-1)!，但是n = 4时是一种特殊情况，4 = 2*2，而3！中只有一个因子2，所以并不满足n(n-1)!。
简单的参考代码如下。
#include<cstdio> #include<cmath> using namespace std; bool isprime(int m) { if(m == 2) return false; else if(m % 2 == 0) return false; else { for(int i = 3;i <= (int)sqrt(m);++i) { if(m % i == 0) return false; } } return true; } int main() { int n , t; scanf("%d",&t); while(t--) { scanf("%d",&n); if(n == 4) printf("2 "); else if(isprime(n)) printf("%d ",n-1); else printf("0 "); } }
我们再来看一道利用威尔逊定理进行化简运算的题目。(Problem source : hdu 2973)
However, the test turned out difficult for some math PhD students and even for some professors. Therefore, the math department wants to write a helper program which solves this task (it is not irrational, as they are going to make money on selling the program).
The task that is presented to anyone visiting the start page of the math department is as follows: given a natural n, compute
#include<cstdio> using namespace std; const int N = 3000000 + 50; int num[N]; bool notprime[N]; void Init() { for(int i = 2;i < N;i++) { if(notprime[i] == false) { if( i > 7 && (i-7)%3 == 0) num[(i-7)/3]++; for(int j = i + i;j < N ;j += i) notprime[j] = true; } } for(int i = 1;i < N;i++) num[i] += num[i-1]; } int main() { Init(); int t; scanf("%d",&t); while(t--) { int n; scanf("%d",&n); printf("%d ",num[n]); } }