查找最近的较小素数的快速算法

如果输入的大小是有限的，在预先计算的素数表中查找可能是最快的。

这是Daniel Fischer在评论中提到的Baillie-Wagstaff伪最优性测试的伪代码实现。我们从一个简单的埃拉托色尼筛子开始，我们以后会用到它。

function primes(n)
ps := []
sieve := makeArray(2..n, True)
for p from 2 to n step 1
if sieve(p)
ps.append(p)
for i from p * p to n step p
sieve[i] := False
return ps

powerMod函数将基底b提升为指数e，所有计算都以m模进行;它比先求幂，然后求结果的模要快得多，因为中间的计算量会很大。

function powerMod(b, e, m)
x := 1
while e > 0
if e % 2 == 1
x := (b * x) % m
b := (b * b) % m
e := floor(e / 2)
return x

数论中的jacobi函数告诉我们a是否对p取二次余数。

function jacobi(a, p)
a := a % p
t := 1
while a != 0
while a % 2 == 0
a := a / 2
if p % 8 == 3 or p % 8 == 5
t := -t
a, p := p , a # swap
if a % 4 == 3 and p % 4 == 3
t := -t
a := a % p
if p == 1 return t else return 0

Gary Miller的强伪素数检验是基于Pierre de Fermat的小定理，该定理表明如果p是素数，那么对于任何a!= 0，a^ (p- 1) == 1 (modp)。米勒的检验比费马的强一些，因为它不会被卡迈克尔数愚弄。

function isStrongPseudoprime(n, a)
d := n - 1; s := 0
while d % 2 == 0
d := d / 2; s := s + 1
t = powerMod(a, d, n)
if t == 1 return ProbablyPrime
while s > 0
if t == n - 1 return ProbablyPrime
t := (t * t) % n; s := s - 1
return Composite

Miller-Rabin检验执行k强伪素数检验，其中k通常在10到25之间。强伪素数测试可以被愚弄，但如果你执行足够多的测试，被愚弄的可能性非常小。

function isPrime(n) # Miller-Rabin
for i from 1 to k
a := randInt(2 .. n-1)
if not isStrongPseudoprime(n, a)
return Composite
return ProbablyPrime

这个素数检验在大多数情况下是足够的，而且速度也足够快。但如果你想要更强一点、更快一点的东西，可以使用基于卢卡斯链的测试。这是卢卡斯链的计算。

function chain(n, u, v, u2, v2, d, q, m)
k := q
while m > 0
u2 := (u2 * v2) % n; v2 := (v2 * v2 - 2 * q) % n
q := (q * q) % n
if m % 2 == 1
t1 := u2 * v; t2 := u * v2
t3 := v2 * v; t4 := u2 * u * d
u, v := t1 + t2, t3 + t4
if u % 2 == 1 u := u + n
if v % 2 == 1 v := v + n
u, v, k := (u / 2) % n, (v / 2) % n), (q * k) % n
m := floor(m / 2)
return u, v, k

通常使用John Selfridge提出的算法初始化Lucas链。

function selfridge(n)
d, s := 5, 1; ds := d * s
repeat
if gcd(ds, n) > 1 return ds, 0, 0
if jacobi(ds, n) == 1 return ds, 1, (1 - ds) / 4
d, s := d + 2, s * -1; ds := d * s

则卢卡斯伪素数检验确定一个数是素数还是可能是合数。就像费马测试一样，它有两种风格，一种是标准的，一种是强的，和费马测试一样，它也可以被愚弄，尽管费马测试的错误是合数可能被错误地报告为素数，但卢卡斯测试的错误是素数可能被错误地报告为合数。

function isLucasPseudoprime(n) # standard
d, p, q := selfridge(n)
if p == 0 return n == d
u, v, k := chain(n, 0, 2, 1, p, d, q, (n + 1) / 2)
return u == 0
function isLucasPseudoprime(n) # strong
d, p, q := selfridge(n)
if p == 0 return n == d
s, t := 0, n + 1
while t % 2 == 0
s, t := s + 1, t / 2
u, v, k := chain(n, 1, p, 1, p, d, q, t // 2
if u == 0 or v == 0 return Prime
r := 1
while r < s
v := (v * v - 2 * k) % n; k := (K * k) % n
if v == 0 return Prime
return ProbablyComposite

那么Baillie-Wagstaff检验就很简单了。首先检查输入是否小于2或完全平方(检查平方根是否为整数)。然后试除小于100的素数，很快找到大多数合数，最后以2为基数进行强伪素数检验(有些人为了更确定，还以3为基数进行强伪素数检验)，然后进行卢卡斯伪素数检验，最终确定合数。

function isPrime(n) # Baillie-Wagstaff
if n < 2 or isSquare(n) return False
for p in primes(100)
if n % p == 0 return n == p
return isStrongPseudoprime(n, 2) 
and isLucasPseudoprime(n) # standard or strong

Baillie-Wagstaff试验没有已知错误。

一旦你有一个好的素数检验，你可以找到小于n的最大素数，从n开始倒数，直到第一个素数。

如果你对用素数编程感兴趣，我在我的博客上谦虚地推荐这篇文章，或者许多其他与素数相关的博客文章，你可以通过使用博客上的搜索功能找到它们。

除了上面还注意到Bertrand的假设，总存在至少一个素数p，其中n<p<2n-2。这就给出了上界

探讨Miller-Rabin素数检验。这是概率性的，但如果你做几百次，这几乎可以保证在long long范围内的精度。

同样，如果你会使用Java,BigInteger.isProbablePrime也可以提供帮助。C c++似乎没有一个内置的函数来测试素数。

你的任务很有趣!

我决定用纯c++从头开始为您实现一个相当先进，庞大但快速(高效)的解决方案(编译时需要c++ 20标准)。

在我的代码中使用了以下补充算法-埃拉托色尼筛，费马素数检验，米勒拉宾素数检验，试除法(使用素数)，二分查找。

还使用外部c++库Boost Multiprecision来实现大整数运算。你可以使用其他库，如果你想，我的代码是这样一种方式，任何库都可以使用。CLang/GCC编译器也有__int128类型，如果你的任务只在128位范围内，你可以用它来代替大整数。

我的代码支持任何位大小的整数，但我彻底测试了正确性和速度，直到1024位(即测试了8,16,32,64,128,256,512,1024位)。

所有可能的计算我都使用多核，如果你有8个核(实际上是硬件线程)，那么在所有的计算中，我启动一个8个线程池并使用它们。但是这种多线程方法只用于大于128位的整数，这是一种默认启用的启发式算法，但如果需要，您可以通过在主函数PrevPrime()中设置额外的参数来强制任何位大小的单线程或多线程算法。

我的主要算法(PrevPrime()函数)执行以下步骤:

1. 有一次(整个程序运行一次)，我使用埃拉托色尼筛子计算了一个小于2^20的素数表。该表最多可容纳1024位整数。如果您使用2048位或更大的质数可以使用。

2. 当我们给定N时，我开始从N向下，再次使用埃拉托色尼筛子从范围(N - length, N)中过滤出合数。不是所有的2^20素数都被使用，实验上我计算了对于每个输入N位大小的素数的最佳数量(对于32位100个素数，对于64位500个素数，对于128位- 1200个素数，对于256位- 7000个素数，对于512位- 17500个素数，对于1024位- 75000个素数)。

3. Eratosthenes筛需要开始筛选每个特定素数p从一些偏移N - off，其中off = N % p。因为除法是CPU中最昂贵的基本运算，所以我使用多线程(多核)方法来计算偏移量(对于较大的整数，对于较小的单线程使用)。

4. 筛选完成后，我们可以确定范围(N - length, N)中哪些数字是合数，但我们不确定其余的数字是质数还是合数。因此，最后一步是对每个剩余的数字运行米勒-拉宾的素数测试。这些测试也使用多线程方法完成(仅针对较大的整数)。

PrevPrime()的计时遵循我的旧的和缓慢的1.2 GHzCPU，它有8个硬件线程，不同的整数位大小的计时:

32-bit -小于0.00005 sec/1 PrevPrime()整数，64-bit -0.0005 sec，128-bit -0.002秒，256-bit -0.012秒，512-bit -0.07-0.1秒，1024-bit -0.55-1.3秒。

所以你可以看到在我的慢1.2 GHz CPU 1024位整数在1秒内输出前一个素数。

上面的计时不仅取决于整数的位大小，还取决于素数的密度。众所周知，对于位大小为B的数，质数在ln(2^B)中平均出现一次。因此，更大的数不仅具有更大整数的数学难度，而且具有更多的稀有素数，因此更大的数在输出前一个素数之前需要检查更多的数。

为了了解在PrevPrime()中筛选的效率，我提供了Fermat和Miller-Rabin测试本身的计时(F是Fermat,MR是Miller-Rabin):

128-位-复合MR &F -0.0002秒，素数MR -0.0026秒F -0.0054秒;256年-位-复合MR &F -0.0007秒，质数MR -0.0116秒F -0.0234秒;512年-位-复合MR &F -0.0043秒，素数MR -0.0773秒F -0.1558秒;1024年-位-复合MR &F -0.0290秒，MR -0.4518秒F -0.9085秒

从上面的统计数据我们可以得出结论，当目标是相同的失败概率时，费马检验比米勒-拉宾检验花费的时间多两倍。

正如你所看到的，单个米勒-拉宾对1024位素数的测试大约是0.5秒，而整个PrevPrime()在1秒内运行，这意味着PrevPrime()是如此高效，它只比使用米勒-拉宾对素数的单次测试慢两倍。

当然，你可以注意到我使用了概率素数测试，因为1024位数字确定性测试需要非常非常大的时间。

对于概率测试，我选择失败的目标概率等于1/2^32 = 1/4294967296，因此素数测试失败的情况非常非常罕见。你可以选择较小的概率，如果你需要更高的精度，它在我的程序中是可以调整的。

从理论上讲，费马的每个单次测试都有1/2的失败几率(不包括几乎100%失败几率的卡迈克尔数)，而米勒-拉宾的失败几率为1/4。因此，要达到目标概率1/2^32，需要费马32步，米勒-拉宾16步。

每个位大小所需质数的阈值可以自动调整，但需要大量的单次预计算，这就是为什么我计算表并硬编码它，但是计算表的代码仍然存在于我的代码中。如果您需要重新计算阈值，然后转到函数ComputePrevPrimeOptimalPrimesCnt()并注释掉第一行(使用return bits <= 32 ? 128 : .....)并运行程序，它将输出适合您的PC和位大小的最佳大小。但是我的预计算表应该适用于任何CPU，只有新的位大小需要调整。

如果你不知道我唯一依赖的库是Boost。它可以通过sudo apt install libboost-dev-all在Linux下轻松安装。在Windows下稍微困难一点，但仍然可以通过Chocolatey的Boost Package安装(执行choco install boost-msvc-14.3，但先安装Chocolatey)。

默认情况下，当程序运行时，它运行所有测试，测试是代码中以宏TEST(TestName) { .....开头的所有函数。如果需要额外的测试，只需编写像TEST(NewTest) { ...code here... }这样的主体，不需要额外的测试注册。如果您确信在使用我的库时也可以删除所有TEST()函数，则它们不需要正常工作，除非它们对完全测试我的库很有用。

重要通知! !． 这段代码是我在一天内写的，所以它没有经过很好的测试，可能在某个地方包含一些bug。因此，不建议在没有额外注意查看代码和/或额外测试的情况下直接在生产中使用它。它更像是教育代码，可以作为你编写自己的生产就绪代码的一种指导。

代码如下。

上网试试!

(注意，Try it online!在线运行限制为128位(1024可能)，由于GodBolt服务器限制程序总运行时间为3秒，下载后将test_till_bits = 128更改为1024)。

源代码在这里． 托管在Github Gist源代码。不幸的是，由于StackOverflow的30K符号限制，我不能在这里内联完整的源代码，因为只有我的源代码有30K字节，不包括上面的长英文描述(这是额外的10K)。因此，我在Github GIST(上面的链接)和GodBolt服务器(上面的Try it online!链接)上分享源代码，两个链接都有完整的源代码。

代码控制台输出:

'Test GenPrimes' 0.006 sec
'Test GenPrimes_CntPrimes' 0.047 sec
'Test MillerRabin' 0.006 sec
Fermat_vs_MillerRabin: bits    8 prp_cnt  58 mr_time 0.0000|0.0000 f_time 0.0000|0.0000, total time 0.009 sec
Fermat_vs_MillerRabin: bits   16 prp_cnt  57 mr_time 0.0000|0.0000 f_time 0.0000|0.0000, total time 0.001 sec
Fermat_vs_MillerRabin: bits   32 prp_cnt  46 mr_time 0.0000|0.0000 f_time 0.0000|0.0000, total time 0.004 sec
Fermat_vs_MillerRabin: bits   64 prp_cnt  22 mr_time 0.0000|0.0005 f_time 0.0000|0.0011, total time 0.041 sec
Fermat_vs_MillerRabin: bits  128 prp_cnt  13 mr_time 0.0002|0.0028 f_time 0.0002|0.0056, total time 0.134 sec
Fermat_vs_MillerRabin: bits  256 prp_cnt   6 mr_time 0.0008|0.0120 f_time 0.0008|0.0245, total time 0.321 sec
Fermat_vs_MillerRabin: bits  512 prp_cnt   1 mr_time 0.0042|0.0618 f_time 0.0042|0.1242, total time 0.857 sec
Fermat_vs_MillerRabin: bits 1024 prp_cnt   1 mr_time 0.0273|0.4101 f_time 0.0273|0.8232, total time 3.807 sec
'Test IsProbablyPrime_Fermat_vs_MillerRabin' 5.176 sec
'Test PrevPrime_ReCheckWith_PrimesDiv_and_Fermat' 5.425 sec
PrevPrime    8-bit threads:1  0.0000 sec, avg distance    3.3, total time 0.001 sec
PrevPrime    8-bit threads:8  0.0000 sec, avg distance    3.3, total time 0.000 sec
PrevPrime   16-bit threads:1  0.0000 sec, avg distance   11.2, total time 0.000 sec
PrevPrime   16-bit threads:8  0.0000 sec, avg distance   11.2, total time 0.000 sec
PrevPrime   32-bit threads:1  0.0001 sec, avg distance   10.5, total time 0.001 sec
PrevPrime   32-bit threads:8  0.0001 sec, avg distance   10.5, total time 0.001 sec
PrevPrime   64-bit threads:1  0.0014 sec, avg distance   39.7, total time 0.025 sec
PrevPrime   64-bit threads:8  0.0012 sec, avg distance   39.7, total time 0.023 sec
PrevPrime  128-bit threads:1  0.0110 sec, avg distance   85.2, total time 0.170 sec
PrevPrime  128-bit threads:8  0.0084 sec, avg distance   85.2, total time 0.142 sec
PrevPrime  256-bit threads:1  0.0452 sec, avg distance  207.4, total time 0.570 sec
PrevPrime  256-bit threads:8  0.0331 sec, avg distance  207.4, total time 0.473 sec
PrevPrime  512-bit threads:1  0.1748 sec, avg distance  154.0, total time 2.246 sec
PrevPrime  512-bit threads:8  0.1429 sec, avg distance  154.0, total time 2.027 sec
PrevPrime 1024-bit threads:1  2.1249 sec, avg distance  379.4, total time 15.401 sec
PrevPrime 1024-bit threads:8  1.6030 sec, avg distance  379.4, total time 12.778 sec
'Test PrevPrime' 33.862 sec
'Program run time' 44.524 sec

看起来你正在解决这个问题。

正如@Ziyao Wei所说，你可以简单地使用Miller-Rabin素数测试来解决它。

这是我的解决方案

#include<cstdio>
#include<cstdlib>
#include<ctime>
short T;
unsigned long long n;
inline unsigned long long multi_mod(const unsigned long long &a,unsigned long long b,const unsigned long long &n)
{
unsigned long long exp(a%n),tmp(0);
while(b)
{
if(b&1)
{
tmp+=exp;
if(tmp>n)
tmp-=n;
}
exp<<=1;
if(exp>n)
exp-=n;
b>>=1;
}
return tmp;
}
inline unsigned long long exp_mod(unsigned long long a,unsigned long long b,const unsigned long long &c)
{
unsigned long long tmp(1);
while(b)
{
if(b&1)
tmp=multi_mod(tmp,a,c);
a=multi_mod(a,a,c);
b>>=1;
}
return tmp;
}
inline bool miller_rabbin(const unsigned long long &n,short T)
{
if(n==2)
return true;
if(n<2 || !(n&1))
return false;
unsigned long long a,u(n-1),x,y;
short t(0),i;
while(!(u&1))
{
++t;
u>>=1;
}
while(T--)
{
a=rand()%(n-1)+1;
x=exp_mod(a,u,n);
for(i=0;i<t;++i)
{
y=multi_mod(x,x,n);
if(y==1 && x!=1 && x!=n-1)
return false;
x=y;
}
if(y!=1)
return false;
}
return true;
}
int main()
{
srand(time(NULL));
scanf("%hd",&T);
while(T--)
{
for(scanf("%llu",&n);!miller_rabbin(n,20);--n);
printf("%llun",n);
}
return 0;
}