题目160：阶乘数字的结尾数字

欧拉计划

Factorial trailing digits

For any N, let f(N) be the last five digits before the trailing zeroes in N!.
For example,

      9! = 362880 so f(9)=36288
    10! = 3628800 so f(10)=36288
    20! = 2432902008176640000 so f(20)=17664

Find f(1,000,000,000,000)

题目：

对任意的 N，令定义 f(N) 为 N! ( N 的阶乘）除末尾的 0 以外的最后 5 个数字，比如

      9! = 362880  所以  f(9)=36288
    10! = 3628800  所以  f(10)=36288
    20! = 2432902008176640000  所以  f(20)=17664

求 f(1,000,000,000,000）

guosl

我们用CUDA进行死算。
答案：16576

1. /*
   
2. 答案：16576
   
3. 耗时：184.763907s(cmp30hx，最近矿卡非常便宜 )
   
4. */
   
5. #include <cstdio>
   
6. #include <cuda_runtime.h>
   
7. #include <device_launch_parameters.h>
   
8. #include <helper_cuda.h>
   
9. #include <omp.h>
   
10. using namespace std;
    
11. 
    
12. const long long N = 1000000000000LL;
    
13. const int nStep = 46080 * 16; //分段的区间长度
    
14. const int nBlocksSize = 128;
    
15. const int nBlocks = (nStep + nBlocksSize - 1) / nBlocksSize;
    
16. int* result_d = NULL;//存放结果的显卡内存
    
17. int* result = NULL;//存放结果的主机内存
    
18. 
    
19. int qpow(long long k)
    
20. {
    
21.   long long p;
    
22.   if (k == 0)
    
23.     return 1;
    
24.   if (k == 1)
    
25.     return 2;
    
26.   if ((k & 1) == 0)
    
27.   {
    
28.     p = qpow(k >> 1);
    
29.     p *= p;
    
30.     p %= 100000;
    
31.   }
    
32.   else
    
33.   {
    
34.     p = qpow(k - 1);
    
35.     p *= 2;
    
36.     p %= 100000;
    
37.   }
    
38.   return (int)p;
    
39. }
    
40. 
    
41. __global__ void getpr(long long b, int *result)
    
42. {
    
43.   extern __shared__ long long nVal[];
    
44.   nVal[threadIdx.x] = 1;
    
45.   long long k = blockIdx.x * nBlocksSize + threadIdx.x + b;//当前的计算序号
    
46.   if (k <= N)
    
47.   {
    
48.     long long p = k;
    
49.     while ((p & 1) == 0)
    
50.       p = p >> 1;
    
51.     while (p % 5 == 0)
    
52.       p /= 5;
    
53.     p %= 100000;
    
54.     nVal[threadIdx.x] = p;
    
55.   }
    
56.   __syncthreads();
    
57.   //累乘本block内得到的值
    
58.   for (int k = nBlocksSize >> 1; k > 0; k = k >> 1)
    
59.   {
    
60.     if (threadIdx.x < k)
    
61.     {
    
62.       nVal[threadIdx.x] *= nVal[threadIdx.x + k];
    
63.       nVal[threadIdx.x] %= 100000;
    
64.     }
    
65.     __syncthreads();
    
66.   }
    
67.   if (threadIdx.x == 0)
    
68.     result[blockIdx.x] = (int)nVal[threadIdx.x];//放入返回的变量中
    
69. }
    
70. 
    
71. int main(void)
    
72. {
    
73.   long long k = 1, ppr = 1;
    
74.   long long c2 = 0, c5 = 0, n = N;
    
75.   double t = omp_get_wtime();
    
76.   cudaDeviceProp deviceProp;
    
77.   int devID = 0;
    
78.   cudaError_t cudaStatus = cudaSetDevice(devID);//设置计算卡
    
79.   if (cudaStatus != cudaSuccess)
    
80.     goto end;
    
81.   cudaStatus = cudaGetDeviceProperties(&deviceProp, devID);
    
82.   if (cudaStatus != cudaSuccess)
    
83.     goto end;
    
84.   printf("GPU Device %d: "%s" with compute capability %d.%d\n\n", devID, deviceProp.name, deviceProp.major, deviceProp.minor);
    
85.   while (n > 1)
    
86.   {
    
87.     c2 += n / 2;
    
88.     n = n >> 1;
    
89.   }
    
90.   n = N;
    
91.   while (n > 1)
    
92.   {
    
93.     c5 += n / 5;
    
94.     n /= 5;
    
95.   }
    
96.   c2 -= c5;
    
97.   ppr = qpow(c2);
    
98.   result = new int[nBlocks];
    
99.   if (result == NULL)
    
100.     goto end;
     
101.   cudaStatus = cudaMalloc((void**)&result_d, nBlocks * sizeof(int));
     
102.   if (cudaStatus != cudaSuccess)
     
103.     goto end;
     
104.   while (k <= N)
     
105.   {
     
106.     cudaStatus = cudaMemset(result_d, 0, nBlocks * sizeof(int));//将返回的数组清零
     
107.     if (cudaStatus != cudaSuccess)
     
108.       goto end;
     
109.     getpr << <nBlocks, nBlocksSize, nBlocksSize * sizeof(long long) >> > (k, result_d);
     
110.     cudaError_t cudaStatus = cudaGetLastError();
     
111.     if (cudaStatus != cudaSuccess)
     
112.       goto end;
     
113.     cudaDeviceSynchronize();//等待计算完成
     
114.                             //将区间和复制回主机内存
     
115.     cudaStatus = cudaMemcpy(result, result_d, nBlocks * sizeof(int), cudaMemcpyDeviceToHost);
     
116.     if (cudaStatus != cudaSuccess)
     
117.       goto end;
     
118.     for (int i = 0; i < nBlocks; ++i)
     
119.     {
     
120.       ppr *= (long long)result[i];
     
121.       ppr %= 100000;
     
122.     }
     
123.     k += nStep;
     
124.   }
     
125.   cudaDeviceReset();
     
126.   t = omp_get_wtime() - t;
     
127.   printf("%lld\n%lf\n", ppr, t);
     
128. end:
     
129.   if (result != NULL)
     
130.     delete[] result;
     
131.   if (result_d != NULL)
     
132.     cudaFree(result_d);
     
133.   return 0;
     
134. }

复制代码

guosl

我们再给一种解法：
当n%5==0时有公式：n!=n-!*5^(n/5)*(n/5)!这个可以应用数学归纳法来证明。
其中n-!为普通阶乘中去掉所有5的倍数，例如：10-！= 1*2*3*4*6*7*8*9
由于10^(n/5)只产生一串零，故在计算中舍去。只需计算：n-!/(2^(n/5))*(n/5)!
对于n-!/(2^(n/5))的计算
我们把n个数分成若干段，每段有如下形式
(kM+0)*(kM+1)*...*(kM+4)*(kM+6)*...*(kM+9)*(kM+11)*...*(kM+i)*...*(kM+M-1)
这样的段有n/M段，而每段在%M的作用下都变成了
(1*...*4*6*...*9*11*...*i*...*(M-1))
整个计算变成了
(1*...*4*6*...*9*11*...*i*...*(M-1))^(n/M)
而根据扩展欧拉定理：i^(n/M)=i^(Phi(M)+(n/M)%Phi(M))，其中Phi()为欧拉函数
所以：f(n)=(g(n)*f(n/5))%100000
我们可以根据这个得到下面的递归计算
答案：16576，计算耗时可以忽略不记。

1. #include <iostream>
   
2. using namespace std;
   
3. 
   
4. const long long N = 1000000000000LL;
   
5. const int M = 100000, Phi_M = 40000;//Phi_M为M的欧拉函数值
   
6. 
   
7. int qpow(int b, int n)//快速幂
   
8. {
   
9.   long long x = b;
   
10.   if (n == 0)
    
11.     return 1;
    
12.   if (n == 1)
    
13.     return int(x % M);
    
14.   if ((n & 1) == 0)
    
15.   {
    
16.     x = qpow(b, n >> 1);
    
17.     x *= x;
    
18.     return int(x % M);
    
19.   }
    
20.   else
    
21.   {
    
22.     x = qpow(b, n - 1);
    
23.     x *= b;
    
24.     return int(x % M);
    
25.   }
    
26. }
    
27. 
    
28. int g(long long n)//计算n-!/2^(n/5)的值
    
29. {
    
30.   //下面计算的依据是扩展欧拉定理
    
31.   long long m1 = n / M;
    
32.   if (m1 > Phi_M)
    
33.     m1 = m1 % Phi_M + Phi_M;
    
34.   long long m2 = (m1 + 1) / 2;
    
35.   if (m2 > Phi_M)
    
36.     m2 = m2 % Phi_M + Phi_M;
    
37.   long long m3 = m1 / 2;
    
38.   if (m3 > Phi_M)
    
39.     m3 = m3 % Phi_M + Phi_M;
    
40.   long long r = 1;
    
41.   for (int i = 0; i < M; ++i)
    
42.   {
    
43.     if (i % 5 == 0)
    
44.       continue;
    
45.     //我们将2^(-n/5)分配到乘积中，具体就是个位数是2或4的数除个2
    
46.     if (i % 10 == 2)
    
47.     {
    
48.       r *= qpow(int(i / 2), (int)m2);//处理其中除2的因子，这个要具体计算
    
49.       r %= M;
    
50.       r *= qpow(int((M + i) / 2), (int)m3);
    
51.     }
    
52.     else if (i % 10 == 4)
    
53.     {
    
54.       r *= qpow(int(i / 2), (int)m2);//处理其中除2的因子
    
55.       r %= M;
    
56.       r *= qpow(int((M + i) / 2), (int)m3);
    
57.     }
    
58.     else
    
59.       r *= qpow(int(i), (int)m1);
    
60.     r %= M;
    
61.   }
    
62.   return (int)r;
    
63. }
    
64. 
    
65. int f(long long n)
    
66. {
    
67.   long long r = 1;
    
68.   if (n % M == 0)
    
69.   {
    
70.     r = g(n);
    
71.     r *= f(n / 5);//递归
    
72.     r %= M;
    
73.     return (int)r;
    
74.   }
    
75.   //n的值已经减小到可以直接计算了
    
76.   int c2 = 0;
    
77.   {
    
78.     int k = (int)n, c5 = 0;
    
79.     while (k >= 2)
    
80.     {
    
81.       c2 += k / 2;
    
82.       k /= 2;
    
83.     }
    
84.     k = (int)n;
    
85.     while (k >= 5)
    
86.     {
    
87.       c5 += k / 5;
    
88.       k /= 5;
    
89.     }
    
90.     c2 -= c5;
    
91.   }
    
92.   r = qpow(2, c2);
    
93.   for (int i = 1; i <= n; ++i)
    
94.   {
    
95.     int k = i;
    
96.     while ((k & 1) == 0)
    
97.       k /= 2;
    
98.     while (k % 5 == 0)
    
99.       k /= 5;
    
100.     k %= M;
     
101.     r *= k;
     
102.     r %= M;
     
103.   }
     
104.   return (int)r;
     
105. }
     
106. 
     
107. int main(void)
     
108. {
     
109.   cout << f(N) << endl;//调用递归计算
     
110.   return 0;
     
111. }

复制代码