10^67以上10^12个整数的相邻素数间隔暴力搜索

无心人 · 发表于 2019-9-11 17:31:55

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？欢迎注册

×

实际从10^67-10开始搜索
因为10^67-10能分解成2*3^3*5*7*11^2*13*23*37*67*4093*8779*21649
30030=2*3*5*7*11*13
然后，对所有候选，筛掉10000以内素因子
再miller-rabin测试1次
每次搜索30030个整数，每线程进行333000次暴力搜索，保存大于等于2000的差的相邻结果

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
#include <time.h>
#define step (2*3*5*7*11*13)
#define slice 333000
#define i64 long long
#define u64 unsigned long long
#define i32 long
#define u32 unsigned long
char temp[step];
char work[step];
char n4[10000];
u32 prime[1300];
u32 primeCount = 0;
i64 prev, gap, cnt;
i64 base12, offset, curr;
mpz_t base67, start, rem, test;
void prepareTemp( void )
{
u32 smallPrime[] = {2, 3, 5, 7, 11, 13};
for (i32 i = 0; i < step; i += 2) temp[i] = 0;
for (i32 i = 1; i < step; i += 2) temp[i] = 1;
for (i32 j = 1; j < 6; j ++)
for (i32 i = smallPrime[j]; i < step; i += 2* smallPrime[j])
temp[i] = 0;
}
void copyToWork( void )
{
memcpy(work, temp, step);
}
void preparePrime( void )
{
char n2[100];
for (i32 i = 0; i < 100; i+=2) n2[i] = 0;
for (i32 i = 3; i < 100; i+=2) n2[i] = 1;
n2[2] = 1; n2[9] = 0;
for (i32 i = 0; i < 10; i ++) if (n2[i] == 1) prime[primeCount++] = i;
for (i32 i = 1; i < primeCount; i ++)
{
i32 p = prime[i];
i32 t = p * p;
while (t < 100)
{
n2[t] = 0;
t += 2 * p;
}
}
for (i32 i = 11; i < 100; i += 2) if (n2[i] == 1) prime[primeCount++] = i;
for (i32 i = 4; i < 10000; i+=2) n4[i] = 0;
for (i32 i = 3; i < 10000; i+=2) n4[i] = 1;
for (i32 i = 1; i < primeCount; i ++)
{
i32 p = prime[i];
i32 t = p * p;
while (t < 10000)
{
n4[t] = 0;
t += 2 * p;
}
}
for (i32 i = 101; i < 10000; i += 2) if (n4[i] == 1) prime[primeCount++] = i;
}
void init( u32 index )
{
//base67 = 10^67 - 10
mpz_init(base67);
mpz_ui_pow_ui(base67, 10, 67);
mpz_sub_ui(base67, base67, 10);
mpz_init2(start, 256);
mpz_init2(test, 256);
mpz_init2(rem, 256);
preparePrime();
prepareTemp();
base12 = index * slice * step;
mpz_add_ui(base67, base67, base12);
mpz_sub_ui(test, base67, 1);
for (i32 i = 1; i < step; i +=2)
{
if (mpz_probab_prime_p(test, 10) >= 1)
{
mpz_sub(test, test, base67);
prev = mpz_get_si(test);
break;
}
mpz_sub_ui(test, test, 2);
}
cnt = 0;
gap = 0;
printf("from 10^67 - 10 + %lld ...\n", base12);
}
void check( i64 i )
{
{
curr = i + offset;
gap = curr - prev;
if ( gap >= 2000 )
{
//10^67 - 10 + base12 + ...
printf("%lld, %lld, %lld\n", prev, curr, gap);
}
prev = curr;
cnt ++;
}
}
void oneStep( mpz_t start )
{
copyToWork();
for (i32 i = 6; i < primeCount; i ++)
{
u32 r = mpz_mod_ui(rem, start, prime[i]);
r = prime[i] - r;
if (r%2==0) r += prime[i];
while (r < step)
{
work[r] = 0;
r += 2 * prime[i];
}
}
for (i32 i = 1; i < step; i += 2)
if (work[i] == 1)
{
mpz_add_ui(test, start, i);
if (mpz_millerrabin(test, 1) == 0)
work[i] = 0;
}
/*
int b = 0;
printf("----------------------------\nCheck work..........\n");
for (i32 i = 1; i < step; i += 2)
{
if (work[i] == 1)
{
printf("(%d, %d)", b, i);
b ++;
if (b >= 10) break;
}
}
printf("\n--------------------------\n");
*/
for (i32 i = 1; i < step; i += 2)
if (work[i] == 1)
check( i );
}
void circle( u32 index )
{
for (u32 i = 0; i < slice; i ++)
{
offset = i * step;
mpz_add_ui(start, base67, offset);
//printf("i = %d, offset = %lld\n", i, offset);
//gmp_printf("start = %Zd\n", start);
oneStep( start );
}
}
i32 main( int argc, char * argv[] )
{
u32 index = 0;
time_t timer;
char buffer[80];
struct tm * timeinfo;
if (argc > 1)
index = strtol(argv[1], NULL, 10);
else
exit(1);
init(index);
//gmp_printf("start from %Zd\n", base67);
time(&timer);
timeinfo = localtime (&timer);
strftime (buffer,80,"start time: %D %T.\n", timeinfo);
puts(buffer);
circle(index);
time(&timer);
timeinfo = localtime (&timer);
strftime (buffer,80,"\nend time: %D %T.\n", timeinfo);
puts(buffer);
//mpz_add_ui(test, base67, curr);
//gmp_printf("end with %Zd\n", test);
return 1;
}

复制代码

无心人 · 发表于 2019-9-12 17:36:45

最后结果里相邻素数差大于等于3000的

base67 = 10^67 - 10
base12 = index * 333000 * 30030
base67 + base12 + ...

index = 10
9905037943, 9905040953, 3010
index = 20
8431176409, 8431179511, 3102
index = 39
3210833959, 3210837001, 3042
index = 4
7309153469, 7309156469, 3000
index = 41
3330732749, 3330735757, 3008
index = 53
6540340721, 6540343861, 3140
index = 64
2997043067, 2997046211, 3144
index = 7
8918617231, 8918620237, 3006
index = 78
8040089569, 8040092633, 3064
index = 87
3378529073, 3378532187, 3114
index = 90
3920277503, 3920280509, 3006
index = 94
1470578173, 1470581179, 3006
index = 99
9648273853, 9648276893, 3040

wayne · 发表于 2019-9-12 21:43:02

#define i32 long
#define u32 unsigned long

复制代码

现在基本上都是64位的天下, 很少有人玩32位机器吧, 于是64位机器下,是不是命名有误导性质了,

无心人 · 发表于 2019-9-15 09:35:05

wayne 发表于 2019-9-12 21:43
现在基本上都是64位的天下, 很少有人玩32位机器吧, 于是64位机器下,是不是命名有误导性质了,

要节约资源啊，能用32还是用32吧，虽然还是占64位寄存器，这叫绿色环保意识

熊一兵广义概率 · 发表于 2020-2-10 11:12:56

在《概率素数论》中有一个x处相邻素数最大间距d(x)定理：
d(x)=(lnx)^2，请计算实际值验证，

账号		自动登录	找回密码
密码			欢迎注册

[分享] 10^67以上10^12个整数的相邻素数间隔暴力搜索

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。