题目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

import time
st=time.time()

from math import factorial as f
def C(n,k):
        return f(n)//f(k)//f(n-k)

# from functools import lru_cache
# @lru_cache(maxsize=None)
# def C(n,k):
#         if k==0 or k==n: return 1
#         return C(n-1,k-1)+C(n-1,k)

def gen_prime_sqr(n):
        primes = [True]*n
        primes[:2] = [False]*2
        for i, prime in enumerate(primes):
                if prime:
                        for j in range(i*i, n, i):
                                primes[j] = False
        return [k*k for k, trueprime in enumerate(primes) if trueprime]

prime_sqr = gen_prime_sqr(int(C(50,25)**0.25)+1)

def is_sqrfree(n, prime_sqr=prime_sqr):
        for p in prime_sqr:
                if n % p == 0:
                        return False
        return True

sqrfreelist = set()
for n in range(0, 51):
        for k in range(0, (n+1)//2+1):
                if is_sqrfree(C(n,k)):
                        sqrfreelist.add(C(n,k))
print(sum(sqrfreelist))
print(time.time()-st)

34029210557338
0.015608787536621094

guosl

/*
答案：34029210557338
耗时：0.000266566秒
*/
#include <iostream>
#include <vector>
#include <unordered_set>
#include <algorithm>
#include <cmath>
#include <omp.h>
using namespace std;

const int N = 50;

struct np
{
  int v;
  int k;
  bool operator==(const np &n) const
  {
    return (v == n.v);
  }
};

int p[N + 1];//p[i]记录i的最小素因数
vector<np> vP[N + 1];//vP[k]记录k!的素数分解
unordered_set<long long> uS;//记录无素因数的组合数

void NDiv(vector<np> &p1, const vector<np> &p2)
{
  for (auto itr = p2.begin(); itr != p2.end(); ++itr)
  {
    np tmp = *itr;
    auto itr1 = find(p1.begin(), p1.end(), tmp);
    if (itr1->k == tmp.k)
      p1.erase(itr1);
    else
      itr1->k -= tmp.k;
  }
}

int main(void)
{
  double t = omp_get_wtime();
  //应用筛法求一个数的最小素因数
  for (int i = 1; i <= N; ++i)
    p[i] = i;
  for (int i = 2; i <= (int)sqrt((double)N); ++i)
  {
    if (p[i] < i)
      continue;
    for (int j = i * i; j <= N; j += i)
    {
      if (p[j] > i)
        p[j] = i;
    }
  }
  //求k!的素因数分解
  for (int i = 2; i <= N; ++i)
  {
    vP[i] = vP[i - 1];
    int x = i;
    while (x > 1)
    {
      int y = p[x];
      int z = 1;
      x /= y;
      while (x % y == 0)
      {
        ++z;
        x /= y;
      }
      np tmp;
      tmp.v = y;
      tmp.k = z;
      auto itr = find(vP[i].begin(), vP[i].end(), tmp);
      if (itr == vP[i].end())
        vP[i].push_back(tmp);
      else
        itr->k += tmp.k;
    }
  }
  uS.insert(1);
  for (int i = 2; i <= N; ++i)//对行进行循环
  {
    for (int j = 1; j <= (i + 1) / 2; ++j)//对列循环
    {
      //计算组合数
      vector<np> vB = vP[i];
      NDiv(vB, vP[j]);
      NDiv(vB, vP[i - j]);
      //分析是否有平方的素因数
      bool bFlag = true;
      long long a = 1;
      for (auto itr = vB.begin(); itr != vB.end(); ++itr)
      {
        if (itr->k >= 2)
        {
          bFlag = false;
          break;
        }
        else
          a *= itr->v;
      }
      if (bFlag)
        uS.insert(a);
    }
  }
  long long nS = 0;
  for (auto itr = uS.begin(); itr != uS.end(); ++itr)
    nS += *itr;
  t = omp_get_wtime() - t;
  cout << nS << endl << t << endl;
  return 0;
}