题目154：研究帕斯卡金字塔

欧拉计划

Exploring Pascal's pyramid

A triangular pyramid is constructed using spherical balls so that each ball rests on exactly three balls of the next lower level.

Then, we calculate the number of paths leading from the apex to each position:

A path starts at the apex and progresses downwards to any of the three spheres directly below the current position.

Consequently, the number of paths to reach a certain position is the sum of the numbers immediately above it (depending on the position, there are up to three numbers above it).

The result is Pascal's pyramid and the numbers at each level n are the coefficients of the trinomial expansion (x + y + z)ⁿ.

How many coefficients in the expansion of (x + y + z)²⁰⁰⁰⁰⁰ are multiples of 10¹²?

题目：

三角状的金字塔是用圆球建造的：每个球都与下一层的三个球相接触。

于是，我们计算从顶点到任一点的路径数：

路径是从顶点开始，从下一层相接触的三个球中任选一个一直往下走。

那么，到达某一位置的路径的条数就是到达它上一层的路径数之和（取决于具体位置，一个球，最多可能有三条路径可以从上层到该位置，所以，最多只有三个路径数）。

这个结果就是帕斯卡金字塔，第n层的路径数是三项式乘方 (x + y + z)ⁿ的系数.

请问，(x + y + z)²⁰⁰⁰⁰⁰ 的系数中，有多少个是 10¹²倍数？

gunjang

#（x + y + z)^n = ∑(n!/(r!* s!* t!) * x^r * y^s * z^t) ,r+s+t=n
所以问题变成 n!/(r!* s!* t!) 里有多少10^12的倍数,r+s+t=n
10^12分解成因子就是12个2和12个5
问题变成n!/r!s!t!有几个的2和5的因子数都大于等于12
算法很简单
只所以采用pascal实在是因为python计算太慢了（也可能我的算法不够优化）
delphi7运行的话大概1分半，python目测起码10个小时。。。

1. program p154;
   
2. 
   
3. {$APPTYPE CONSOLE}
   
4. 
   
5. uses
   
6.   SysUtils, Types, Windows,Math;
   
7. const
   
8.   maxn = 200000;
   
9. type
   
10.   nlist = array[0..maxn+5] of integer;
    
11. var
    
12.   st: DWord;
    
13.   twolist, fivelist: nlist;
    
14.   cnt: integer;
    
15.   r,s,t: integer;
    
16. 
    
17.   function getXcountin1n(const x: word): nlist;
    
18.   var
    
19.     cnt: integer;
    
20.     i,i1: integer;
    
21.   begin
    
22.     cnt := 0;
    
23.     for i:=0 to x-1 do
    
24.       result[i] := 0;
    
25.     i := x;
    
26.     while i <= maxn do
    
27.     begin
    
28.       i1 := i;
    
29.       while (i1 > 0) and (i1 mod x = 0) do
    
30.       begin
    
31.         inc(cnt);
    
32.         i1 := i1 div x
    
33.       end;
    
34.       for i1:= i to i+x-1 do
    
35.         result[i1] := cnt;
    
36.       inc(i, x);
    
37.     end;
    
38.   end;
    
39. 
    
40. begin
    
41. st := gettickcount();
    
42. twolist := getXcountin1n(2);
    
43. fivelist := getXcountin1n(5);
    
44. cnt := 0;
    
45. for r:= 0 to maxn do
    
46. begin
    
47.   for s:= 0 to (maxn-r) do
    
48.   begin
    
49.     t := maxn-r-s;
    
50.     if (twolist[maxn]-twolist[r]-twolist[s]-twolist[t] >= 12) and (fivelist[maxn]-fivelist[r]-fivelist[s]-fivelist[t] >= 12) then
    
51.        inc(cnt);
    
52.   end;
    
53. end;
    
54. 
    
55. writeln(cnt);  //479742450 87.4s
    
56. writeln(format('cost %n s', [(gettickcount-st)/1000]));
    
57. readln; //wait Enter
    
58. end.

复制代码

guosl

思路参考上面的帖子

1. /*
   
2. 答案：479742450
   
3. 耗时：1.24439秒(R5 5600G)
   
4. */
   
5. #include <iostream>
   
6. #include <omp.h>
   
7. using namespace std;
   
8. 
   
9. const int N = 200000;
   
10. int p2[N + 1], p5[N + 1];
    
11. 
    
12. int main(void)
    
13. {
    
14.   long long nCount = 0;
    
15.   #pragma omp parallel shared (p2,p5) reduction(+:nCount)
    
16.   {
    
17. #pragma omp for
    
18.     for (int i = 2; i <= N; ++i)
    
19.     {
    
20.       int k = 2;
    
21.       while (k <= i)
    
22.       {
    
23.         p2[i] += i / k;
    
24.         k *= 2;
    
25.       }
    
26.       k = 5;
    
27.       while (k <= i)
    
28.       {
    
29.         p5[i] += i / k;
    
30.         k *= 5;
    
31.       }
    
32.     }
    
33.     for (int i = 0; i <= N; ++i)
    
34.     {
    
35. #pragma omp for
    
36.       for (int j = 0; j <= N - i; ++j)
    
37.       {
    
38.         int k = N - i - j;
    
39.         int c2 = p2[N] - p2[i] - p2[j] - p2[k];
    
40.         int c5 = p5[N] - p5[i] - p5[j] - p5[k];
    
41.         if (c2 >= 12 && c5 >= 12)
    
42.           ++nCount;
    
43.       }
    
44.     }
    
45.   }
    
46.   cout << nCount << endl;
    
47.   return 0;
    
48. }

复制代码