数学研发论坛

 找回密码
 欢迎注册
楼主: 无心人

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

[复制链接]
发表于 2008-7-16 12:38:07 | 显示全部楼层
原帖由 无心人 于 2008-7-16 10:28 发表


GxQ尝试过多大的数字的素性测试
具体时间是多少?


没具体测试 HugeCalc 素性判定的能力,
应该说是足够大的,可以容纳当前已知的最大素数。

其算法还可进一步优化,留待下一版吧。。。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-7-16 13:04:33 | 显示全部楼层


我就想知道进行一个10000位数字的一次测试需要多少秒?
是一次测试,即单独一个素数的测试
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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秒做一次素性筛选
大概多长时间能找到一组通过素性筛选的四生素数可能组合?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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
  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <string.h>
  4. #include <math.h>
  5. #include <assert.h>
  6. #include <sys/time.h>
  7. #include "gmp.h"

  8. unsigned long primes[1280];        //小于10000素数,实际上有1229个
  9. unsigned long quadprimes[4800];        //小于10^8四生素数,实际上有4768组
  10. unsigned long sieve[30030];
  11. unsigned long template[30030];
  12. unsigned long A, B, C;

  13. void
  14. InitPrimes (void)
  15. {
  16.   primes[0] = 2;
  17.   primes[1] = 3;
  18.   primes[2] = 5;
  19.   primes[3] = 7;
  20.   primes[4] = 11;
  21.   primes[5] = 13;
  22.   primes[6] = 17;
  23.   primes[7] = 19;
  24.   primes[8] = 23;
  25.   primes[9] = 29;
  26.   primes[10] = 31;
  27.   primes[11] = 37;
  28.   primes[12] = 41;
  29.   primes[13] = 43;
  30.   primes[14] = 47;
  31.   primes[15] = 53;
  32.   primes[16] = 59;
  33.   primes[17] = 61;
  34.   primes[18] = 67;
  35.   primes[19] = 71;
  36.   primes[20] = 73;
  37.   primes[21] = 79;
  38.   primes[22] = 83;
  39.   primes[23] = 89;
  40.   primes[24] = 97;
  41. }

  42. void
  43. InitSieve (void)
  44. {
  45.   unsigned long i;
  46.   for (i = 1; i <= 9999; i += 2)
  47.     sieve[i] = 1;
  48. }

  49. void
  50. Sieve (void)
  51. {
  52.   unsigned long i, j, start;
  53.   for (i = 0; i <= 24; i++)
  54.     {
  55.       start = primes[i] * primes[i];
  56.       for (j = start; j <= 9999; j += 2 * primes[i])
  57.         sieve[j] = 0;
  58.     }
  59. }

  60. void
  61. StoreSieveResult (void)
  62. {
  63.   unsigned long i, k;
  64.   k = 25;
  65.   for (i = 101; i <= 9999; i += 2)
  66.     if (sieve[i])
  67.       primes[k++] = i;
  68.   for (; k < 1280; k++)
  69.     primes[k] = 0;
  70.   primes[1229] = 10007;                //以这个素数结束,否则下面的四生素数筛在接近10^8会超越1229素数的范围
  71. }

  72. void
  73. InitQuadPrimes (void)
  74. {
  75.   unsigned long i, k = 0;
  76.   for (i = 0; i < 1229; i++)
  77.     if ((primes[i + 1] - primes[i] == 2) && (primes[i + 2] - primes[i] == 6)
  78.         && (primes[i + 3] - primes[i] == 8))
  79.       quadprimes[k++] = primes[i];        //共12组
  80. }

  81. void
  82. InitTemplate (void)
  83. {
  84.   unsigned long i;
  85.   for (i = 1; i < 30030; i += 2)
  86.     template[i] = 1;
  87.   for (i = 0; i < 30030; i += 2)
  88.     template[i] = 0;
  89.   for (i = 3; i < 30030; i += 6)
  90.     template[i] = 0;
  91.   for (i = 5; i < 30030; i += 10)
  92.     template[i] = 0;
  93.   for (i = 7; i < 30030; i += 14)
  94.     template[i] = 0;
  95.   for (i = 11; i < 30030; i += 22)
  96.     template[i] = 0;
  97.   for (i = 13; i < 30030; i += 26)
  98.     template[i] = 0;
  99. }

  100. void
  101. SearchQuadSieve (void)
  102. {
  103.   unsigned long start, end;
  104.   unsigned long i, j, k, m, s_start;

  105.   InitTemplate ();

  106. //首次筛过程
  107.   start = 0;
  108.   end = 30030;
  109.   memcpy (sieve, template, 30030 * sizeof (unsigned long));

  110.   for (j = 6; primes[j] * primes[j] <= end; j++)
  111.     {
  112.       s_start = primes[j] * primes[j];
  113.       for (k = s_start; k < 30030; k += 2 * primes[j])
  114.         sieve[k] = 0;
  115.     }

  116.   m = 12;
  117.   for (i = 10001; i < 30030; i += 10)
  118.     if ((sieve[i]) && (sieve[i + 2]) && (sieve[i + 6]) && (sieve[i + 8]))
  119.       quadprimes[m++] = i;

  120. //完整筛过程,经测试剩余的数字是99999900-99999999,不存在四生素数
  121.   for (i = 1; i < 100000000 / 30030; i++)
  122.     {
  123.       start = i * 30030;
  124.       end = start + 30030;
  125.       memcpy (sieve, template, 30030 * sizeof (unsigned long));

  126.       for (j = 6; primes[j] * primes[j] <= start; j++)
  127.         {
  128.           s_start = (start / primes[j]) * primes[j] + primes[j];
  129.           if (s_start % 2 == 0)
  130.             s_start += primes[j];
  131.           s_start -= start;
  132.           for (k = s_start; k < 30030; k += 2 * primes[j])
  133.             sieve[k] = 0;
  134.         }

  135.       for (; primes[j] * primes[j] < end; j++)
  136.         {
  137.           s_start = primes[j] * primes[j] - start;
  138.           for (k = s_start; k < 30030; k += 2 * primes[j])
  139.             sieve[k] = 0;
  140.         }

  141.       s_start = (start / 30) * 30 + 11;
  142.       if (s_start < start)
  143.         s_start += 30;
  144.       s_start -= start;

  145.       for (k = s_start; k < 30030; k += 30)
  146.         if ((sieve[k]) && (sieve[k + 2]) && (sieve[k + 6]) && (sieve[k + 8]))
  147.           quadprimes[m++] = start + k;

  148.     }
  149. }

  150. void
  151. Init (void)
  152. {
  153.   struct timeval start, end;
  154.   double timeuse;

  155.   gettimeofday (&start, NULL);
  156.   InitPrimes ();
  157.   InitSieve ();
  158.   Sieve ();
  159.   StoreSieveResult ();

  160.   InitQuadPrimes ();
  161.   SearchQuadSieve ();
  162.   gettimeofday (&end, NULL);

  163.   timeuse =
  164.     1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
  165.   timeuse /= 1000000;

  166.   printf ("Time use:%8.6f:  %\n", timeuse);

  167. }

  168. void
  169. product (mpz_t r, unsigned long p[], unsigned long max)
  170. {
  171.   unsigned long i;
  172.   mpz_set_ui (r, 1);

  173.   if (max < 2)
  174.     return;

  175.   for (i = 0; primes[i] <= max; i++)
  176.     mpz_mul_ui (r, r, p[i]);
  177. }

  178. //状态文件分三行,分别存储A, B, C的数值
  179. void
  180. LoadStatus (void)
  181. {
  182.   FILE *fstatus;
  183.   fstatus = fopen ("quadprimes.status", "r+t");
  184.   if (fstatus == NULL)
  185.     {
  186.       printf ("读取状态文件失败。\n");
  187.       exit (1);
  188.     }
  189.   fscanf (fstatus, "%u%u%u", &A, &B, &C);
  190.   fclose (fstatus);
  191. }

  192. void
  193. SaveStatus (void)
  194. {
  195.   FILE *fstatus;
  196.   fstatus = fopen ("quadprimes.status", "w+t");
  197.   if (fstatus == NULL)
  198.     {
  199.       printf ("存储状态文件失败。\n");
  200.       exit (1);
  201.     }
  202.   fprintf (fstatus, "%u\n%u\n%u\n", A, B, C);
  203.   fclose (fstatus);
  204. }

  205. void
  206. SaveLog (void)
  207. {
  208.   FILE *flog;
  209.   flog = fopen ("quadprimes.log", "a+t");
  210.   if (flog == NULL)
  211.     {
  212.       printf ("存储日志文件失败。\n");
  213.       exit (1);
  214.     }
  215.   fprintf (flog, "%u    %u    %u\n", A, B, C);
  216.   printf ("%u    %u    %u\n", A, B, C);
  217.   fclose (flog);

  218. }

  219. unsigned long SearchQuadPrime(unsigned long p)
  220. {
  221.   unsigned long k;
  222.   for (k = 0; k < 4768; k ++)
  223.      if (quadprimes[k] >= p)
  224.           return (k);
  225.     return (0);
  226. }

  227. void
  228. Test (void)
  229. {
  230.   unsigned long m, k;
  231.   mpz_t b0, b, p1, p2, p6, p8;
  232.   LoadStatus ();
  233.   mpz_init (b0);
  234.   mpz_init(b);
  235.   mpz_init (p1);
  236.   mpz_init (p2);
  237.   mpz_init (p6);
  238.   mpz_init (p8);

  239.   product (b0, primes, B);
  240.   k = SearchQuadPrime(C);
  241.   C = quadprimes[k];
  242.   mpz_mul_ui (b, b0, A);
  243.   mpz_add_ui (p1, b, C);
  244.   m = 0;
  245.   for (;;)
  246.     {
  247.       if (mpz_probab_prime_p (p1, 1))
  248.         {
  249.           mpz_add_ui (p2, p1, 2);
  250.           if (mpz_probab_prime_p (p2, 1))
  251.             {
  252.               mpz_add_ui (p6, p1, 6);
  253.               if (mpz_probab_prime_p (p6, 1))
  254.                 {
  255.                   mpz_add_ui (p8, p1, 8);
  256.                   if (mpz_probab_prime_p (p8, 1))
  257.                  SaveLog ();
  258.                 }
  259.             }
  260.          }

  261.     //
  262.         m ++;
  263.         if (m > 100)
  264.           {
  265.           SaveStatus ();
  266.           m = 0;
  267.           }
  268.         k ++;
  269.         if (k >= 4768)
  270.           {
  271.           k = SearchQuadPrime(B);
  272.           C = quadprimes[k];
  273.           A ++;
  274.           if (A == 0)
  275.              {
  276.             SaveLog();
  277.             exit(2);
  278.              }
  279.           else
  280.              {
  281.             mpz_add(b, b, b0);
  282.             mpz_add_ui(p1, b, C);
  283.              }
  284.           }
  285.         else
  286.           {
  287.           C =  quadprimes[k];
  288.           mpz_add_ui(p1, b, C);
  289.           }
  290.     }
  291.   mpz_clear (b);
  292.   mpz_clear(b0);
  293.   mpz_clear (p1);
  294.   mpz_clear (p2);
  295.   mpz_clear (p6);
  296.   mpz_clear (p8);
  297. }

  298. int main (void)
  299. {
  300.   unsigned long a, c;

  301.   Init ();

  302.   Test ();
  303. //  for (a = 0; a < 4800; a++)
  304. //    if (quadprimes[a])
  305. //      printf ("%6u ", quadprimes[a]);
  306. //  printf ("\n");
  307.   return 0;
  308. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-1-13 21:22:09 | 显示全部楼层


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

总想拿这个来算下
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-12-9 22:54 , Processed in 0.102112 second(s), 15 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表