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;
char work;
char n4;
u32 prime;
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 = 0;
for (i32 i = 1; i < step; i += 2) temp = 1;
for (i32 j = 1; j < 6; j ++)
for (i32 i = smallPrime; i < step; i += 2* smallPrime)
temp = 0;
}

void copyToWork( void )
{
memcpy(work, temp, step);
}

void preparePrime( void )
{
char n2;
for (i32 i = 0; i < 100; i+=2)n2 = 0;
for (i32 i = 3; i < 100; i+=2)n2 = 1;
n2 = 1; n2 = 0;
for (i32 i = 0; i < 10; i ++)if (n2 == 1) prime = i;
for (i32 i = 1; i < primeCount; i ++)
{
i32 p = prime;
   i32 t = p * p;
   while (t < 100)
{
      n2 = 0;
      t += 2 * p;
}
}

for (i32 i = 11; i < 100; i += 2)if (n2 == 1) prime = i;

for (i32 i = 4; i < 10000; i+=2)n4 = 0;
for (i32 i = 3; i < 10000; i+=2)n4 = 1;
for (i32 i = 1; i < primeCount; i ++)
{
i32 p = prime;
   i32 t = p * p;
   while (t < 10000)
{
n4 = 0;
t += 2 * p;
}
}
for (i32 i = 101; i < 10000; i += 2)if (n4 == 1) prime = 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);
r = prime - r;
if (r%2==0) r += prime;
while (r < step)
{
work = 0;
r += 2 * prime;
}
}

for (i32 i = 1; i < step; i += 2)
if (work == 1)
{
mpz_add_ui(test, start, i);
if (mpz_millerrabin(test, 1) == 0)
work = 0;
}
/*
int b = 0;
printf("----------------------------\nCheck work..........\n");
for (i32 i = 1; i < step; i += 2)
{

if (work == 1)
{
printf("(%d, %d)", b, i);
b ++;
if (b >= 10) break;
}
}
printf("\n--------------------------\n");
*/
for (i32 i = 1; i < step; i += 2)
if (work == 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;
struct tm * timeinfo;

if (argc > 1)
index = strtol(argv, 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位机器下,是不是命名有误导性质了, :lol

无心人 发表于 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，请计算实际值验证，

页: [1]

数学研发论坛's Archiver

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