Ten digit numbers

mathe

上面的公式有个问题，没有考虑到对于越大的数，对应的计算量越大。由于我们实际计算过程都是一个超大整数乘上一个小于\(10^{10}\)的整数
所以实际上每次乘法本身就是O(n)的，由此总体复杂度可以表示为
\(\sum_{x=10^9}^{10^{10}-1}\frac{x\ln(0.1)}{\ln(10^{-10}x)}\approx 10^{20}\ln(0.1)(Ei(2\ln(1-10^{-10}))-Ei(2\ln(0.1)))\).
而计算过程中，如果我们从\(10^9\)已经计算到\(10^{10}-10^7\),那么我们就可以估算出计算过程已经占了整体计算量的
\(\frac{Ei(2\ln(0.999))-Ei(2\ln(0.1)))}{Ei(2\ln(1-10^{-10}))-Ei(2\ln(0.1)))} \approx 0.2088\)
通过这种方法可以估算出完成整个计算大概需要8K CPU日左右, 对于现代多核处理器也不是不能完成，就是时间略微长了点，如果能够优化一下会更好。

1. #include <math.h>
   
2. #include <string.h>
   
3. #include <stdio.h>
   
4. #include <time.h>
   
5. 
   
6. #define SN 501
   
7. #define MP 100000000
   
8. 
   
9. #define MOD 10000000000L
   
10. #define HMOD 100000
    
11. 
    
12. char dc[HMOD][10];
    
13. void initdc()
    
14. {
    
15.   int i;
    
16.   for(i=0;i<HMOD;i++){
    
17.     int j;
    
18.     int x=i;
    
19.     for(j=0;j<5;j++){
    
20.       dc[i][x%10]++;
    
21.       x/=10;
    
22.     }
    
23.   }
    
24. }
    
25. 
    
26. void get_count(long v[],int length, int count[10])
    
27. {
    
28.   int i;
    
29.   for(i=0;i<10;i++)count[i]=0;
    
30.   for(i=0;i<length;i++){
    
31.     int x=v[i]%HMOD;
    
32.     int y=v[i]/HMOD;
    
33.     int j;
    
34.     for(j=0;j<10;j++){
    
35.       count[j]+=dc[x][j]+dc[y][j];
    
36.     }
    
37.   }
    
38. }
    
39. 
    
40. bool mulby(long v[], int length, long m)
    
41. {
    
42.   int i;
    
43.   long carry=0L;
    
44.   for(i=0;i<length;i++){
    
45.     __int128_t prod = (__int128_t)v[i]*m+carry;
    
46.     v[i]=prod%MOD;
    
47.     carry=prod/MOD;
    
48.   }
    
49.   v[length]=carry;
    
50.   return carry>=(MOD/10);
    
51. }
    
52. 
    
53. int main()
    
54. {
    
55.     double sv=pow(0.1,1.0/SN)*MOD;
    
56.     long is = (long)ceil(sv);
    
57.     initdc();
    
58.     time_t start=time(NULL);
    
59. #pragma omp parallel
    
60.     {
    
61.       long *data = (long *)malloc((MP+1)*sizeof(long));
    
62. #pragma omp for schedule(dynamic)
    
63.     for(long i=is;i<MOD;i++){
    
64.             int count[10];
    
65.       data[0]=i;
    
66.         int j;
    
67.         for(j=2;j<=MP;j++){
    
68.     if(!mulby(data,j-1,i))break;
    
69.     if(j<SN)continue;
    
70.     get_count(data,j,count);
    
71.             int k;
    
72.             for(k=0;k<10;k++)if(count[k]!=j)break;
    
73.             if(k==10){
    
74. #pragma omp critical
    
75.                     {
    
76.                         printf("%ld^%d\n",i,j);
    
77.                               fflush(stdout);
    
78.                     }
    
79.             }
    
80.         }
    
81. #pragma omp critical
    
82.         if(i%1000000==0)
    
83.         {
    
84.                 time_t now=time(NULL);
    
85.                 fprintf(stderr,"%ld cost %ld\n",i,now-start);
    
86.                 fflush(stderr);
    
87.         }
    
88.     }
    
89.     }
    
90.     printf("done!\n");
    
91.     return 0;
    
92. }
    

复制代码