估计下找到10^100以上的一个四生素数的工作量

gxqcn · 发表于 2008-7-16 12:38:07

原帖由 无心人 于 2008-7-16 10:28 发表 GxQ尝试过多大的数字的素性测试具体时间是多少？

没具体测试 HugeCalc 素性判定的能力，应该说是足够大的，可以容纳当前已知的最大素数。其算法还可进一步优化，留待下一版吧。。。

无心人 · 发表于 2008-7-16 13:04:33

我就想知道进行一个10000位数字的一次测试需要多少秒？是一次测试，即单独一个素数的测试

gxqcn · 发表于 2008-7-16 13:25:30

HugeCalc 内部针对不同的数据特点有不同的加速算法，比如判定 $10^10000 +- \delta \ ( \delta\ "is a small integer")$ 将比普通的 10000 位整数快很多。

无心人 · 发表于 2008-7-16 14:44:49

那就按判定普通数字测试下

无心人 · 发表于 2008-7-16 15:46:07

我评估了下素性测试的时间因为是linux系统，用的GMP Input:67 Begin Test 2^67 - 1. Composite Number. Used Time:0.000097 Input:31 Begin Test 2^31 - 1. Probably prime Number. Used Time:0.000086 Input:127 Begin Test 2^127 - 1. Probably prime Number. Used Time:0.000119 Input:61 Begin Test 2^61 - 1. Probably prime Number. Used Time:0.000070 Input:257 Begin Test 2^257 - 1. Composite Number. Used Time:0.000249 Input:4423 Begin Test 2^4423 - 1. Probably prime Number. Used Time:1.412644 Input:9941 Begin Test 2^9941 - 1. Probably prime Number. Used Time:11.602822 Input:9689 Begin Test 2^9689 - 1. Probably prime Number. Used Time:10.859898 Input:0

无心人 · 发表于 2008-7-16 15:54:08

可以看到平均每10秒能测试出一个3000位十进制整数的素性现在问题是3000位的四生素数如果用 a * b# + C, a = 1..999999999999, b = 5000..9999, C是所有10^7以下的四生素数的集合请问，如果以10秒做一次素性筛选大概多长时间能找到一组通过素性筛选的四生素数可能组合？

gxqcn · 发表于 2008-7-16 16:19:26

GMP 采用的是5次随机 Miller-Rabin 测试，强度并不够， HugeCalc 采用了更高强度的素性检测。注：65# 的测试数据都比较特殊，HugeCalc 可以自动切换高速算法检测。

无心人 · 发表于 2008-7-16 16:56:26

我不要强度，我要速度成万的测试，速度第一另外，不要考虑特殊形式因为就是为了测试普通数字

无心人 · 发表于 2008-7-17 16:37:01

带状态保存的搜索大四生素数程序最大到(2^31-1) * 9999# + 10^8

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <sys/time.h>
#include "gmp.h"
unsigned long primes[1280]; //小于10000素数，实际上有1229个
unsigned long quadprimes[4800]; //小于10^8四生素数，实际上有4768组
unsigned long sieve[30030];
unsigned long template[30030];
unsigned long A, B, C;
void
InitPrimes (void)
{
primes[0] = 2;
primes[1] = 3;
primes[2] = 5;
primes[3] = 7;
primes[4] = 11;
primes[5] = 13;
primes[6] = 17;
primes[7] = 19;
primes[8] = 23;
primes[9] = 29;
primes[10] = 31;
primes[11] = 37;
primes[12] = 41;
primes[13] = 43;
primes[14] = 47;
primes[15] = 53;
primes[16] = 59;
primes[17] = 61;
primes[18] = 67;
primes[19] = 71;
primes[20] = 73;
primes[21] = 79;
primes[22] = 83;
primes[23] = 89;
primes[24] = 97;
}
void
InitSieve (void)
{
unsigned long i;
for (i = 1; i <= 9999; i += 2)
sieve[i] = 1;
}
void
Sieve (void)
{
unsigned long i, j, start;
for (i = 0; i <= 24; i++)
{
start = primes[i] * primes[i];
for (j = start; j <= 9999; j += 2 * primes[i])
sieve[j] = 0;
}
}
void
StoreSieveResult (void)
{
unsigned long i, k;
k = 25;
for (i = 101; i <= 9999; i += 2)
if (sieve[i])
primes[k++] = i;
for (; k < 1280; k++)
primes[k] = 0;
primes[1229] = 10007; //以这个素数结束，否则下面的四生素数筛在接近10^8会超越1229素数的范围
}
void
InitQuadPrimes (void)
{
unsigned long i, k = 0;
for (i = 0; i < 1229; i++)
if ((primes[i + 1] - primes[i] == 2) && (primes[i + 2] - primes[i] == 6)
&& (primes[i + 3] - primes[i] == 8))
quadprimes[k++] = primes[i]; //共12组
}
void
InitTemplate (void)
{
unsigned long i;
for (i = 1; i < 30030; i += 2)
template[i] = 1;
for (i = 0; i < 30030; i += 2)
template[i] = 0;
for (i = 3; i < 30030; i += 6)
template[i] = 0;
for (i = 5; i < 30030; i += 10)
template[i] = 0;
for (i = 7; i < 30030; i += 14)
template[i] = 0;
for (i = 11; i < 30030; i += 22)
template[i] = 0;
for (i = 13; i < 30030; i += 26)
template[i] = 0;
}
void
SearchQuadSieve (void)
{
unsigned long start, end;
unsigned long i, j, k, m, s_start;
InitTemplate ();
//首次筛过程
start = 0;
end = 30030;
memcpy (sieve, template, 30030 * sizeof (unsigned long));
for (j = 6; primes[j] * primes[j] <= end; j++)
{
s_start = primes[j] * primes[j];
for (k = s_start; k < 30030; k += 2 * primes[j])
sieve[k] = 0;
}
m = 12;
for (i = 10001; i < 30030; i += 10)
if ((sieve[i]) && (sieve[i + 2]) && (sieve[i + 6]) && (sieve[i + 8]))
quadprimes[m++] = i;
//完整筛过程，经测试剩余的数字是99999900-99999999，不存在四生素数
for (i = 1; i < 100000000 / 30030; i++)
{
start = i * 30030;
end = start + 30030;
memcpy (sieve, template, 30030 * sizeof (unsigned long));
for (j = 6; primes[j] * primes[j] <= start; j++)
{
s_start = (start / primes[j]) * primes[j] + primes[j];
if (s_start % 2 == 0)
s_start += primes[j];
s_start -= start;
for (k = s_start; k < 30030; k += 2 * primes[j])
sieve[k] = 0;
}
for (; primes[j] * primes[j] < end; j++)
{
s_start = primes[j] * primes[j] - start;
for (k = s_start; k < 30030; k += 2 * primes[j])
sieve[k] = 0;
}
s_start = (start / 30) * 30 + 11;
if (s_start < start)
s_start += 30;
s_start -= start;
for (k = s_start; k < 30030; k += 30)
if ((sieve[k]) && (sieve[k + 2]) && (sieve[k + 6]) && (sieve[k + 8]))
quadprimes[m++] = start + k;
}
}
void
Init (void)
{
struct timeval start, end;
double timeuse;
gettimeofday (&start, NULL);
InitPrimes ();
InitSieve ();
Sieve ();
StoreSieveResult ();
InitQuadPrimes ();
SearchQuadSieve ();
gettimeofday (&end, NULL);
timeuse =
1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
timeuse /= 1000000;
printf ("Time use:%8.6f: %\n", timeuse);
}
void
product (mpz_t r, unsigned long p[], unsigned long max)
{
unsigned long i;
mpz_set_ui (r, 1);
if (max < 2)
return;
for (i = 0; primes[i] <= max; i++)
mpz_mul_ui (r, r, p[i]);
}
//状态文件分三行，分别存储A, B, C的数值
void
LoadStatus (void)
{
FILE *fstatus;
fstatus = fopen ("quadprimes.status", "r+t");
if (fstatus == NULL)
{
printf ("读取状态文件失败。\n");
exit (1);
}
fscanf (fstatus, "%u%u%u", &A, &B, &C);
fclose (fstatus);
}
void
SaveStatus (void)
{
FILE *fstatus;
fstatus = fopen ("quadprimes.status", "w+t");
if (fstatus == NULL)
{
printf ("存储状态文件失败。\n");
exit (1);
}
fprintf (fstatus, "%u\n%u\n%u\n", A, B, C);
fclose (fstatus);
}
void
SaveLog (void)
{
FILE *flog;
flog = fopen ("quadprimes.log", "a+t");
if (flog == NULL)
{
printf ("存储日志文件失败。\n");
exit (1);
}
fprintf (flog, "%u %u %u\n", A, B, C);
printf ("%u %u %u\n", A, B, C);
fclose (flog);
}
unsigned long SearchQuadPrime(unsigned long p)
{
unsigned long k;
for (k = 0; k < 4768; k ++)
if (quadprimes[k] >= p)
return (k);
return (0);
}
void
Test (void)
{
unsigned long m, k;
mpz_t b0, b, p1, p2, p6, p8;
LoadStatus ();
mpz_init (b0);
mpz_init(b);
mpz_init (p1);
mpz_init (p2);
mpz_init (p6);
mpz_init (p8);
product (b0, primes, B);
k = SearchQuadPrime(C);
C = quadprimes[k];
mpz_mul_ui (b, b0, A);
mpz_add_ui (p1, b, C);
m = 0;
for (;;)
{
if (mpz_probab_prime_p (p1, 1))
{
mpz_add_ui (p2, p1, 2);
if (mpz_probab_prime_p (p2, 1))
{
mpz_add_ui (p6, p1, 6);
if (mpz_probab_prime_p (p6, 1))
{
mpz_add_ui (p8, p1, 8);
if (mpz_probab_prime_p (p8, 1))
SaveLog ();
}
}
}
//
m ++;
if (m > 100)
{
SaveStatus ();
m = 0;
}
k ++;
if (k >= 4768)
{
k = SearchQuadPrime(B);
C = quadprimes[k];
A ++;
if (A == 0)
{
SaveLog();
exit(2);
}
else
{
mpz_add(b, b, b0);
mpz_add_ui(p1, b, C);
}
}
else
{
C = quadprimes[k];
mpz_add_ui(p1, b, C);
}
}
mpz_clear (b);
mpz_clear(b0);
mpz_clear (p1);
mpz_clear (p2);
mpz_clear (p6);
mpz_clear (p8);
}
int main (void)
{
unsigned long a, c;
Init ();
Test ();
// for (a = 0; a < 4800; a++)
// if (quadprimes[a])
// printf ("%6u ", quadprimes[a]);
// printf ("\n");
return 0;
}

复制代码

无心人 · 发表于 2009-1-13 21:22:09

上面这个程序我忘记是否是调通的了总想拿这个来算下

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

[讨论] 估计下找到10^100以上的一个四生素数的工作量