估计下找到10^100以上的一个四生素数的工作量 - 第7页 - 编程擂台 - 数学研发论坛

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

原帖由无心人于 2008-7-16 10:28 发表 http://bbs.emath.ac.cn/images/common/back.gif
:)

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; //小于10000素数，实际上有1229个
unsigned long quadprimes; //小于10^8四生素数，实际上有4768组
unsigned long sieve;
unsigned long template;
unsigned long A, B, C;

void
InitPrimes (void)
{
primes = 2;
primes = 3;
primes = 5;
primes = 7;
primes = 11;
primes = 13;
primes = 17;
primes = 19;
primes = 23;
primes = 29;
primes = 31;
primes = 37;
primes = 41;
primes = 43;
primes = 47;
primes = 53;
primes = 59;
primes = 61;
primes = 67;
primes = 71;
primes = 73;
primes = 79;
primes = 83;
primes = 89;
primes = 97;
}

void
InitSieve (void)
{
unsigned long i;
for (i = 1; i <= 9999; i += 2)
sieve = 1;
}

void
Sieve (void)
{
unsigned long i, j, start;
for (i = 0; i <= 24; i++)
{
   start = primes * primes;
   for (j = start; j <= 9999; j += 2 * primes)
sieve = 0;
}
}

void
StoreSieveResult (void)
{
unsigned long i, k;
k = 25;
for (i = 101; i <= 9999; i += 2)
if (sieve)
   primes = i;
for (; k < 1280; k++)
primes = 0;
primes = 10007; //以这个素数结束，否则下面的四生素数筛在接近10^8会超越1229素数的范围
}

void
InitQuadPrimes (void)
{
unsigned long i, k = 0;
for (i = 0; i < 1229; i++)
if ((primes - primes == 2) && (primes - primes == 6)
&& (primes - primes == 8))
   quadprimes = primes; //共12组
}

void
InitTemplate (void)
{
unsigned long i;
for (i = 1; i < 30030; i += 2)
template = 1;
for (i = 0; i < 30030; i += 2)
template = 0;
for (i = 3; i < 30030; i += 6)
template = 0;
for (i = 5; i < 30030; i += 10)
template = 0;
for (i = 7; i < 30030; i += 14)
template = 0;
for (i = 11; i < 30030; i += 22)
template = 0;
for (i = 13; i < 30030; i += 26)
template = 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 * primes <= end; j++)
{
   s_start = primes * primes;
   for (k = s_start; k < 30030; k += 2 * primes)
sieve = 0;
}

m = 12;
for (i = 10001; i < 30030; i += 10)
if ((sieve) && (sieve) && (sieve) && (sieve))
   quadprimes = 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 * primes <= start; j++)
{
   s_start = (start / primes) * primes + primes;
   if (s_start % 2 == 0)
   s_start += primes;
   s_start -= start;
   for (k = s_start; k < 30030; k += 2 * primes)
   sieve = 0;
}

   for (; primes * primes < end; j++)
{
   s_start = primes * primes - start;
   for (k = s_start; k < 30030; k += 2 * primes)
   sieve = 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) && (sieve) && (sieve) && (sieve))
   quadprimes = 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 <= max; i++)
mpz_mul_ui (r, r, p);
}

//状态文件分三行，分别存储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 >= 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;
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;
      A ++;
      if (A == 0)
         {
         SaveLog();
         exit(2);
         }
      else
         {
         mpz_add(b, b, b0);
         mpz_add_ui(p1, b, C);
         }
      }
   else
      {
      C =quadprimes;
      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)
//    printf ("%6u ", quadprimes);
//printf ("\n");
return 0;
}

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

:L

上面这个程序
我忘记是否是调通的了

总想拿这个来算下

页: 1 2 3 4 5 6 [7]

数学研发论坛's Archiver