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 #include #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。