数学研发论坛

 找回密码
 欢迎注册
楼主: medie2005

[讨论] K-Smith numbers

[复制链接]
发表于 2008-12-1 12:44:42 | 显示全部楼层
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.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-12-1 12:49:35 | 显示全部楼层
下面的代码段 筛出 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.         np=primeTab[i];
  9.         dsum=getDigSum2(np);
  10.        
  11.         k= (base + np-1)/np;
  12.         if (  k<np)
  13.                 k=np;
  14.         m= k*np;
  15.         while (m <=end)
  16.         {
  17.                 g_valueArray[m-base ] /= np;
  18.                 g_digSumArray[m-base] += dsum;
  19.                 m+= np;
  20.         }
  21.        
  22.         npx=np;
  23.         npx_end=end/npx;
  24.         while (npx<=npx_end)
  25.         {
  26.                 npx *= np;
  27.                 m= (base + npx-1)/npx*npx;
  28.                 while (m <=end)
  29.                 {
  30.                         g_valueArray[m-base ] /= np;
  31.                         g_digSumArray[m-base] += dsum;
  32.                         m+= npx;
  33.                 }
  34.         }
  35.         i++;
  36. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-12-1 12:54:18 | 显示全部楼层
为了便于阅读,将输出 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. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-12-1 17:38:48 | 显示全部楼层
大哥厉害啊!学习了!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2021-2-25 21:16 , Processed in 0.056449 second(s), 15 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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