数学研发论坛

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

[讨论] K-Smith numbers

[复制链接]
发表于 2008-11-25 18:57:55 | 显示全部楼层
ksmith100000000.TXT (46.03 KB, 下载次数: 3)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-25 18:59:10 | 显示全部楼层


1小时14分钟
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-26 09:29:37 | 显示全部楼层
看来无心人你这个程序还是难敌C程序呀。到今日凌晨5点,我也完成了我的C程序。这个程序的复杂度近似于O(n),在PM 1.7上,计算10^8以内的smith数(不输出),不到10秒,好像是4.2秒,可是U盘忘带了。等明天吧,明天我贴出程序。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-26 09:40:44 | 显示全部楼层
不是不敌
是我初学
写不出高效的算法啊

如果C用的时间是1
Haskell通常在2-5
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-26 15:00:59 | 显示全部楼层


用筛法做,是最佳的选择了,呵呵
不过啊
haskell实现筛法,似乎有点挑战啊
不是haskell不行
是俺笨哦
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-26 15:42:14 | 显示全部楼层

回复 25# 无心人 的帖子

使用晒法,我在5楼就提到了,关于复杂度我也在5楼提到了。只是在实现的时候发现不需使用内存池和链表。
  另外,可以采用一些优化手段。我的做法是事先算出10000以下所有的数的数字和,用的时候查表。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-26 17:14:39 | 显示全部楼层
呵呵

这里我并不奢望多高等的算法
当学习haskell了
目前的估计是
30分钟
用的是逐个分解的算法
另外,算100000000个数字的数字和
是5秒
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-27 10:01:45 | 显示全部楼层
贴出我的解法:
  1.欲求n以内的所有smith数,首先需要求出sqrt(n)以内的所有素数。
  2.为了快速的求出数x的数字和,事先计算出10-10000的数字和,在需要的时候直接查表。
  3.需要对 n以内的每个数分解质因数,对x分解质因数,最笨的方法是用sqrt(x)以内的素数逐个试除,
我这里采用类似用筛法求素数的方法,平均下来,对n以内的每个数分解质因数的复杂度为log(log(n).
  1. #include "stdlib.h"
  2. #include "stdio.h"
  3. #include "math.h"
  4. #include "memory.h"

  5. typedef unsigned __int64 UINT64;
  6. typedef unsigned long DWORD;
  7. typedef unsigned short WORD;
  8. typedef unsigned char BYTE;


  9. //文件描述:使用最简单的方法筛素数

  10. // 一定范围内质数的个数
  11. //                           10                            4
  12. //                          100                           25
  13. //                        1,000                          168
  14. //                       10,000                        1,229
  15. //                      100,000                        9,592
  16. //                    1,000,000                       78,498
  17. //                   10,000,000                      664,579
  18. //                  100,000,000                    5,761,455
  19. //                1,000,000,000                   50,847,534
  20. //               10,000,000,000                  455,052,511
  21. //              100,000,000,000                4,118,054,813
  22. //            1,000,000,000,000               37,607,912,018
  23. //           10,000,000,000,000              346,065,536,839
  24. //          100,000,000,000,000            3,204,941,750,802
  25. //        1,000,000,000,000,000           29,844,570,422,669
  26. //       10,000,000,000,000,000          279,238,341,033,925
  27. //      100,000,000,000,000,000        2,623,557,157,654,233
  28. //    1,000,000,000,000,000,000       24,739,954,287,740,860
  29. //  估计x 范围内的质数的个数, 欧拉公式 pi(x)=x /( ln(x));

  30. // 经过测试,
  31. // 在1-1百万之间,相邻2个质数的最大差为114,这两个质数分别为492113,492227
  32. // 在1-1千万之间,相邻2个质数的最大差为154,这两个质数分别为4652353,4652507
  33. // 在1-3千万之间,相邻2个质数的最大差为210,这两个质数分别为20831323,20831533
  34. // 在1-1亿之间,  相邻2个质数的最大差为220,这两个质数分别为47326693,47326913


  35. //计算n以内素数个数的上限
  36. UINT64 maxPrimeCount(UINT64 n)
  37. {
  38.         double f= (double)n / log((double)n) * 10.0/9.0 +32;
  39.         return (UINT64)f;
  40. }

  41. DWORD Sift32(DWORD n,DWORD *primeTab)
  42. // 对0到n 之间(不包括n)的数进行筛选,返回质数的个数
  43. // n 不能大于L2 cache,否则效能很低
  44. //**********************************************
  45. {
  46.         BYTE *pBuff= new BYTE[n];
  47.         if (pBuff==NULL)
  48.                 return 0;

  49.         DWORD i,sqrt_n,p;
  50.         DWORD pn;
  51.         //----------------------------
  52.         sqrt_n=(DWORD)(sqrt(n)+1);
  53.        
  54.         while (sqrt_n * sqrt_n >=n)
  55.                 sqrt_n--;
  56.        
  57.         memset(pBuff,1,n*sizeof(BYTE));
  58.         *pBuff=0;
  59.         *(pBuff+1)=0; //1不是质数
  60.         *(pBuff+2)=1; //2是质数

  61.         for ( i=4;i<n;i+=2) //筛去2的倍数,不包括2
  62.                 *(pBuff+i)=0;
  63.        
  64.         for (i=3;i<=sqrt_n;) //i 是质数
  65.         {
  66.                 for (p=i*i;p<n; p+=2*i)
  67.                         *(pBuff+p)=0; //筛去i 的奇数倍
  68.                 i+=2;
  69.                 while ( *(pBuff+i)==0 )
  70.                         i++; //前进到下一个质数
  71.         }
  72.         pn=0; //素数个数
  73.         for (i=2;i<n;i++)
  74.         {
  75.                 if (pBuff[i])
  76.                 {
  77.                         primeTab[pn++]=i;
  78.                 }
  79.         }
  80.         delete[] pBuff;
  81.         return pn;
  82. }
复制代码
  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <math.h>
  4. #include <memory.h>
  5. #include <windows.h>
  6. #include <assert.h>

  7. typedef unsigned __int64 UINT64;

  8. typedef unsigned long DWORD;
  9. typedef unsigned short WORD;
  10. typedef unsigned char BYTE;


  11. #define SIFT_BLOCK_LEN        (64*1024)

  12. typedef struct fac_struct
  13. {
  14.         DWORD n;
  15.         WORD digSum;
  16. }FAC_ST;

  17. BYTE        g_digSumTab[10000];
  18. FAC_ST        g_facArray[SIFT_BLOCK_LEN+1];

  19. extern UINT64 maxPrimeCount(UINT64 n);
  20. extern DWORD Sift32(DWORD n,DWORD *primeTab);

  21. WORD getDigSum(DWORD n)
  22. {
  23.         DWORD s=0;
  24.         while (n>0)
  25.         {
  26.                 s+= (n % 10);
  27.                 n/=10;
  28.         }
  29.         return (WORD)s;
  30. }

  31. WORD getDigSum2(DWORD n)
  32. {
  33.         assert(g_digSumTab[1]==1);

  34.         DWORD s=0;
  35.         if (n<10000)
  36.                 return g_digSumTab[n];
  37.         else if (n<100000000)
  38.         {
  39.                 DWORD t=n % 10000;
  40.                 n/=10000;
  41.                 return g_digSumTab[n]+g_digSumTab[t];
  42.         }
  43.         else
  44.         {
  45.                 DWORD t0,t1;
  46.                 t0=n % 10000;
  47.                 n/=10000;
  48.                 t1=n % 10000;
  49.                 n/=10000;
  50.                 return g_digSumTab[n]+g_digSumTab[t0]+g_digSumTab[t1];
  51.         }
  52. }

  53. void initTable(BYTE *tab,DWORD n)
  54. {
  55.         DWORD i;
  56.         assert(n==10 || n==100 || n==1000 || n==10000);
  57.         for (i=0;i<n;i++)
  58.         {
  59.                 tab[i]=(BYTE)getDigSum(i);
  60.         }
  61. }


  62. void searchSmithNumber(DWORD n,const char *fileName)
  63. {
  64.         DWORD *primeTab=NULL;
  65.         DWORD base;
  66.         DWORD s_n;
  67.         DWORD primeCount;
  68.         DWORD i,count;
  69.         DWORD time1,time2;
  70.         DWORD c;
  71.         DWORD numBuff[1024];
  72.         int   countInBuff;
  73.         FILE *fp=NULL;

  74.         if (fileName!=NULL)
  75.         {
  76.                 fp=fopen(fileName,"wb");
  77.         }

  78.         time1=GetTickCount();
  79.         initTable(g_digSumTab,10000 );

  80.         s_n =(DWORD)sqrt((double)n)+2;
  81.         while (s_n * s_n >n)
  82.                 s_n--;

  83.         //估算s_n以内素数个数的上限;
  84.         c=(DWORD)maxPrimeCount(s_n);
  85.        
  86.         primeTab=new DWORD[c];
  87.         if (primeTab==NULL )
  88.         {
  89.                 goto thisExit;
  90.         }
  91.        
  92.         primeCount=Sift32(s_n+1,primeTab);
  93.        
  94.         count=0;        base=0;
  95.         countInBuff=0;
  96.         while (base<=n)
  97.         {
  98.                 DWORD end=base+SIFT_BLOCK_LEN-1;
  99.                
  100.                 for (i=0;i<SIFT_BLOCK_LEN;i++)
  101.                 {
  102.                         g_facArray[i].n=base+i;
  103.                         g_facArray[i].digSum=0;
  104.                 }
  105.                
  106.                 i=0;
  107.                 while (primeTab[i]*primeTab[i]<=end && i<primeCount)
  108.                 {
  109.                         DWORD k,m;
  110.                        
  111.                         k= base/primeTab[i];
  112.                         if (k<primeTab[i])
  113.                                 k=primeTab[i];
  114.                         m= k*primeTab[i];
  115.                         if (m<base)
  116.                                 m+= primeTab[i];
  117.                        
  118.                         while (m <=end)
  119.                         {
  120.                                 while ( g_facArray[m-base].n  % primeTab[i] ==0 )
  121.                                 {
  122.                                         g_facArray[m-base].n /= primeTab[i];
  123.                                         g_facArray[m-base].digSum += getDigSum2(primeTab[i]);
  124.                                 }
  125.                                 m+= primeTab[i];
  126.                         }
  127.                         i++;
  128.                 }
  129.                
  130.                 for (i=0;i<SIFT_BLOCK_LEN && i+base<=n;i++)
  131.                 {
  132.                         DWORD x=i+base;
  133.                         WORD s1= getDigSum2(x);
  134.                         WORD s2= g_facArray[i].digSum;
  135.                          
  136.                         if ( g_facArray[i].n !=1 && g_facArray[i].n != x)
  137.                                 s2+=getDigSum2(g_facArray[i].n);
  138.                         if (s1==s2 && s1>0)
  139.                         {
  140.                                 if (fp!=NULL)
  141.                                 {
  142.                                         numBuff[countInBuff++]=x;
  143.                                         if ( countInBuff*sizeof(DWORD)>=sizeof(numBuff) )
  144.                                         {
  145.                                                 fwrite(numBuff,sizeof(DWORD),countInBuff,fp);
  146.                                                 countInBuff=0;
  147.                                         }
  148.                                 }
  149.                                 count++;
  150. #ifdef _DEBUG
  151.                                 printf("%u\n",x);
  152. #endif
  153.                         }
  154.                 }

  155.                 base+= SIFT_BLOCK_LEN;
  156.         }
  157.        
  158.         if (fp!=NULL && countInBuff!=0 )
  159.         {
  160.                 fwrite(numBuff,sizeof(DWORD),countInBuff,fp);
  161.         }
  162.        
  163.         time2=GetTickCount()-time1;
  164.         printf("It take %d ms,total %d smith num within %d\n",time2,count,n);

  165. thisExit:
  166.         if (primeTab!=NULL)
  167.         {
  168.                 delete[] primeTab;        primeTab=NULL;
  169.         }
  170.         if (fp!=NULL)
  171.         {
  172.                 fclose(fp);        fp=NULL;
  173.         }
  174. }
复制代码
  1. #include "stdlib.h"
  2. #include "searchSmithNum.h"

  3. void searchSmithNumber(DWORD n,const char *fileName);
  4. int main(int argc, char* argv[])
  5. {
  6.         DWORD n=100000000;       
  7.         searchSmithNumber(n,NULL); //不输出
  8.         //searchSmithNumber(n,"smith.dat"); //输出到2进制文件
  9.        
  10.         //计算10^8以内的smith数,不输出,在 PM1.7 需 23.766秒,若输出到话,时间将轻微增加,为24.109秒
  11.         //计算10^9以内的smith数,不输出,在 PM1.7 需 252.578秒
  12.         return 0;
  13. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-27 10:21:36 | 显示全部楼层


你这个没达到帖子目的
要求4个以上连续是smith Number的数字的例子
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-27 10:24:51 | 显示全部楼层
该程序在 PIV2.6G 电脑上运行,计算10^8以内的smith数,需 22.128秒。这个CPU
为Northwood 系列,总线 200M,13倍频,512M L2 chche. 需要调整SIFT_BLOCK_LEN    为   (32*1024),可获得最好的性能。
  如果这个程序在主流CPU上运行,速度将更快。大家可给出测试结果。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2021-2-25 04:35 , Processed in 0.065725 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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