数学研发论坛

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

[讨论] 10秒内计算出1亿素数表, 不输出

[复制链接]
 楼主| 发表于 2008-3-21 16:37:15 | 显示全部楼层

80字节拷贝函数, 速度4.3GB/s, P4Xeon2.0 DDR266

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
        }
}
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-22 09:16:53 | 显示全部楼层

P4 Xeon2/1G/DDR 570ms/如何优化?

#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 15:36:26 | 显示全部楼层
已优化到530ms
在做汇编化
:)
争取到400
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-22 16:53:43 | 显示全部楼层
呵呵
还是对调用栈糊涂
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-3-24 12:59:46 | 显示全部楼层
看了一下楼主的代码,使用的是基于BYTE数组的筛法,对于每一个sqrt(n)以内的小素数(2,3等除外),将筛去素数的所有奇数倍。
未测试之前,估计楼主程序的速度应该大于我在47#楼给出的第一个版本,而慢于第二个版本。但实际测试结果却令人大跌眼镜。楼主的版本竟然慢于我的没有使用任何技巧的第一个版本。
    因为楼主的程序功能比较简单,在筛完后没有统计素数的个数,为了便于比较,我对我的程序作了一点修改(删除统计素数个数的过程,并改用了高精度计时器),以下是测试结果(PIV2.6G, 768M RAM, 取5次做小值)。

n      Sieve8          siftPrime1           siftPrime2
1亿     502ms           375ms                129ms
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-3-24 13:16:52 | 显示全部楼层
楼主的程序中,汇编代码占用了很大的比重,但效果却不太理想。

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

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

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

2. 对最核心部分采用汇编化,主要手段是尽量使用寄存器而不是内存,降低指令数,降低内存读写(特别是写)次数,预计速度至少可提高15%以上。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-24 13:21:18 | 显示全部楼层
可是我有把素数保存的
如果不保存是300ms多
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-3-24 13:31:50 | 显示全部楼层
怪我没看仔细,怪不得测试结果另我不得其解呢?下面的代码保存了素数
  1.         for (i = 1; i < LCM; i +=2)
  2.         if (WorkCache[i])
  3.                         Prime[++CurrPos] = base + i;
  4.         LowSqrt = HighSqrt;
复制代码
对你的代码修改后(仅仅首次调用OneSieve时,保存素数)进行测试,多次测试的最小运行时间是245ms.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-24 15:54:57 | 显示全部楼层



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

那个最小时间是可以通过置进程为高优先, 停掉任何前台程序得到的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-3-24 17:07:40 | 显示全部楼层
原帖由 无心人 于 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。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2020-11-30 12:28 , Processed in 0.072430 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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