3.1 素数运算的实验范例
素数又称质数，指的是在大于1的自然数中除了1和自身外无法被其他自然数整除的数。简言之，只有两个正因数（1和自身）的自然数即为素数。比1大但不是素数的数称为合数；1和0既非素数也非合数。合数是由若干个素数相乘而得到的，所以素数是合数的基础，没有素数就没有合数。整数的基本元素是素数，数论作为研究整数性质的一门理论，其本质就是对素数性质的研究。
下面，我们展开素数运算的两个实验：
1）计算整数区间［2..n］中的所有素数。
2）大整数的素数测试。
3.1.1 使用筛法生成素数的实验范例
如何计算整数区间［2..n］中的所有素数？我们先介绍一种最为简便的筛法——埃拉托斯特尼筛法：
设筛子u［］，初始时区间中的所有数都在筛中。按递增顺序搜索筛中的最小数，将其倍数从筛中筛去，最终筛中留下的数即为素数。int i,j,k;
for(i=2;i<=n;i++)u［i］=true; // 初始时所有数都在筛中
for(i=2;i<=n;i++)// 顺序搜索筛中的最小数
if(u［i］){
for(j=2;j*i<=n;j++) // 将i的倍数从筛中筛去
u［j*i］=false;
}
for(i=2;i<=n;i++)if(u［i］){// 将筛中的所有素数放入su［］中
su［++num］=i;
}
上述算法的时间复杂度为O(n*long logn)。算法中合数是作为素数的倍数被筛去的。显然，如果每个合数仅被它最小的质因数筛去，则算法效率可以大幅提升。由此引出一种优化的算法——欧拉筛法：
int i,j,num=1;
memset(u,true,sizeof(u));
for(i=2;i<=n;i++){// 顺序分析整数区间的每个数
if(u［i］) su［num++］=i;// 将筛中最小数送入素数表
for(j=1;j<num;j++) {// 搜索素数表的每个数
if (i*su［j］>n)break;// 若i与当前素数的乘积超出范围，则分析下一个整数i
u［i*su［j］］=false;// 将i与当前素数的乘积从筛子中筛去
if (i% su［j］==0) break;// 若当前素数为i的最小素因子，则分析下一个整数i
}
}欧拉筛法的时间复杂度可优化至O(n)。证明如下：
设合数n的最小素因子为p，它的另一个大于p的素因子为p′，令n=pm=p′m′。观察上面的程序片段，可以发现j循环到素因子p时合数n第一次被标记（若循环到p之前已经跳出循环，说明n有更小的素因子）。若也被p′标记，则是在这之前（因为m′素数筛选法通常作为数论运算的核心子程序。我们不妨来看几个实例。
【3.1.1.1 Goldbach’s Conjecture】
【问题描述】
在1742年，德国的一个业余数学家Christian Goldbach给Leonhard Euler写信，在信中给出如下猜想（哥德巴赫猜想）：
每个大于4的偶数都可以写成两个奇素数的和。例如：8=3+5，3和5都是奇素数；而20=3+17=7+13；42=5+37=11+31=13+29=19+23。
现在哥德巴赫猜想仍然没有被证明是否正确。现在请证明对所有小于1000000的偶数，哥德巴赫猜想成立。
输入：
输入包含一个或多个测试用例。每个测试用例给出一个偶整数n，6≤n<1000000。输入以0结束。
输出：
对每个测试用例，输出形式为n=a+b，其中a和b是奇素数，数字和操作符要用一个空格分开，如样例输出所示。如果有多于一对的奇素数的和为n，就选择b-a最大的一对。如果没有这样的数对，则输出"Goldbach’s conjecture is wrong."。
试题来源：Ulm Local 1998
在线测试地址：POJ 2262，ZOJ 1951，UVA 543
试题解析
先离线计算［2..1000000］的素数表su［］和素数筛u［］。然后每输入一个偶整数n，顺序搜索su中的每个素数(2*su［i］≤n)：若整数n-su［i］亦是素数（u［n-su［i］］=true），则su［i］和n-su［i］为满足条件的数对。
程序清单
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<cstdio>
using namespace std;
bool u［1111111］; // 筛子
int su［1111111］,num;// 素数表为su［］，该表长度为num
void prepare(){// 使用筛选法构建素数表su［］
int i,j,num;
memset(u,true,sizeof(u));
for(i=2;i<=1000000;i++){// 顺序分析整数区间的每个数
if(u［i］) su［++num］=i;// 将筛中最小数送入素数表
for(j=1;j<=num;j++) {// 搜索素数表的每个数
if (i*su［j］>n)break;
u［i*su［j］］=false;// 将i与当前素数的乘积从筛子中筛去
if (i% su［j］==0) break;// 若当前素数为i的素因子，则分析下一个整数
}
}
}
int main () {
prepare();// 使用筛选法构建素数表su［］
int i,j,k,n;
while(scanf("%d",&n)>0&&n)// 反复输入偶整数，直至输入0为止
{
bool ok=false;
for(i=2;i<=num;i++)// 按照递增顺序搜索素数表中的每个素数
{
if(su［i］*2>n)break;// 搜索完所有素数和的形式
if(u［n-su［i］］){// 若n能够拆分出两个素数和的形式，则成功退出
ok=true;
break;
}
}
if(!ok)puts("Goldbach's conjecture is wrong.");// 输出结果
else printf("%d = %d + %d＼n",n,su［i］,n-su［i］);
}
return 0;
}【3.1.1.2 Summation of Four Primes】
【问题描述】
Euler证明的经典定理之一是素数在数量上是无限的。但每个数字是否可以表示成4个素数的总和呢？我不知道答案，请你来帮我。
输入：
在输入的每行中给出一个整数n (n≤10000000)，请将这一整数表示为4个素数的总和。输入以EOF为结束。
输出：
对于输入的每一行，给出一行4个素数的输出。如果给出的数字不能表示为4个素数的总和，则输出一行“Impossible.”。存在多个解的情况，任何成立的解答都会被接受。
试题来源：Regionals 2001 Warmup Contest
在线测试地址：UVA 10168
试题解析
本题是哥德巴赫猜想的一个扩展。计算方法如下：
先采用筛选法计算［2..9999999］的素数表su［］，表长为num。在离线计算出su［］的基础上，通过直接查表计算每个n的分解方案：
1)先直接推算出n≤12的分解方案：
n<8，则n不能表示为4个素数的总和。
n==8，则n分解出2 2 2 2。
n==9，则n分解出2 2 2 3。
n==10，则n分解出2 2 3 3。
n==11，则n分解出2 3 3 3。
n==12，则n分解出3 3 3 3。
2)在n>12的情况下，先分解出头两个素数：
若n为偶数(n%2==0)，则两个素数为2 2 ，n-=4；
若n为奇数，则两个素数为2 3 ，n-=5。
显然，剩余的n为大于4的偶数。依据哥德巴赫猜想（每个大于4的偶数可以写成两个奇素数的和），顺序搜索su［］中的每个素数（1≤i≤num，2*su［i］≤n）：若su［n-su［i］］=true，则说明剩余的n又分解出后两个素数su［i］和n-su［i］，成功退出。
程序清单
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
using namespace std;
bool u［10000001］; // 筛子
int su［5000000］,num;// 素数表及其表长
void prepare(){// 使用筛选法构建［2..9999999］的素数表su［］
int i,j,num;
memset(u,true,sizeof(u));// 初始时所有数在筛中
for(i=2;i<=9999999;i++){// 顺序分析整数区间的每个数
if(u［i］) su［++num］=i;// 将筛中最小数送入素数表
for(j=1;j<=num;j++) {// 搜索素数表的每个数
if (i*su［j］>n)break;// 若i与当前素数的乘积超出范围，则分析下一个整数
u［i*su［j］］=false;// 将i与当前素数的乘积从筛子中筛去
if (i% su［j］==0) break;// 若当前素数为i的素因子，则分析下一个整数
}
}
}
int main ()
{
prepare();// 使用筛选法构建素数表su［］
int n,i,j,k;
while(scanf("%d",&n)>0){// 输入整数n
if(n==8){puts("2 2 2 2");continue;}// 输出特例
if(n==9){puts("2 2 2 3");continue;}
if(n==10){puts("2 2 3 3");continue;}
if(n==11){puts("2 3 3 3");continue;}
if(n==12){puts("3 3 3 3");continue;}
if(n<8){puts("Impossible.");continue;}
if(n%2==0){printf("2 2 ");n-=4;}// 先分离出头两项
else{printf("2 3 ");n-=5;}
for(i=1;i<=num;i++)// 按照递增顺序枚举素数
{
if(su［i］*2>n)break;// 若无法产生另两项素数，则退出循环
if(u［n-su［i］］){// 若su［i］与n-su［i］为另两项素数，则输出
printf("%d %d＼n",su［i］,n-su［i］);
break;
}
}
}
}【3.1.1.3 Digit Primes】
【问题描述】
一个素数是一个能被两个不同的整数整除的正整数。一个位素数（digit prime）是一个所有的位数相加的总和也是素数的素数。例如素数41是一个位素数，因为4+1=5，而5是一个素数。17不是位素数，因为1+7=8，8不是素数。本题要求在1000000的范围内，给出位素数的数量。
输入：
输入的第一行给出一个整数N (0输出：
除了第一行外，对输入的每一行输出一行，给出一个整数，表示在t1和t2之间（包含t1和t2）的位素数的个数。
注意：本题输入和输出要用scanf()和printf()，cin和cout太慢，会导致超时。
试题来源：The Diamond Wedding Contest: Elite Panel’s 1st Contest 2003
在线测试地址：UVA 10533
试题解析
设：
［2..1100001］的素数筛为u［］；
前缀的位素数个数为u2［］，其中u2［i］为［2..i］中位素数的个数(2≤i≤1100001)。显然，区间［i，j］中位素数的个数为u2［j］-u2［i-1］ (2≤i≤j≤1100001)。
我们先离线计算出u2［］，计算方法如下:
1)采用筛选法计算出［2，1100001］的素数筛u［］。
2)计算［2，1100001］中的每个位素数i，u2［i］=1∣u［i］&&u［i的数位和］=true。
3)递推前缀的位素数个数u2［］：u2［i］+=u2［i-1］ ( 2≤i≤1100001)。
借助u2［］表，便可以直接计算任意区间［i，j］中的位素数个数（u2［j］-u2［i-1］）。
程序清单
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
using namespace std;
bool u［1100001］;// 素数筛
int u2［1100001］;// 前缀的位素数个数
void prepare(){// 计算［2..1100001］的素数筛u［］
int i,j,k;
for(i=2;i<1100001;i++)u［i］=1;// 初始时所有数在素数筛
for(i=2;i<1100001;i++)// 取出筛中最小数，将其倍数从筛中筛去
if(u［i］)
for(j=i+i;j<1100001;j+=i)
u［j］=false;
}
bool ok(int x){// 返回x的所有数位和为素数的标志
int i,j,k=0;
while(x){// 计算x的数位和
k+=x%10;x/=10;
}
return u［k］;
}
int main (){
int i,j,k;
prepare();// 计算［2‥1100001］的素数筛
for(i=2;i<1100001;i++)// 计算［2..109999］中的所有位素数
if(u［i］)&&(ok(i)) u2［i］=1;
for(i=2;i<1100001;i++)u2［i］+=u2［i-1］;// u2［i］为［2..i］中位素数的个数
scanf("%d",&k);// 输入区间的个数
while(k--){
scanf("%d %d",&i,&j);// 输入当前区间［i,j］
printf("%d＼n",u2［j］-u2［i-1］);// 输出［i,j］中位素数的个数
}
}【3.1.1.4 Prime Gap】
【问题描述】
在两个相继的素数p和p+n之间，n-1个连续合数（composite number，不是素数且不等于1的正整数）组成的序列，被称为长度为n的素数间隔（prime gap）。例如，在23和29之间长度为6的素数间隔是<24, 25, 26, 27, 28>。
给出一个正整数k，请编写一个程序，计算包含k的素数间隔的长度。如果没有包含k的素数间隔，则长度为0。
输入：
输入由一个行序列组成，每行一个正整数，每个正整数大于1、小于或等于第100000个素数，也就是1299709。以包含一个0的一行标志输入结束。
输出：
输出有若干行，每行给出一个非负整数，如果相应的输入整数是一个合数，则输出素数间隔的长度；否则输出0。输出中没有其他字符出现。
试题来源：ACM Japan 2007
在线测试地址：POJ 3518，UVA 3883
试题解析
设ans［k］为包含k的素数间隔的长度。显然，若k为素数，则ans［k］=0；若k为合数且k位于素数p1和p2的中间，则k所在合数区间内每个合数的ans值都为同一个数，即ans［p1+1］=ans［p1+2］…=ans［p2-1］=(p2-1)-(p1+1)+2。由此得出算法：
1)通过下述方法计算ans［］：
①采用筛选法计算［2..1299709］的素数筛u［］。
②顺序枚举［2..max n-1］中的每个数i：若i是素数(u［i］=true)，则ans［i］=0；否则计算i右邻的素数j（u［i］=u［i+1］=…=u［j-1］=false，u［j］=true），置ans［i］=ans［i+1］=…=ans［j-1］=j-i+2，并设i=j，以提高枚举效率。
2）在离线计算出ans［］基础上，每输入一个整数k，则包含k的素数间隔的长度即为ans［k］。
程序清单
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
using namespace std;
const int maxn=1299710;// 整数值的上限
bool u［maxn］;// 素数筛
int ans［maxn］;// 包含每个整数的素数间隔长度
void prepare(){
int i,j,k;
for(i=2;i<maxn;i++)u［i］=1;// 使用筛选法计算［2..1299709］的素数筛u［］
for(i=2;i<maxn;i++)
if(u［i］)// 若i为素数，则将其倍数从筛中筛去
for(j=2;j*i<maxn;j++) u［i*j］=0;
for(i=2;i<maxn;i++)// 枚举［2.. maxn-1］中的每个数
if(!u［i］){// 若i是合数，则计算合数区间［i，j］
j=i;
while(j<maxn&&!u［j］) j++;
j--;
for(k=i;k<=j;k++) ans［k］=j-i+2;// 置合数区间内每个合数的ans值
i=j;
}else ans［i］=0;// 素数的ans值为0
}
int main ()
{
int i,j,k;
prepare();// 使用筛选法计算［2..1299709］的素数筛u［］
while(scanf("%d",&k)>0&&k>0){// 反复输入整数k，直至输入0为止
printf("%d＼n",ans［k］);// 输出包含k的素数间隔的长度
}
}
3.1.2 测试大素数的实验范例
解决素数测试问题最简便的方法有试除法，即试用2..n去除n。n是素数则当且仅当没有一个试用的除数能被n整除。但缺点是该算法的时效取决于n。如果n很小，试除法可以在短时间内出解。如果n较大或者很大时，判断n是否为素数则需要花费很多时间。有没有减少测试时间的办法呢？有。下面介绍两种优化算法：
1）筛选法+试除法。
2）Miller_Rabin方法。
（1）筛选法+试除法
如果整数x的上限n比较大，不妨采用筛选法+试除法来提高运算时效：
先通过筛选法计算［2..n］的素数筛u［］和素数表su［1］…su［num］。x是素数当且仅当x为［2..n］中的一个素数（u［x］=1）或者x不能被su［］表中的任何素数整除（x%su［1］≠0，…，x%su［num］≠0）。时间复杂度为O(n)。
【3.1.2.1 Primed Subsequence】
【问题描述】
给出一个长度为n 的正整数序列，一个素序列（primed subsequence）是一个长度至少为2 的连续的子序列，总和是大于或等于2的一个素数。例如给出序列：3 5 6 3 8，存在两个长度为2的素序列（5+6=11以及3+8=11）、一个长度为3的素序列（6+3+8=17）和一个长度为4的素序列（3+5+6+3=17）。
输入：
输入包含若干测试用例。第一行给出一个整数t（1输出：
对每个序列，输出“Shortest primed subsequence is length x:”，其中x是最短的素序列的长度，然后给出最短素序列，用空格分开。如果操作多个这样的序列，输出第一个出现的序列；如果没有这样的序列，输出“This sequence is anti-primed.”。
试题来源：June 2005 Monthly Contest
在线测试地址：UVA 10871
试题解析
由于序列长度n的上限为10000，而序列中每个非负整数的上限为10000，因此首先需要解决的一个问题是如何快捷地判断素序列，即判断序列中若干元素的和x为素数。
我们首先使用筛选法，离线计算出［2..10010］的素数表su［］和素数筛u［］，su［］表的长度为num。若x为［2..10010］中的一个素数（u［x］=1）或者x不能被su［］表中的任何素数整除（x%su［1］≠0，…，x%su［num］ ≠0），则表明x为素数。
在离线计算出区间［2..10010］中素数的基础上，展开素序列的计算：输入长度为n的序列，递推序列中前i个正整数的和S［i］（1≤i≤n，s［i］+=s［i-1］）；
使用动态规划方法计算最短的素序列：
枚举长度i（2≤i≤n）：
枚举首指针j(1≤j≤n-i+1)：
If（s［i+j-1］-s［j-1］）为素数）
输出序列中第j.. j+i-1个整数（s［j+k-1］-s［j+k-2］，1≤k≤i）并退出程序；
输出失败信息；
程序清单
#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
using namespace std;
bool u［10010］; // 素数筛
int su［10010］,num;// 素数表及其长度
void prepare(){// 使用筛选法构建［2..10010］的素数表su［］
int i,j,num;
memset(u,true,sizeof(u));
for(i=2;i<=10010;i++){// 顺序分析整数区间的每个数
if(u［i］) su［++num］=i;// 将筛中最小数送入素数表
for(j=1;j<=num;j++) {// 搜索素数表的每个数
if (i*su［j］>n)break;
u［i*su［j］］=false;// 将i与当前素数的乘积从筛子中筛去
if (i% su［j］==0) break;// 若当前素数为i的素因子，则分析下一个整数
}
}
}
bool pri(int x){// 若x在小于10010的情况下为素数，或者x在不
// 小于10010的情况下不能被su［］表中的任何素数整除，则返回true；否则返回false
int i,j,k;
if(x<10010)return u［x］;
for(i=1;i<=num;i++)
if(x%su［i］==0)return false;
return true;
}
int n,s［10010］;// 序列中前i个正整数的和为S［i］
int main()
{
int i,j,k;
prepare();// 离线计算素数表su［］
int te;
scanf("%d",&te);// 输入测试用例数
while(te--){
scanf("%d",&n);// 输入序列长度
s［0］=0;
for(i=1;i<=n;i++)// 输入n个整数，计算前缀和
{
scanf("%d",&s［i］);
s［i］+=s［i-1］;
}
bool ok=false;
for(i=2;i<=n;i++){// 枚举长度
for(j=1;j+i-1<=n;j++)// 枚举首指针
{
k=s［i+j-1］-s［j-1］;// 计算第j～j+i-1个整数的和
if(pri(k)){// 若k为素数或者k不能被任何素数整除，则第
// j～i+j-1个整数组成素序列，输出并成功退出
ok=true;
printf("Shortest primed subsequence is length %d:",i);
for(k=1;k<=i;k++)printf(" %d",s［j+k-1］-s［j+k-2］);
puts("");
break;
}
}
if(ok)break;
}
if(!ok)puts("This sequence is anti-primed."); // 若不存在素序列，则返回失败信息
}
// system("pause");
}如果O(n)的时间复杂度仍未达到预期，则不妨采用另一种简便的素数测试方法——Miller_Rabin方法。
(2)Miller_Rabin方法
测试大素数的Miller_Rabin方法基于下述原理：
如果n是素数，且与a 互质，则an-1≡1 (mod n) （1≤a≤n）。
由于a≡a(mod n)，因此根据同余式相乘的性质，an≡a mod n。显然，上述原理与费马小定理an≡a mod n如出一辙（见3.3.1节）。
证明：
要证明an-1≡1(mod n)，只要证明费马小定理an≡a mod n即可。使用数学归纳法证明：
1）当a=1时，由于1n-1≡1，结论显然成立。
2）假设当a为某正整数时, an≡a mod n，我们考虑(a+1)n是否与(a+1)同余。(a+1)n=an+n
1an-1+n
2an－2+…+n
n－1a1+1等号右边中n
1an-1+n
2an－2+…+n
n－1a1的每一项都是 n的倍数，因此(a+1)n≡(an+1)mod n成立。又由于我们假设an≡a mod n ，所以由该式可得 (a+1)n≡(a+1)mod n，即an-1≡1(mod n)成立。
但这里需要指出的是，费马小定理并不是n为素数的充分必要条件，即若n是素数且a与之互质的话，则一定满足an-1≡1(mod n)；但反过来，满足an-1≡1(mod n)的n并不一定是素数。
若想知道n是否为素数，还需不断选取a，看看上述同余式是否成立：如果有足够多的不同a能使同余式成立，则可以说n可能是素数，也可能是伪素数。之所以称为伪素数，是因为虽然选取的所有a能让等式成立，但事实上n却是一个合数。这时，同余式an-1≡1(mod n) 称为Fermat liar。
如果出现了任一个a使同余式不成立，即an-11(mod n)，则判定n是合数。a称为对于n的合数判定的Fermat witness。
基于上述结论，产生了测试素数的Miller-Rabin方法，这是一种随机化算法。设n为待检验的整数；k为选取a的次数。
重复k次计算，每次在［1, n-1］范围内随机选取一个a，若an -1 mod n ≠1 ，则返回n为合数；若随机选取的k个a都使an-1≡1(mod n)成立，则返回n为素数或伪素数的信息。
如果Miller-Rabin方法中的每次计算使用模指数运算的快速算法，则运行时间可提升至O(k×log32n)。下面，我们来证明这个方法的有效性。
证明：
如果a2≡1(mod n)，则必有a≡1(mod n)或a≡n-1(mod n) 。
现在，如果有一个大于2的质数n，令n-1=2s*d，其中d为奇数，则根据费马小定理，如果a不能被素数n整除，则 an-1≡1(mod n)。由上述性质可推出:
ad≡1(mod n) 或者a2≡1(mod n)，a2r*d≡n-1(mod n)，0≤r这样，我们只要找到一个a，满足ad1(mod n)或者存在一个自然数r(0≤r事实上经过反复检验，可以得出这样一个结论：
对于32位内的任一个整数n，如果其通过了以2、7、61为底的Miller-Rabin测试，那么n一定是素数；反之则一定不是。
下面，给出Miller-Rabin测试的算法模板。
int64 qpow(int64 a,int64 b,int64 M){ // 通过快速幂取模计算ab mod M
int64 ans=1;
while(b){
if(b&1) ans*=a,ans%=M;
a*=a;a%=M;b>>=1;
}
return ans;
}
bool MillerRabinTest(int64 x,int64 n){// 选取x为底，判定n是否可能为素数
int64 y=n-1;
while(!(y&1))y>>=1;// 略去n-1（=d*2s）右端连续的0，将其调整为d
x=qpow(x,y,n);// x=xd mod n
while(y<n-1&&x!=1&&x!=n-1)// 将x反复平方，直到模数值出现n-1或1为止
x=(x*x)%n,y<<=(int64)1;
return x==n-1||y&1==1;// 若x为n-1或y为奇数，则n可能是素数，否则一定不是
}
bool isprime32(int64 n){// 判断32位内的整数n是否为素数
if(n==2||n==7||n==61) return 1;// 若n为［2,7,61］中的元素，则n为素数
if(n==1||(n&1)==0) return 0;// 若n是1或是非2偶数，则n为合数
return MillerRabinTest(2,n)&&MillerRabinTest(7,n)&&MillerRabinTest(61,n);
// 若n通过以2、7、61为基准的Miller-Rabin测试，则n一定是素数；否则为合数
}
