K-Smith numbers

liangbch

55楼代码虽然有bug,但依然有正确的思想，我用其思想，改进了程序，速度提高1倍多，在计算10^8以内的smith素数时，用时从23秒提高到11秒。

55楼代码中，“j *= p;" 这句话隐藏了一个较难发现的bug. 如果 j<=end 且j*p>end. 但执行了j*=p后，j仍然有可能小于end,

例如：
    j=3251*3251= 10569001;
    j*=p 后，本应该是34359822251， 但由于整数溢出，j的值变为83883.

liangbch

下面的代码段 筛出 base 到end之间所有smith数

g_valueArray[  i  ]的初值为 base+i;
g_digSumArray [ i ]的初值为 0;

1. i=0;
   
2. while (primeTab[i]*primeTab[i]<=end && i<primeCount)
   
3. {
   
4.         DWORD k,m,np;
   
5.         DWORD npx;
   
6.         DWORD npx_end;
   
7.         BYTE dsum;
   
8. 
   
9.         np=primeTab[i];
   
10.         dsum=getDigSum2(np);
    
11.        
    
12.         k= (base + np-1)/np;
    
13.         if (  k<np)
    
14.                 k=np;
    
15.         m= k*np;
    
16.         while (m <=end)
    
17.         {
    
18.                 g_valueArray[m-base ] /= np;
    
19.                 g_digSumArray[m-base] += dsum;
    
20.                 m+= np;
    
21.         }
    
22.        
    
23.         npx=np;
    
24.         npx_end=end/npx;
    
25.         while (npx<=npx_end)
    
26.         {
    
27.                 npx *= np;
    
28.                 m= (base + npx-1)/npx*npx;
    
29.                 while (m <=end)
    
30.                 {
    
31.                         g_valueArray[m-base ] /= np;
    
32.                         g_digSumArray[m-base] += dsum;
    
33.                         m+= npx;
    
34.                 }
    
35.         }
    
36.         i++;
    
37. }

复制代码

liangbch

为了便于阅读，将输出 base到end之间的smith数的代码一并贴出

1. for (i=0;i<SIFT_BLOCK_LEN && i+base<=n;i++)
   
2. {
   
3.         DWORD x=i+base;
   
4.         WORD s1= getDigSum2(x);
   
5.         WORD s2= g_digSumArray[i];
   
6.         if (s1==0)
   
7.                continue;
   
8.                        
   
9.         if ( g_valueArray[i] !=1 && g_valueArray[i] != x)
   
10.                 s2+=getDigSum2(g_valueArray[i]);
    
11.         if (s1==s2)
    
12.         {
    
13.                 if (fp!=NULL)
    
14.                 {
    
15.                         numBuff[countInBuff++]=x;
    
16.                         if ( countInBuff*sizeof(DWORD)>=sizeof(numBuff) )
    
17.                         {
    
18.                                 fwrite(numBuff,sizeof(DWORD),countInBuff,fp);
    
19.                                 countInBuff=0;
    
20.                         }
    
21.                 }
    
22.                 count++;
    
23.         }
    
24. }

复制代码

suan1024

大哥厉害啊！学习了！