求10的10次方以内的素数，怎么能在几小时内算出来结果？

glasslion

原帖由 ll3021359 于 2009-4-22 06:19 发表
有个蒙特卡罗算法专门素数测试，我们学校的OJ上有一道题要求计算10的9次方之内的素数，一个大牛仅在7秒内搞定，该题限时是30秒

估计是是判断一个数是素数
要找你说的计算10的9次方之内的素数，输出大概有上百MB,光一个IO就超时了
而且估计没一个学校的OJ经得起这种输出

tyc611

主要是只处理奇数，另外，我发现在计算倍数时用乘法代替我先前使用的加法会加些，主要代码如下（虽然用C++写的，但也还是容易看懂的）

1. 
   
2.     IntType root = IntType(std::sqrt(double(m_maxNum)));
   
3.     for (IntType i = 3; i <= root; i += 2) { // 只处理奇数
   
4.         if (!GetTag(i)) {
   
5.             continue;
   
6.         }
   
7. 
   
8.         IntType limit = m_maxNum / i;  // 防止无符号类型的溢出
   
9.         for (IntType k = i; k <= limit; k += 2) { // 只处理奇数
   
10.             IntType m = i * k;
    
11.             if (GetTag(m)) {
    
12.                 ClearTag(m);
    
13.                 --m_primesCount;
    
14.             }
    
15.         }
    
16.     }
    

复制代码

在你的代码中，下面这段我没看明白

1. 
   
2.         for (j=(3*p)/2; j<(N)/2; j+=p)
   
3.         {
   
4.             if (primes[j/32] & (0x1UL)<<(j%32))
   
5.                 primes[j/32] &= ~((0x1UL)<<(j%32));
   
6.         }
   

复制代码

这里的两个除2是什么用途？

tyc611

查表更慢了=

F:\tmp\Primes>..\timing.exe EratosthenesPrimes_example.exe
2009-04-29 00:09:31 running...
Total number of primes: 203280221
2009-04-29 00:11:15 ...finished

Kernel Time: 140ms
User   Time: 1m43s140ms
Total  Time: 1m43s281ms

log271828182845

1) 首先偶数都不是素数;
2)所有的奇数可以写成 p = (2i+1)形式
3) 任何素数当中的非素数都可以写成至少两个奇数的乘积:
    3.1  p = m * n
    3.2  p > m > 2,  p > n > 2

    令  
          p  = 2i+1
          m = 2j+1
          n  = kj+1
         i > j >= k >= 1
           
=>  2i+1 = (2j+1)*(2k+1)
=>  2i +1 = 4jk +  2j + 2k + 1
=> i = 2jk + j + k                                       -----------------(I)  
=> i+1 = 2jk + j + k + 1 = (j+1)(k+1)
    令 p < 10^10
        则 2i + 1 < 10^10
             i <=  5*10^9
=>        1 <= k <=j < i  < 5*10^9   
=>           (j+1)^2         < 5*10^9  
=>           j+1 <

由 (I) 和 (II)  

睡觉...............

dreamice

原帖由 ll3021359 于 2009-4-22 13:19 发表
有个蒙特卡罗算法专门素数测试，我们学校的OJ上有一道题要求计算10的9次方之内的素数，一个大牛仅在7秒内搞定，该题限时是30秒

7秒内搞定的，肯定是以前算过的，不然不可能

handsomezhi

我也学习学习

redspider

楼太高，有没有人总结一下，究竟讨论出方法了没。概率的不算，不靠谱。

retuor

原帖由 tyc611 于 2009-4-28 23:49 发表

在你的代码中，下面这段我没看明白


        for (j=(3*p)/2; j<(N)/2; j+=p)
        {
            if (primes[j/32] & (0x1UL)<<(j%32))
                primes[j/32] &= ~((0x1UL)<<(j%32));
        }

这里的两个除2是什么用途？
  ...

如果要处理 N 个数，则位图大小是 N bit，但如果只处理奇数，则那些偶数可以不在位图中出现，位图的大小可为  N/2 bit，里面的信息都是关于奇数的。

位置：0 1 2 3 4
数值：1 3 5 7 9 ... ...

若当前素数为 p，则我要滤掉 2p, 3p, 4p, 5p ......

其中 2p, 4p 在我的位图里不出现，所以从 3p 开始，然后是 5p , 7p

对应位置为：3p/2, 5p/2, 7p/2

若 j=3p/2, 则 5p/2=j+p, 7p/2=j+2p

tyc611

got it

retuor

小改进：

若当前素数为 p, 过滤可从 p*p 开始，因为 3p, 5p, .... (p-2)p 其实都滤过的。