10秒内计算出1亿素数表, 不输出

无心人

void CopyCache(void)
{
        __asm
        {
        mov ecx, LCM
        sub ecx, 80
loop1:
        pshufd xmm0, BaseCache[ecx], 11100100b
        pshufd xmm1, BaseCache[ecx+16], 11100100b
        pshufd xmm2, BaseCache[ecx+32], 11100100b
        pshufd xmm3, BaseCache[ecx+48], 11100100b
        pshufd xmm4, BaseCache[ecx+64], 11100100b
        movdqa WorkCache[ecx], xmm0
        movdqa WorkCache[ecx+16], xmm1
        movdqa WorkCache[ecx+32], xmm2
        movdqa WorkCache[ecx+48], xmm3
        movdqa WorkCache[ecx+64], xmm4
        sub ecx, 80
        jnc loop1
        }
}

无心人

#include <stdio.h>
#include <stdlib.h>
#include < wtypes.h >

#define LCM (16*3*5*7*11*13)

char BaseCache[LCM], WorkCache[LCM];
unsigned int Prime[6000000];
unsigned int CurrPos; //当前素数表位置
unsigned int LowSqrt, HighSqrt; //起始位置和结束位置的平方根

unsigned int __fastcall intSieveSqrt(unsigned int n)
{
        __asm
        {
        mov ebx, n
        cvtsi2sd xmm0, ebx
        sqrtsd xmm0, xmm0
        cvtsd2si eax, xmm0
        mov ecx, eax
        mul ecx
        cmp eax, ebx
        jbe exit1
        sub ecx, 1
exit1:
    mov eax, ecx
    }
}

void initCache(void)
{
        __asm
        {
            mov eax, 01000100h
                movd xmm0, eax
        pshufd xmm0, xmm0, 00000000b
        mov ecx, LCM
                sub ecx, 16
loop1:
                movdqa BaseCache[ecx], xmm0
            sub ecx, 16
                jnc loop1

        mov ecx, LCM
                sub ecx, 3
loop2:
                mov byte ptr BaseCache[ecx], 0
            sub ecx, 6
                jnc loop2

        mov ecx, LCM
                sub ecx, 5
loop3:
                mov byte ptr BaseCache[ecx], 0
            sub ecx, 10
                jnc loop3
        
        mov ecx, LCM
                sub ecx, 7
loop4:
                mov byte ptr BaseCache[ecx], 0
            sub ecx, 14
                jnc loop4

        mov ecx, LCM
                sub ecx, 11
loop5:
                mov byte ptr BaseCache[ecx], 0
            sub ecx, 22
                jnc loop5

        mov ecx, LCM
                sub ecx, 13
loop6:
                mov byte ptr BaseCache[ecx], 0
            sub ecx, 26
                jnc loop6
        }
}

void CopyCache(void)
{
        __asm
        {
        mov ecx, LCM
        sub ecx, 80
loop1:
        pshufd xmm0, BaseCache[ecx], 11100100b
        pshufd xmm1, BaseCache[ecx+16], 11100100b
        pshufd xmm2, BaseCache[ecx+32], 11100100b
        pshufd xmm3, BaseCache[ecx+48], 11100100b
        pshufd xmm4, BaseCache[ecx+64], 11100100b
        movdqa WorkCache[ecx], xmm0
        movdqa WorkCache[ecx+16], xmm1
        movdqa WorkCache[ecx+32], xmm2
        movdqa WorkCache[ecx+48], xmm3
        movdqa WorkCache[ecx+64], xmm4
        sub ecx, 80
        jnc loop1
        }
}

void firstSieve(void)
{
        unsigned int i, j, k, k2, l, SieveIndex, p, p2;
        LowSqrt = intSieveSqrt(LCM - 1);
        CopyCache();
        Prime[0] = 2;
        Prime[1] = 3;
        Prime[2] = 5;
        Prime[3] = 7;
        Prime[4] = 11;
        Prime[5] = 13;
    j = 5;
//得到已筛出的初始素数， 小于17 * 17
        for (i = 17; i < (17 * 17); i += 2)
        {
      if (WorkCache[i])
                Prime[++j] = i;
        }
   
//利用新筛出的素数筛选, 最大应该283
        for (i = 6; i <= j; i ++)
        {
                p = Prime[i];
                p2 = p << 1;
        for (k = p*p; k < LCM; k+=p2)
           WorkCache[k] = 0;
        }
   
// <= Prime[j] * Prime[j]已筛
//从SieveIndex下个素数筛
        SieveIndex = j + 1;
    k = Prime[j] * Prime[j];
        for (i = Prime[SieveIndex - 1] + 2; i <= k; i += 2)
        {
      if (WorkCache[i])
                Prime[++j] = i;
        }

//最终筛，到sqrt(LCM)结束
    for (i = SieveIndex; Prime[i] <= LowSqrt; i ++)
        {
          k = Prime[i] * Prime[i];
          k2 = k << 1;
          for (l = k; l < LCM; l += k2)
                  WorkCache[l] = 0;
        }

//输出最后筛结果
        k = Prime[j] * Prime[j];
        for (i = k + 2; i < LCM; i += 2)
        {
      if (WorkCache[i])
                  Prime[++j] = i;
        }
        CurrPos = j;
}

void __fastcall OneSieve(unsigned int base)
{
    unsigned int i, j, k, l, p2;
        CopyCache();
        HighSqrt = intSieveSqrt(base + LCM - 1);
   
        for (i = 6; Prime[i] <= LowSqrt; i ++)
        {
                p2 = Prime[i] << 1;
                j = (base / Prime[i]) * Prime[i];
                if (j < base) j+= Prime[i];
                if (j % 2 == 0) j += Prime[i];
                j -= base;
                for (k = j; k < LCM; k += p2)
                        WorkCache[k] = 0;
        }

        l = i;
        for (i = l ;  Prime[i] <= HighSqrt; i ++)
        {
                p2 = Prime[i] << 1;
                j = Prime[i] * Prime[i] - base;
                for (k = j; k < LCM; k += p2)
                        WorkCache[k] = 0;
        }

        for (i = 1; i < LCM; i +=2)
        if (WorkCache[i])
                        Prime[++CurrPos] = base + i;
        LowSqrt = HighSqrt;
}

int main(void)
{
        UINT64 s_u64Frequency = 1;
    UINT64 s_u64Start, s_u64End;
        unsigned int i;
        QueryPerformanceFrequency((LARGE_INTEGER *)&s_u64Frequency );
        QueryPerformanceCounter((LARGE_INTEGER *)&s_u64Start );

        initCache();

        firstSieve();

        for (i = 1; LCM * i <= 100000000; i ++)
                OneSieve(i * LCM);

        QueryPerformanceCounter((LARGE_INTEGER *)&s_u64End );

        printf( "Elapsed time: %.3f ms\n",
                (double)(( s_u64End - s_u64Start ) * 1000.0 / (double)s_u64Frequency ));

        return 0;
}

代码见附件
Sieve8.rar (2.62 KB, 下载次数: 14)

2008-3-22 09:20 上传

点击文件名下载附件

无心人

已优化到530ms
在做汇编化
：）
争取到400

无心人

呵呵
还是对调用栈糊涂

liangbch

看了一下楼主的代码，使用的是基于BYTE数组的筛法，对于每一个sqrt(n)以内的小素数（2，3等除外），将筛去素数的所有奇数倍。
未测试之前，估计楼主程序的速度应该大于我在47#楼给出的第一个版本，而慢于第二个版本。但实际测试结果却令人大跌眼镜。楼主的版本竟然慢于我的没有使用任何技巧的第一个版本。
    因为楼主的程序功能比较简单，在筛完后没有统计素数的个数，为了便于比较，我对我的程序作了一点修改（删除统计素数个数的过程，并改用了高精度计时器），以下是测试结果（PIV2.6G， 768M RAM， 取5次做小值）。

n      Sieve8          siftPrime1           siftPrime2
1亿     502ms           375ms                129ms

liangbch

楼主的程序中，汇编代码占用了很大的比重，但效果却不太理想。

   优化的原则是有所为有所不为。通常情况下，首先应该优化算法，最后才是考虑使用汇编优化。实践表明，程序中20%代码占用了80%的运行时间，因此优化应该首先找出热点部分（比如多重循环中的最内层循环），然后将这少数的代码用汇编改写。当然对其他部分使用汇编代码，也能提高运行速度，但通常不这样做，这很不划算，恰如一首歌所言，"付出太多，得到太少，整天烦恼".

  受楼主启发，我决定对我的程序做进一步优化，主要从两个方面进行。

1. 将L1 prime table 采用 压缩的bit数组（对n1 到 n2 之间的数进行筛时，需要筛去2 to sqrt(n2)之间的所有素数的某些倍数,我这里把2 to sqrt(n2)之间的素数称为L1 Prime table），这样减少了内存的占用，从而提高了L2 cache的命中率，这个方法应该不会提速太多。

2. 对最核心部分采用汇编化，主要手段是尽量使用寄存器而不是内存，降低指令数，降低内存读写（特别是写）次数，预计速度至少可提高15%以上。

无心人

可是我有把素数保存的
如果不保存是300ms多

liangbch

怪我没看仔细，怪不得测试结果另我不得其解呢？下面的代码保存了素数

1.         for (i = 1; i < LCM; i +=2)
   
2.         if (WorkCache[i])
   
3.                         Prime[++CurrPos] = base + i;
   
4.         LowSqrt = HighSqrt;

复制代码

对你的代码修改后（仅仅首次调用OneSieve时，保存素数）进行测试，多次测试的最小运行时间是245ms.

无心人

是否能得到结论
如果二级缓存足够大
内存足够大
预先筛的素数越大越好？

那个最小时间是可以通过置进程为高优先， 停掉任何前台程序得到的

liangbch

原帖由 无心人 于 2008-3-24 15:54 发表



是否能得到结论
如果二级缓存足够大
内存足够大
预先筛的素数越大越好？

那个最小时间是可以通过置进程为高优先， 停掉任何前台程序得到的

预先筛的素数不是越大越好，它主要受限于L2 cache的大小。

   如果预先筛（2，3，5，7，11，13）这几个素数，  则采用BYTE数组，模板大小为30030，
   如果预先筛（2，3，5，7，11，13,17）这几个素数，则采用BYTE数组，模板大小为510510， 如果考虑小素数表所占的空间，这已经不适用 具有512K 或者更小的 L2 cache的电脑。
   如果预先筛（2，3，5，7，11，13,17，19）这几个素数，则采用BYTE数组，模板大小为9699690，这已经超出所有的现代计算机L2 cache的大小。

我的第三个版本，用的就是预先筛掉（2，3，5，7，11，13,17）这几个素数的算法，由于采用压缩的bit数组，实际模板的大小仅为17017 byte。