题目203：squarefree的二项式系数

欧拉计划

Squarefree Binomial Coefficients

The binomial coefficients ⁿC_k can be arranged in triangular form, Pascal's triangle, like this:

It can be seen that the first eight rows of Pascal's triangle contain twelve distinct numbers: 1, 2, 3, 4, 5, 6, 7, 10, 15, 20, 21 and 35.

A positive integer n is called squarefree if no square of a prime divides n. Of the twelve distinct numbers in the first eight rows of Pascal's triangle, all except 4 and 20 are squarefree. The sum of the distinct squarefree numbers in the first eight rows is 105.

Find the sum of the distinct squarefree numbers in the first 51 rows of Pascal's triangle.

题目：

二项式系数 ⁿC_k 可以排成帕斯卡三角的形式，如下图：

可知，前 8 行中，总共有 12 个不同的数字 1, 2, 3, 4, 5, 6, 7, 10, 15, 20, 21 和 35。对于正整数 n，如果没有任何一个质数的平方可以整除 n，则 n 被叫做 squarefree 。
那么。帕斯卡三角的前 8 行得到的那 12 个不同的数字中，除了 4 和 20，其余都是 squarefree 的。所有 squarefree 数之和是 105。

请给出帕斯卡三角前 51 行中不同的 squarefree 数之和。

jerryxjr1220

1. import time
   
2. st=time.time()
   
3. 
   
4. from math import factorial as f
   
5. def C(n,k):
   
6.         return f(n)//f(k)//f(n-k)
   
7. 
   
8. # from functools import lru_cache
   
9. # @lru_cache(maxsize=None)
   
10. # def C(n,k):
    
11. #         if k==0 or k==n: return 1
    
12. #         return C(n-1,k-1)+C(n-1,k)
    
13. 
    
14. def gen_prime_sqr(n):
    
15.         primes = [True]*n
    
16.         primes[:2] = [False]*2
    
17.         for i, prime in enumerate(primes):
    
18.                 if prime:
    
19.                         for j in range(i*i, n, i):
    
20.                                 primes[j] = False
    
21.         return [k*k for k, trueprime in enumerate(primes) if trueprime]
    
22. 
    
23. prime_sqr = gen_prime_sqr(int(C(50,25)**0.25)+1)
    
24. 
    
25. def is_sqrfree(n, prime_sqr=prime_sqr):
    
26.         for p in prime_sqr:
    
27.                 if n % p == 0:
    
28.                         return False
    
29.         return True
    
30. 
    
31. sqrfreelist = set()
    
32. for n in range(0, 51):
    
33.         for k in range(0, (n+1)//2+1):
    
34.                 if is_sqrfree(C(n,k)):
    
35.                         sqrfreelist.add(C(n,k))
    
36. print(sum(sqrfreelist))
    
37. print(time.time()-st)

复制代码

34029210557338
0.015608787536621094

guosl

1. /*
   
2. 答案：34029210557338
   
3. 耗时：0.000266566秒
   
4. */
   
5. #include <iostream>
   
6. #include <vector>
   
7. #include <unordered_set>
   
8. #include <algorithm>
   
9. #include <cmath>
   
10. #include <omp.h>
    
11. using namespace std;
    
12. 
    
13. const int N = 50;
    
14. 
    
15. struct np
    
16. {
    
17.   int v;
    
18.   int k;
    
19.   bool operator==(const np &n) const
    
20.   {
    
21.     return (v == n.v);
    
22.   }
    
23. };
    
24. 
    
25. int p[N + 1];//p[i]记录i的最小素因数
    
26. vector<np> vP[N + 1];//vP[k]记录k!的素数分解
    
27. unordered_set<long long> uS;//记录无素因数的组合数
    
28. 
    
29. void NDiv(vector<np> &p1, const vector<np> &p2)
    
30. {
    
31.   for (auto itr = p2.begin(); itr != p2.end(); ++itr)
    
32.   {
    
33.     np tmp = *itr;
    
34.     auto itr1 = find(p1.begin(), p1.end(), tmp);
    
35.     if (itr1->k == tmp.k)
    
36.       p1.erase(itr1);
    
37.     else
    
38.       itr1->k -= tmp.k;
    
39.   }
    
40. }
    
41. 
    
42. int main(void)
    
43. {
    
44.   double t = omp_get_wtime();
    
45.   //应用筛法求一个数的最小素因数
    
46.   for (int i = 1; i <= N; ++i)
    
47.     p[i] = i;
    
48.   for (int i = 2; i <= (int)sqrt((double)N); ++i)
    
49.   {
    
50.     if (p[i] < i)
    
51.       continue;
    
52.     for (int j = i * i; j <= N; j += i)
    
53.     {
    
54.       if (p[j] > i)
    
55.         p[j] = i;
    
56.     }
    
57.   }
    
58.   //求k!的素因数分解
    
59.   for (int i = 2; i <= N; ++i)
    
60.   {
    
61.     vP[i] = vP[i - 1];
    
62.     int x = i;
    
63.     while (x > 1)
    
64.     {
    
65.       int y = p[x];
    
66.       int z = 1;
    
67.       x /= y;
    
68.       while (x % y == 0)
    
69.       {
    
70.         ++z;
    
71.         x /= y;
    
72.       }
    
73.       np tmp;
    
74.       tmp.v = y;
    
75.       tmp.k = z;
    
76.       auto itr = find(vP[i].begin(), vP[i].end(), tmp);
    
77.       if (itr == vP[i].end())
    
78.         vP[i].push_back(tmp);
    
79.       else
    
80.         itr->k += tmp.k;
    
81.     }
    
82.   }
    
83.   uS.insert(1);
    
84.   for (int i = 2; i <= N; ++i)//对行进行循环
    
85.   {
    
86.     for (int j = 1; j <= (i + 1) / 2; ++j)//对列循环
    
87.     {
    
88.       //计算组合数
    
89.       vector<np> vB = vP[i];
    
90.       NDiv(vB, vP[j]);
    
91.       NDiv(vB, vP[i - j]);
    
92.       //分析是否有平方的素因数
    
93.       bool bFlag = true;
    
94.       long long a = 1;
    
95.       for (auto itr = vB.begin(); itr != vB.end(); ++itr)
    
96.       {
    
97.         if (itr->k >= 2)
    
98.         {
    
99.           bFlag = false;
    
100.           break;
     
101.         }
     
102.         else
     
103.           a *= itr->v;
     
104.       }
     
105.       if (bFlag)
     
106.         uS.insert(a);
     
107.     }
     
108.   }
     
109.   long long nS = 0;
     
110.   for (auto itr = uS.begin(); itr != uS.end(); ++itr)
     
111.     nS += *itr;
     
112.   t = omp_get_wtime() - t;
     
113.   cout << nS << endl << t << endl;
     
114.   return 0;
     
115. }
     

复制代码