我试图解决 SPOJ 上的问题。我们需要计算第 n 个孪生素数对(相差 2 的素数)。n 可以大到 10^5。我尝试使用筛子进行预先计算,我必须筛分到 10^8 才能获得最大的 n 个孪生素数,但是时间限制很严格(2 秒)并且超时。我注意到人们已经在 0.00 秒内解决了它,所以我在谷歌上四处寻找一个公式,但找不到任何有用的东西。有人可以指导我吗?
提前致谢!!
出于好奇,我使用 Eratosthenes 筛的两种变体解决了这个问题。第一个变体在 0.93 秒内在测试机上完成,第二个在 0.24 秒内完成。作为比较,在我的电脑上,第一个在 0.08 秒内完成,第二个在 0.04 秒内完成。
第一个是对奇数的标准筛子,第二个是一个稍微复杂的筛子,除了偶数之外还省略了 3 的倍数。
SPOJ 的测试机器又旧又慢,所以程序在它们上运行的时间比在典型的新机器上运行的时间要长得多;而且它们的缓存很小,因此保持计算量很小很重要。
这样做,Eratosthenes 的筛子就足够快了。但是,保持较小的内存使用量非常重要。第一个变体,每个数字使用一个字节,在 SPOJ 上给出了“超出时间限制”,但在我的盒子上运行了 0.12 秒。因此,鉴于 SPOJ 试验机的特性,使用位筛在给定时间内解决它。
在 SPOJ 机器上,通过进一步将筛子的空间减少一半,我获得了显着的加速(运行时间 0.14 秒)。因为 - 除了第一对 (3,5) - 所有素数双胞胎都有形式,如果不产生双素数对,(6*k-1, 6*k+1)
你不需要知道这两个数字中的哪个是合数,只筛选就足够了k
指数k
。
(6*k + 1
可被 5 整除当且仅当k = 5*m + 4
对于某些m
,并且6*k - 1
可被 5 整除当且仅当k = 5*m+1
对于某些m
,所以 5 将标记5*m ± 1, m >= 1
为不产生孪生素数。类似地,6*k+1
可被 13 整除当且仅当k = 13*m + 2
对于某些m
且6*k - 1
当且仅当k = 13*m - 2
对于 some m
,所以 13 会标记出来13*m ± 2
。)
这不会改变标记的数量,因此在缓存足够大的情况下,运行时间的变化很小,但对于小缓存来说,这是一个显着的加速。
不过,还有一件事。你10 8的限制太高了。我使用了一个下限(2000 万),它不会高估第 100,000个孪生素数对。限制为 10 8,第一个变体肯定不会及时完成,第二个可能不会。
随着限制的减少,需要对阿特金筛子进行一些优化以击败省略偶数和 3 倍数的 Eratosthenes 变体,一个简单的实现将明显变慢。
关于您的(维基百科的伪代码)阿特金筛子的一些评论:
#define limit 100000000
int prime1[MAXN];
int prime2[MAXN];
您不需要第二个数组,素数双胞胎中较大的配对可以很容易地从较小的配对中计算出来。您正在浪费空间并破坏从两个数组读取的缓存位置。(不过,与筛分所需的时间相比,这微不足道。)
int root = ceil(sqrt(limit));
bool sieve[limit];
在当今的许多操作系统上,即使限制有所减少,这也是即时的段错误。堆栈大小通常限制为 8MB 或更小。应该在堆上分配该大小的数组。
如上所述,bool
每个数字使用一个会使程序运行得比必要的慢得多。您应该自己使用std::bitset
或std::vector<bool>
或旋转这些位。此外,建议至少省略偶数。
for (int x = 1; x <= root; x++)
{
for (int y = 1; y <= root; y++)
{
//Main part of Sieve of Atkin
int n = (4*x*x)+(y*y);
if (n <= limit && (n % 12 == 1 || n % 12 == 5)) sieve[n] ^= true;
n = (3*x*x)+(y*y);
if (n <= limit && n % 12 == 7) sieve[n] ^= true;
n = (3*x*x)-(y*y);
if (x > y && n <= limit && n % 12 == 11) sieve[n] ^= true;
}
}
这是非常低效的。它尝试了太多的 xy 组合,对于每个组合,它会进行三到四次除法以检查余数模 12,并且它在数组中来回跳跃。
分离不同的二次方。
因为4*x^2 + y^2
,很明显,你只需要考虑x < sqrt(limit)/2
和奇数y
。那么余数模 12 是 1、5 或 9。如果余数是 9,则4*x^2 + y^2
实际上是 9 的倍数,所以这样的数将被排除为不是无平方数。然而,最好从筛子中完全省略 3 的倍数,n % 12 == 1
并n % 12 == 5
分别处理这些情况。
对于3*x^2 + y^2
,很明显,您只需要考虑x < sqrt(limit/3)
并且稍微思考一下就会发现它x
必须是奇数和y
偶数(并且不能被 3 整除)。
对于3*x^2 - y^2
with y < x
,显然您只需要考虑y < sqrt(limit/2)
。查看余数模 12,您会看到它y
不能被 3 整除,x
并且y
必须具有不同的奇偶性。
我在 0.66 秒内获得了交流电。因为,有 0.0s 的解决方案,我认为可以进行更好的优化,但是,我在这里描述了我的方法。
我在Sieve of Eratosthenes
. 您知道这2
是唯一的偶数素数,使用它可以将计算素数的计算时间和内存减少一半。
其次,所有作为孪生素数的数字都不会是2
和的倍数3
(因为它们是素数!)。所以,这些数字将是形式6N+1
和6N+5
(其余的肯定不是素数)。6N+5 = 6N+6-1 = 6(N+1)-1
. 所以可以看出6N+1
并且6N-1
可能是N
> = 1的孪生素数。因此,您使用之前计算的素数预先计算所有这些值。(小例是 3 5)
注意:在 10^8 之前不需要计算素数,上限要低得多。[编辑:如果你愿意,我可以分享我的代码,但如果你自己想出一个解决方案会更好。:)]
因此,根据 Wolfram Alpha 的说法,基本上,筛选多达 20,000,000 个就足够了。vector<bool>
在 C++ 中使用 Eratosthenes 的普通筛子(顺便说一句,您使用的是什么语言?)。
跟踪筛环内的孪生素数。当你找到双胞胎时,将一对的下质数存储在一个单独的向量中,如果请求一个无序(小于前一个)索引(它们是,与描述页面上显示的示例相反),只需从此存储中获取素数:
size_t n = 10000000, itop=2236;
vector<bool> s;
vector<int> twins;
s.resize(n, true);
int cnt, k1, k2, p1=3, p2, k=0;
cin >> cnt;
if( cnt-- > 0 )
{
cin >> k1;
for( size_t i=1; i < n; ++i ) // p=2i+1
{
if( s[i] )
{
p2 = 2*i+1;
if( p2-p1 == 2 ) { ++k; twins.push_back(p1); }
if( k==k1 )
{
cout << p1 << " " << p2 << endl;
......
等 1.05 秒(在 Ideone 上为 0.18 秒)接受。或者解开逻辑 - 只需立即预先计算 100,000 个孪生素数对,然后在单独的循环中访问它们(0.94 秒)。
可以在此处找到解决此问题的有效算法的描述@Programming Praxis entry此外,还提供了 Scheme 和 Perl 示例代码。
我使用埃拉托色尼筛法预先计算了一个大的素数列表,然后遍历列表,计算比它们的后继少 2 个的项目,直到找到其中的 n 个。在http://ideone.com/vYjuC上运行 1.42 秒。我也想知道如何在零秒内计算答案。
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define ISBITSET(x, i) (( x[i>>3] & (1<<(i&7)) ) != 0)
#define SETBIT(x, i) x[i>>3] |= (1<<(i&7));
#define CLEARBIT(x, i) x[i>>3] &= (1<<(i&7)) ^ 0xFF;
typedef struct list {
int data;
struct list *next;
} List;
List *insert(int data, List *next)
{
List *new;
new = malloc(sizeof(List));
new->data = data;
new->next = next;
return new;
}
List *reverse(List *list) {
List *new = NULL;
List *next;
while (list != NULL)
{
next = list->next;
list->next = new;
new = list;
list = next;
}
return new;
}
int length(List *xs)
{
int len = 0;
while (xs != NULL)
{
len += 1;
xs = xs->next;
}
return len;
}
List *primes(int n)
{
int m = (n-1) / 2;
char b[m/8+1];
int i = 0;
int p = 3;
List *ps = NULL;
int j;
ps = insert(2, ps);
memset(b, 255, sizeof(b));
while (p*p < n)
{
if (ISBITSET(b,i))
{
ps = insert(p, ps);
j = (p*p - 3) / 2;
while (j < m)
{
CLEARBIT(b, j);
j += p;
}
}
i += 1; p += 2;
}
while (i < m)
{
if (ISBITSET(b,i))
{
ps = insert(p, ps);
}
i += 1; p += 2;
}
return reverse(ps);
}
int nth_twin(int n, List *ps)
{
while (ps->next != NULL)
{
if (n == 0)
{
return ps->data - 1;
}
if (ps->next->data - ps->data == 2)
{
--n;
}
ps = ps->next;
}
return 0;
}
int main(int argc, char *argv[])
{
List *ps = primes(100000000);
printf("%d\n", nth_twin(100000, ps));
return 0;
}
这就是我所尝试的。我有一串 TLE。
bool mark [N];
vector <int> primeList;
void sieve ()
{
memset (mark, true, sizeof (mark));
mark [0] = mark [1] = false;
for ( int i = 4; i < N; i += 2 )
mark [i] = false;
for ( int i = 3; i * i <= N; i++ )
{
if ( mark [i] )
{
for ( int j = i * i; j < N; j += 2 * i )
mark [j] = false;
}
}
primeList.clear ();
primeList.push_back (2);
for ( int i = 3; i < N; i += 2 )
{
if ( mark [i] )
primeList.push_back (i);
}
//printf ("%d\n", primeList.size ());
}
int main ()
{
sieve ();
vector <int> twinPrime;
for ( size_t i = 1; i < primeList.size (); i++ )
{
if ( primeList [i] - primeList [i - 1] == 2 )
twinPrime.push_back (primeList [i - 1]);
}
int t;
scanf("%d",&t);
int s;
while ( t-- )
{
scanf("%d",&s);
printf ("%d %d\n", twinPrime [s - 1], twinPrime [s - 1] + 2);
}
return 0;
}
这是一个可以回答您的问题的程序:
素数除以 3 时,修正为十进制 0(零)后具有相等的商数是孪生素数。
这可以写成
对于任何一对素数 Px, Py, 如果 [Px/3, 0] = [Py/3, 0] 那么 Px 和 Py 是素数孪晶。
其基础是,如果素数相差 2,则当商被修正为十进制零时,除以所有感兴趣的素数将产生唯一的等商。未被 2 分隔的质数在修正为十进制零时将不具有等商。
例如:
• 11, 13 除以 3 将产生唯一的唯一商,当商被修正为十进制零时。
• 17、19 除以 3 时,商被修正为十进制零时,将得出唯一的商 6。
• 29, 31 除以 3 时,将得到唯一商 10,此时商被修正为十进制零。
等等。
下面是一个使用 Excel 的简单过程:
• 从任何素数列表中找出素数孪生 • 在任何素数范围内找出孪生素数 • 找出最大的素数孪生素数 • 找出孪生素数之间的间隙
要找到最大的孪生素数,请使用上述程序,将已知最大素数的范围放入第 1 列(例如,最高的 10k 素数)。
如果在此范围内未找到素数孪生,则转到下一个最低范围,直到找到孪生素数。这将是最大的孪生素数。
希望这可以帮助。