回归数

medie2005

设k位数n的十进制形式为：$bar{d_{1}d_{2}...d_{k}}$。 若n满足： $bar{d_{1}d_{2}...d_{k}}=\sum_{i=1}^{k}{d_{i}^k}$ 则称n为回归数。 比如，153就是一个回归数：$153=1^{3}+5^{3}+3^{3}$。 最大的回归数有39位，40位以上（含40位）就不存在回归数了。 请写一个代码求出所有回归数，看谁的用时最少。 [所有回归数是已知的，我只是想看看是否有很好的方法来求出这些数]

无心人

似乎是我提出问题的子集吧？？

无心人

我猜测最小时间是1分钟内

gxqcn

环比不动点麻烦得多。 但在1分钟内还是蛮有挑战性的。

medie2005

1分钟内?太夸张了吧.我倒是很希望看到这样的算法.

无心人

我算17的环是10秒

medie2005

我现在的算法求20位以下的回归数大概要9秒, 求出所有的回归数要1个小时.

无心人

抱歉，我估计的确实有问题 10分钟都不可能算出来 程序存在问题 修改后再说

winxos

应该可以先进行下数学处理吧？ 比如有那些数是不可能的。

medie2005

发一下代码:

1. #include <cstdlib>
2. #include <iostream>
3. #include <windows.h>
4. #include <cmath>
5. #include "gmp.h"
7. using namespace std;
9. const int L=40;
10. mpz_t Hash[10][L];
11. mpz_t dif[10][10][L];
12. mpz_t POW10[L];
13. char str[L<<1];
15. void Init( ){
16. unsigned i, j;
17. for( i=0; i<10; ++i ){
18. for( j=0; j<L; ++j )
19. mpz_init( Hash[i][j] );
20. }
21. mpz_set_ui( Hash[0][0], 0 );
22. mpz_init_set_ui( POW10[0], 1 );
23. for( i=1; i<10; ++i )
24. mpz_set_ui( Hash[i][0], 1 );
25. for( j=1; j<L; ++j ){
26. mpz_init( POW10[j] );
27. mpz_mul_ui( POW10[j], POW10[j-1], 10 );
28. }
29. for( i=0; i<10; ++i ){
30. for( j=1; j<L; ++j ){
31. mpz_mul_ui( Hash[i][j], Hash[i][j-1], i );
32. }
33. }
34. for( int k=0; k<L; ++k ){
35. for( i=0; i<10; ++i ){
36. for( j=0; j<10; ++j ){
37. mpz_init( dif[i][j][k] );
38. mpz_sub( dif[i][j][k], Hash[i][k], Hash[j][k] );
39. }
40. }
41. }
42. }
44. void Clean( ){
45. unsigned i, j;
46. for( i=0; i<10; ++i )
47. for( j=0; j<L; ++j ){
48. mpz_clear( Hash[i][j] );
49. }
50. for( j=0; j<L; ++j )
51. mpz_clear( POW10[j] );
52. for( int k=0; k<L; ++k ){
53. for( i=0; i<10; ++i ){
54. for( j=0; j<10; ++j ){
55. mpz_clear( dif[i][j][k] );
56. }
57. }
58. }
59. }
61. bool check( mpz_t &n, int *c ){
62. int i, d[10]={0};
63. mpz_get_str( str, 10, n );
64. for( i=0; str[i]!='\0'; ++i ){
65. d[str[i]-'0']++;
66. if( d[str[i]-'0']>c[str[i]-'0'] ) return false;
67. }
68. return true;
69. }
71. void combination( unsigned n, unsigned m ){
72. int *c=new int[m+2], j, *num=new int[m+2];
73. int dig[10]={0};
74. dig[0]=m;
75. for( j=1; j<=m; ++j ) c[j]=j-1, num[j]=c[j]-j+1;
76. c[m+1]=n;
77. mpz_t number;
78. mpz_init_set_ui( number, 0 );
79. R2:
80. if( mpz_cmp( number, POW10[m-1] )>0 && mpz_cmp( number, POW10[m] )<0 && check( number, dig ) ){
81. mpz_out_str( stdout, 10, number );
82. printf("\n");
83. }
84. if( m&1 )
85. if( c[1]+1<c[2] ){
86. --dig[num[1]];
87. mpz_add( number, number, dif[num[1]+1][num[1]][m] );
88. ++c[1], ++num[1];
89. ++dig[num[1]];
90. goto R2;
91. }
92. else{
93. j=2; goto R4;
94. }
95. else
96. if( c[1]>0 ){
97. --dig[num[1]];
98. mpz_sub( number, number, dif[num[1]][num[1]-1][m] );
99. --c[1], --num[1];
100. ++dig[num[1]];
101. goto R2;
102. }
103. else{
104. j=2; goto R5;
105. };
106. R4: if( c[j]>=j ){
107. if( j<=m ){
108. --dig[num[j]];
109. mpz_sub( number, number, dif[num[j]][c[j-1]-j+1][m] );
110. }
111. if( j-1<=m ){
112. --dig[num[j-1]];
113. mpz_sub( number, number, dif[num[j-1]][0][m] );
114. }
115. c[j]=c[j-1], num[j]=c[j]-j+1; c[j-1]=j-2, num[j-1]=c[j-1]-(j-1)+1;
116. if( j<=m )
117. ++dig[num[j]];
118. if( j-1<=m )
119. ++dig[num[j-1]];
120. goto R2;
121. }
122. else
123. ++j;
124. R5: if( c[j]+1<c[j+1] ){
125. if( j<=m ){
126. --dig[num[j]];
127. mpz_sub( number, number, dif[num[j]][c[j]+1-j+1][m] );
128. }
129. if( j-1<=m ){
130. --dig[num[j-1]];
131. mpz_sub( number, number, dif[num[j-1]][c[j]-j+2][m] );
132. }
133. c[j-1]=c[j], num[j-1]=c[j-1]-(j-1)+1; c[j]=c[j]+1, num[j]=c[j]-j+1;
136. if( j<=m )
137. ++dig[num[j]];
138. if( j-1<=m )
139. ++dig[num[j-1]];
140. goto R2;
141. }
142. else{
143. ++j;
144. if( j<=m ) goto R4;
145. }
146. delete []c;
147. delete []num;
148. mpz_clear( number );
149. }
151. int main(int argc, char *argv[])
152. {
153. int time=GetTickCount();
154. Init();
155. for( unsigned i=3; i<L; ++i ){
156. cout<<i<<" digits search start!\n";
157. combination( i+9, i );
158. }
159. Clean( );
160. cout<<"time : "<<GetTickCount()-time<<" ms\n";
161. cout<<"press any key to continue...\n";
162. cin.get();
163. return EXIT_SUCCESS;
164. }

复制代码