K-Smith numbers

无心人

ksmith100000000.TXT (46.03 KB, 下载次数: 3)

2008-11-25 18:59 上传

点击文件名下载附件

无心人

1小时14分钟

liangbch

看来无心人你这个程序还是难敌C程序呀。到今日凌晨5点，我也完成了我的C程序。这个程序的复杂度近似于O(n)，在PM 1.7上,计算10^8以内的smith数(不输出)，不到10秒，好像是4.2秒，可是U盘忘带了。等明天吧，明天我贴出程序。

无心人

不是不敌 是我初学 写不出高效的算法啊 如果C用的时间是1 Haskell通常在2-5

无心人

用筛法做，是最佳的选择了，呵呵 不过啊 haskell实现筛法，似乎有点挑战啊 不是haskell不行 是俺笨哦

liangbch

使用晒法，我在5楼就提到了，关于复杂度我也在5楼提到了。只是在实现的时候发现不需使用内存池和链表。 另外，可以采用一些优化手段。我的做法是事先算出10000以下所有的数的数字和，用的时候查表。

无心人

呵呵 这里我并不奢望多高等的算法 当学习haskell了 目前的估计是 30分钟 用的是逐个分解的算法 另外，算100000000个数字的数字和 是5秒

liangbch

贴出我的解法： 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"
6. typedef unsigned __int64 UINT64;
7. typedef unsigned long DWORD;
8. typedef unsigned short WORD;
9. typedef unsigned char BYTE;
12. //文件描述：使用最简单的方法筛素数
14. // 一定范围内质数的个数
15. // 10 4
16. // 100 25
17. // 1,000 168
18. // 10,000 1,229
19. // 100,000 9,592
20. // 1,000,000 78,498
21. // 10,000,000 664,579
22. // 100,000,000 5,761,455
23. // 1,000,000,000 50,847,534
24. // 10,000,000,000 455,052,511
25. // 100,000,000,000 4,118,054,813
26. // 1,000,000,000,000 37,607,912,018
27. // 10,000,000,000,000 346,065,536,839
28. // 100,000,000,000,000 3,204,941,750,802
29. // 1,000,000,000,000,000 29,844,570,422,669
30. // 10,000,000,000,000,000 279,238,341,033,925
31. // 100,000,000,000,000,000 2,623,557,157,654,233
32. // 1,000,000,000,000,000,000 24,739,954,287,740,860
33. // 估计x 范围内的质数的个数， 欧拉公式 pi(x)=x /( ln(x));
35. // 经过测试，
36. // 在1-1百万之间，相邻2个质数的最大差为114,这两个质数分别为492113,492227
37. // 在1-1千万之间，相邻2个质数的最大差为154,这两个质数分别为4652353,4652507
38. // 在1-3千万之间，相邻2个质数的最大差为210,这两个质数分别为20831323,20831533
39. // 在1-1亿之间， 相邻2个质数的最大差为220,这两个质数分别为47326693,47326913
42. //计算n以内素数个数的上限
43. UINT64 maxPrimeCount(UINT64 n)
44. {
45. double f= (double)n / log((double)n) * 10.0/9.0 +32;
46. return (UINT64)f;
47. }
49. DWORD Sift32(DWORD n,DWORD *primeTab)
50. // 对0到n 之间(不包括n)的数进行筛选，返回质数的个数
51. // n 不能大于L2 cache，否则效能很低
52. //**********************************************
53. {
54. BYTE *pBuff= new BYTE[n];
55. if (pBuff==NULL)
56. return 0;
58. DWORD i,sqrt_n,p;
59. DWORD pn;
60. //----------------------------
61. sqrt_n=(DWORD)(sqrt(n)+1);
63. while (sqrt_n * sqrt_n >=n)
64. sqrt_n--;
66. memset(pBuff,1,n*sizeof(BYTE));
67. *pBuff=0;
68. *(pBuff+1)=0; //1不是质数
69. *(pBuff+2)=1; //2是质数
71. for ( i=4;i<n;i+=2) //筛去2的倍数，不包括2
72. *(pBuff+i)=0;
74. for (i=3;i<=sqrt_n;) //i 是质数
75. {
76. for (p=i*i;p<n; p+=2*i)
77. *(pBuff+p)=0; //筛去i 的奇数倍
78. i+=2;
79. while ( *(pBuff+i)==0 )
80. i++; //前进到下一个质数
81. }
82. pn=0; //素数个数
83. for (i=2;i<n;i++)
84. {
85. if (pBuff[i])
86. {
87. primeTab[pn++]=i;
88. }
89. }
90. delete[] pBuff;
91. return pn;
92. }

复制代码

1. #include <stdlib.h>
2. #include <stdio.h>
3. #include <math.h>
4. #include <memory.h>
5. #include <windows.h>
6. #include <assert.h>
8. typedef unsigned __int64 UINT64;
10. typedef unsigned long DWORD;
11. typedef unsigned short WORD;
12. typedef unsigned char BYTE;
15. #define SIFT_BLOCK_LEN (64*1024)
17. typedef struct fac_struct
18. {
19. DWORD n;
20. WORD digSum;
21. }FAC_ST;
23. BYTE g_digSumTab[10000];
24. FAC_ST g_facArray[SIFT_BLOCK_LEN+1];
26. extern UINT64 maxPrimeCount(UINT64 n);
27. extern DWORD Sift32(DWORD n,DWORD *primeTab);
29. WORD getDigSum(DWORD n)
30. {
31. DWORD s=0;
32. while (n>0)
33. {
34. s+= (n % 10);
35. n/=10;
36. }
37. return (WORD)s;
38. }
40. WORD getDigSum2(DWORD n)
41. {
42. assert(g_digSumTab[1]==1);
44. DWORD s=0;
45. if (n<10000)
46. return g_digSumTab[n];
47. else if (n<100000000)
48. {
49. DWORD t=n % 10000;
50. n/=10000;
51. return g_digSumTab[n]+g_digSumTab[t];
52. }
53. else
54. {
55. DWORD t0,t1;
56. t0=n % 10000;
57. n/=10000;
58. t1=n % 10000;
59. n/=10000;
60. return g_digSumTab[n]+g_digSumTab[t0]+g_digSumTab[t1];
61. }
62. }
64. void initTable(BYTE *tab,DWORD n)
65. {
66. DWORD i;
67. assert(n==10 || n==100 || n==1000 || n==10000);
68. for (i=0;i<n;i++)
69. {
70. tab[i]=(BYTE)getDigSum(i);
71. }
72. }
75. void searchSmithNumber(DWORD n,const char *fileName)
76. {
77. DWORD *primeTab=NULL;
78. DWORD base;
79. DWORD s_n;
80. DWORD primeCount;
81. DWORD i,count;
82. DWORD time1,time2;
83. DWORD c;
84. DWORD numBuff[1024];
85. int countInBuff;
86. FILE *fp=NULL;
88. if (fileName!=NULL)
89. {
90. fp=fopen(fileName,"wb");
91. }
93. time1=GetTickCount();
94. initTable(g_digSumTab,10000 );
96. s_n =(DWORD)sqrt((double)n)+2;
97. while (s_n * s_n >n)
98. s_n--;
100. //估算s_n以内素数个数的上限;
101. c=(DWORD)maxPrimeCount(s_n);
103. primeTab=new DWORD[c];
104. if (primeTab==NULL )
105. {
106. goto thisExit;
107. }
109. primeCount=Sift32(s_n+1,primeTab);
111. count=0; base=0;
112. countInBuff=0;
113. while (base<=n)
114. {
115. DWORD end=base+SIFT_BLOCK_LEN-1;
117. for (i=0;i<SIFT_BLOCK_LEN;i++)
118. {
119. g_facArray[i].n=base+i;
120. g_facArray[i].digSum=0;
121. }
123. i=0;
124. while (primeTab[i]*primeTab[i]<=end && i<primeCount)
125. {
126. DWORD k,m;
128. k= base/primeTab[i];
129. if (k<primeTab[i])
130. k=primeTab[i];
131. m= k*primeTab[i];
132. if (m<base)
133. m+= primeTab[i];
135. while (m <=end)
136. {
137. while ( g_facArray[m-base].n % primeTab[i] ==0 )
138. {
139. g_facArray[m-base].n /= primeTab[i];
140. g_facArray[m-base].digSum += getDigSum2(primeTab[i]);
141. }
142. m+= primeTab[i];
143. }
144. i++;
145. }
147. for (i=0;i<SIFT_BLOCK_LEN && i+base<=n;i++)
148. {
149. DWORD x=i+base;
150. WORD s1= getDigSum2(x);
151. WORD s2= g_facArray[i].digSum;
153. if ( g_facArray[i].n !=1 && g_facArray[i].n != x)
154. s2+=getDigSum2(g_facArray[i].n);
155. if (s1==s2 && s1>0)
156. {
157. if (fp!=NULL)
158. {
159. numBuff[countInBuff++]=x;
160. if ( countInBuff*sizeof(DWORD)>=sizeof(numBuff) )
161. {
162. fwrite(numBuff,sizeof(DWORD),countInBuff,fp);
163. countInBuff=0;
164. }
165. }
166. count++;
167. #ifdef _DEBUG
168. printf("%u\n",x);
169. #endif
170. }
171. }
173. base+= SIFT_BLOCK_LEN;
174. }
176. if (fp!=NULL && countInBuff!=0 )
177. {
178. fwrite(numBuff,sizeof(DWORD),countInBuff,fp);
179. }
181. time2=GetTickCount()-time1;
182. printf("It take %d ms,total %d smith num within %d\n",time2,count,n);
184. thisExit:
185. if (primeTab!=NULL)
186. {
187. delete[] primeTab; primeTab=NULL;
188. }
189. if (fp!=NULL)
190. {
191. fclose(fp); fp=NULL;
192. }
193. }

复制代码

1. #include "stdlib.h"
2. #include "searchSmithNum.h"
4. void searchSmithNumber(DWORD n,const char *fileName);
5. int main(int argc, char* argv[])
6. {
7. DWORD n=100000000;
8. searchSmithNumber(n,NULL); //不输出
9. //searchSmithNumber(n,"smith.dat"); //输出到2进制文件
11. //计算10^8以内的smith数,不输出，在 PM1.7 需 23.766秒,若输出到话，时间将轻微增加，为24.109秒
12. //计算10^9以内的smith数,不输出，在 PM1.7 需 252.578秒
13. return 0;
14. }

复制代码

无心人

你这个没达到帖子目的 要求4个以上连续是smith Number的数字的例子

liangbch

该程序在 PIV2.6G 电脑上运行，计算10^8以内的smith数，需 22.128秒。这个CPU 为Northwood 系列，总线 200M，13倍频，512M L2 chche. 需要调整SIFT_BLOCK_LEN 为 (32*1024)，可获得最好的性能。 如果这个程序在主流CPU上运行，速度将更快。大家可给出测试结果。