数学研发论坛

 找回密码
 欢迎注册
查看: 301|回复: 3

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

[复制链接]
发表于 2019-9-11 17:31:55 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
实际从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的差的相邻结果

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <gmp.h>
  5. #include <time.h>

  6. #define step (2*3*5*7*11*13)
  7. #define slice 333000
  8. #define i64 long long
  9. #define u64 unsigned long long
  10. #define i32 long
  11. #define u32 unsigned long

  12. char temp[step];
  13. char work[step];
  14. char n4[10000];
  15. u32 prime[1300];
  16. u32 primeCount = 0;
  17. i64 prev, gap, cnt;
  18. i64 base12, offset, curr;

  19. mpz_t base67, start, rem, test;

  20. void prepareTemp( void )
  21. {
  22.         u32 smallPrime[] = {2, 3, 5, 7, 11, 13};
  23.         for (i32 i = 0; i < step; i += 2) temp[i] = 0;
  24.         for (i32 i = 1; i < step; i += 2) temp[i] = 1;
  25.         for (i32 j = 1; j < 6; j ++)
  26.                 for (i32 i = smallPrime[j]; i < step; i += 2* smallPrime[j])  
  27.                         temp[i] = 0;
  28. }

  29. void copyToWork( void )
  30. {
  31.         memcpy(work, temp, step);
  32. }

  33. void preparePrime( void )
  34. {
  35.         char n2[100];
  36.         for (i32 i = 0; i < 100; i+=2)  n2[i] = 0;
  37.         for (i32 i = 3; i < 100; i+=2)  n2[i] = 1;
  38.         n2[2] = 1; n2[9] = 0;
  39.         for (i32 i = 0; i < 10; i ++)  if (n2[i] == 1) prime[primeCount++] = i;
  40.         for (i32 i = 1; i < primeCount; i ++)
  41.         {
  42.                 i32 p = prime[i];
  43.         i32 t = p * p;
  44.         while (t < 100)
  45.                 {
  46.           n2[t] = 0;
  47.           t += 2 * p;
  48.                 }
  49.         }
  50.        
  51.         for (i32 i = 11; i < 100; i += 2)  if (n2[i] == 1) prime[primeCount++] = i;

  52.         for (i32 i = 4; i < 10000; i+=2)  n4[i] = 0;
  53.         for (i32 i = 3; i < 10000; i+=2)  n4[i] = 1;
  54.         for (i32 i = 1; i < primeCount; i ++)
  55.         {
  56.                 i32 p = prime[i];
  57.         i32 t = p * p;
  58.         while (t < 10000)
  59.                 {
  60.                         n4[t] = 0;
  61.                         t += 2 * p;
  62.                 }
  63.         }
  64.         for (i32 i = 101; i < 10000; i += 2)  if (n4[i] == 1) prime[primeCount++] = i;
  65. }

  66. void init( u32 index )
  67. {
  68.         //base67 = 10^67 - 10
  69.         mpz_init(base67);
  70.         mpz_ui_pow_ui(base67, 10, 67);
  71.         mpz_sub_ui(base67, base67, 10);
  72.        
  73.         mpz_init2(start, 256);
  74.         mpz_init2(test, 256);
  75.         mpz_init2(rem, 256);
  76.         preparePrime();
  77.         prepareTemp();
  78.         base12 = index * slice * step;
  79.        
  80.         mpz_add_ui(base67, base67, base12);
  81.         mpz_sub_ui(test, base67, 1);
  82.         for (i32 i = 1; i < step; i +=2)
  83.         {
  84.                 if (mpz_probab_prime_p(test, 10) >= 1)
  85.                 {
  86.                         mpz_sub(test, test, base67);
  87.                         prev = mpz_get_si(test);
  88.                         break;
  89.                 }
  90.                 mpz_sub_ui(test, test, 2);
  91.         }
  92.         cnt = 0;
  93.         gap = 0;
  94.         printf("from 10^67 - 10 + %lld ...\n", base12);
  95. }

  96. void check( i64 i )
  97. {
  98.         {
  99.                 curr = i + offset;
  100.                 gap = curr - prev;
  101.                 if ( gap >= 2000 )
  102.                 {
  103.                         //10^67 - 10 + base12 + ...
  104.                         printf("%lld, %lld, %lld\n", prev, curr, gap);
  105.                 }
  106.                 prev = curr;
  107.                 cnt ++;
  108.         }
  109. }

  110. void oneStep( mpz_t start )
  111. {
  112.         copyToWork();
  113.         for (i32 i = 6; i < primeCount; i ++)
  114.         {
  115.                 u32 r = mpz_mod_ui(rem, start, prime[i]);
  116.                 r = prime[i] - r;
  117.                 if (r%2==0) r += prime[i];
  118.                 while (r < step)
  119.                 {
  120.                         work[r] = 0;
  121.                         r += 2 * prime[i];
  122.                 }
  123.         }
  124.        
  125.         for (i32 i = 1; i < step; i += 2)
  126.                 if (work[i] == 1)
  127.                 {
  128.                         mpz_add_ui(test, start, i);
  129.                         if (mpz_millerrabin(test, 1) == 0)
  130.                                 work[i] = 0;
  131.                 }
  132. /*       
  133.         int b = 0;
  134.         printf("----------------------------\nCheck work..........\n");
  135.         for (i32 i = 1; i < step; i += 2)
  136.         {
  137.                
  138.                 if (work[i] == 1)
  139.                 {
  140.                         printf("(%d, %d)", b, i);
  141.                         b ++;
  142.                         if (b >= 10) break;
  143.                 }
  144.         }
  145.         printf("\n--------------------------\n");
  146. */       
  147.         for (i32 i = 1; i < step; i += 2)
  148.                 if (work[i] == 1)
  149.                         check( i );
  150.                
  151.        
  152.                
  153. }

  154. void circle( u32 index )
  155. {
  156.        
  157.         for (u32 i = 0; i < slice; i ++)
  158.         {
  159.                 offset = i * step;
  160.                 mpz_add_ui(start, base67, offset);
  161.                 //printf("i = %d, offset = %lld\n", i, offset);
  162.                 //gmp_printf("start = %Zd\n", start);
  163.                 oneStep( start );
  164.         }
  165. }

  166. i32 main( int argc, char * argv[] )
  167. {
  168.         u32 index = 0;
  169.         time_t timer;
  170.         char buffer[80];
  171.         struct tm * timeinfo;
  172.        
  173.         if (argc > 1)
  174.                 index = strtol(argv[1], NULL, 10);
  175.         else
  176.                 exit(1);
  177.         init(index);
  178.         //gmp_printf("start from %Zd\n", base67);
  179.        
  180.         time(&timer);
  181.         timeinfo = localtime (&timer);
  182.         strftime (buffer,80,"start time: %D %T.\n", timeinfo);
  183.         puts(buffer);
  184.        
  185.         circle(index);

  186.         time(&timer);
  187.         timeinfo = localtime (&timer);
  188.         strftime (buffer,80,"\nend time: %D %T.\n", timeinfo);
  189.         puts(buffer);
  190.        
  191.         //mpz_add_ui(test, base67, curr);
  192.         //gmp_printf("end with %Zd\n", test);
  193.         return 1;
  194. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-9-12 21:43:02 | 显示全部楼层
  1. #define i32 long
  2. #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位寄存器,这叫绿色环保意识
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2019-10-19 18:41 , Processed in 0.095702 second(s), 16 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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