标签：11 right return Network LibreOJ sum uint ull left

Description

给定 \(N,k\)，求：

\[\sum_{i=1}^N\sum_{j=1}^N (f(\gcd(i,j)))^k \]

其中 \(f(x)\) 表示 \(x\) 质因子分解结果中次大的质因子，重复的质因数计算多次。

Constraints

\(1\le N,k\le 2\cdot 10^9\)

Solution

记 \(f_k(x)=(f(x))^k\)。推式子：

\[\begin{aligned} \sum_{i}^N\sum_j^N f_k(\gcd(i,j)) &= \sum_g^N f_k(g)\sum_{i}^{N/g}\sum_{j}^{N/g} [\gcd(i,j)=1] \\ &= \sum_g^N f_k(g)\sum_d^{N/g} \mu(d)\left\lfloor\frac{N}{dg}\right\rfloor^2 \\ &= \sum_g^N f_k(g)G\left(\left\lfloor\frac{N}{g}\right\rfloor\right) \\ G(N) &= \sum_{d}^N \mu(d)\lfloor N/d\rfloor^2 \end{aligned} \]

随便什么筛出 \(\mu\)，min-25 用可以用 UOJ Sanrd 的套路筛出 \(f_k\)​。

整除分块套整除分块复杂度 \(O\left( \int_1^\sqrt{N}\sqrt i\right)+O\left( \int_1^\sqrt{N}\sqrt \frac{N}{i} \right) = O(N^{3/4})\)。题中 \(N\le 2\cdot 10^9\)​​，算比较小的，可以直接过。

#include <algorithm>
#include <cmath>
#include <cstdio>
using uint = unsigned int;
using ull = unsigned long long;

inline uint fastpow(uint a, uint b) {
  uint r = 1;
  for (; b; b >>= 1, a *= a)
    if (b & 1) r *= a;
  return r;
}

const int N = 1e5;
uint n, K;

namespace min_25 {
  uint n, m, tot, cnt, K;
  uint D[N], p[N], pK[N];
  bool vis[N];

  inline uint id(uint x) {
    return x <= m ? cnt - x + 1 : n / x;
  }

  void sieve(uint n) {
    vis[1] = 1;
    for (uint i = 1; i <= n; i++) {
      if (!vis[i]) p[++tot] = i;
      for (uint j = 1; j <= tot && i * p[j] <= n; j++) {
        vis[i * p[j]] = 1;
        if (i % p[j] == 0) break;
      }
    }
    for (uint i = 1; i <= tot; i++)
      pK[i] = fastpow(p[i], K);
  }

  uint s[N];
  void init_s() {
    for (uint i = 1; i <= cnt; i++)
      s[i] = D[i] - 1;
    for (uint j = 1; j <= tot; j++)
      for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
        s[i] -= s[id(D[i] / p[j])] - (j - 1);
  }

  namespace fk {
    uint g[N];
    void init_g() {
      for (uint i = 1; i <= cnt; i++)
        g[i] = s[i];
      for (uint j = tot; j; --j)
        for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
          for (uint pk = p[j]; (ull)pk * p[j] <= D[i]; pk *= p[j])
            g[i] += g[id(D[i] / pk)] - s[id(D[i] / pk)]
                  + pK[j] * (s[id(D[i] / pk)] - j + 1);
    }
    inline uint get(uint x) {
      if (!x) return 0u;
      return g[id(x)];
    } 
  }
  
  namespace mu {
    uint g[N];
    void init_g() {
      for (uint i = 1; i <= cnt; i++)
        g[i] = -s[i];
      for (uint j = tot; j; --j)
        for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
          g[i] -= g[id(D[i] / p[j])] + j;
    }
    inline uint get(uint x) {
      if (!x) return 0u;
      return g[id(x)] + 1;
    } 
  }

  void load(uint _n, uint _K) {
    n = _n, K = _K, m = sqrt(n);
    for (uint x = 1; x <= n; x = n / (n / x) + 1)
      D[++cnt] = n / x;
    sieve(m);
    
    init_s();
    mu::init_g();
    fk::init_g();
  }
}

signed main() {
  scanf("%u%u", &n, &K);
  min_25::load(n, K);

  auto calc_S = [&](uint n) {
    using min_25::mu::get;
    uint res = 0;
    for (uint l = 1, r; l <= n; l = r + 1) {
      r = n / (n / l);
      res += (get(r) - get(l - 1)) * (n / l) * (n / l);
    }
    return res;
  };

  auto calc_F = [&](uint n) {
    using min_25::fk::get;
    uint res = 0;
    for (uint l = 1, r; l <= n; l = r + 1) {
      r = n / (n / l);
      res += calc_S(n / l) * (get(r) - get(l - 1));
    }
    return res;
  };

  printf("%u\n", calc_F(n));
  return 0;
}

Solution Plus 1

做完发现有更优做法！考虑迪利克雷卷积：\(h=f_k*\mu\)。那么答案成了：

\[\sum_i^N h(i)\left\lfloor\frac{N}{i}\right\rfloor^2 \]

式子非常简洁，尝试直接求 \(h\) 的前缀和。由于不是积性函数，而且要比 \(f_k\) 的魔改方法复杂得多，min-25 筛直接求 \(h\) 有些困难。而 \(h\) 和迪利克雷卷积关系较大，由人类智慧可得：

\[h*I = f_k*\mu *I=f_k*\varepsilon=f_k \]

凑出了卷积形式，考虑杜教筛。记 \(F, H\) 分别为 \(f_k,h\) 的前缀和。那么：

\[F(N)=\sum_{i=1}^N H(\lfloor N/i\rfloor) \Leftrightarrow H(N)=F(N)-\sum_{i=2}^N H(\lfloor N/i\rfloor) \]

杜教筛复杂度 \(O(N^{2/3})\)，总复杂度 \(O\left(\frac{N^{3/4}}{\log N}\right)\)，瓶颈在 min-25 求 \(F\)。

Solution Plus 2

经过 Mr_Spade 指导，可以不要杜教筛！

考虑目标在于快速求 \(H\)​​​。这里有一个船新科技：设阈值 \(n^{2/3}\)​，并以此分治：

• 对于 \(\le n^{2/3}\)​​ 的点值，线性筛出 \(f_k,\mu\)​ 的点值，暴力迪利克雷卷积求出 \(h\)​ 在 \(n^{2/3}\)​​ 内所有的点值，最后前缀和求出 \(H\)。复杂度 \(O(n^{2/3}\log n)\)​。
• 对于 \(>n^{2/3}\)​​ 的点值​，使用之前整除分块的做法算出 \(H\)​ 在较大整除点的点值。复杂度 \(O\left(\int_1^\sqrt[3]N \sqrt\frac{N}{i}\right)=O(n^{2/3})\)。

最后一次整除分块就做完了。复杂度仍然是 \(O\left(\frac{N^{3/4}}{\log N}\right)\)。常数比杜教筛小很多！

#include <algorithm>
#include <cmath>
#include <cstdio>
using uint = unsigned int;
using ull = unsigned long long;

inline uint fastpow(uint a, uint b) {
  uint r = 1;
  for (; b; b >>= 1, a *= a)
    if (b & 1) r *= a;
  return r;
}

const int N = 1e5;
const int L = 1e7;
uint n, K, gap;
uint sfk[L], smu[L], sh[L];

namespace min_25 {
  uint n, m, tot, cnt, K;
  uint D[N], p[L / 10], pK[L / 10];
  bool vis[L];

  inline uint id(uint x) {
    return x <= m ? cnt - x + 1 : n / x;
  }

  void sieve(uint n) {
    vis[1] = 1;
    sfk[1] = 0, smu[1] = 1;
    for (uint i = 1; i <= n; i++) {
      if (!vis[i]) {
        p[++tot] = i;
        sfk[i] = 1, smu[i] = -1;
        pK[tot] = fastpow(i, K);
      }
      for (uint j = 1; j <= tot && i * p[j] <= n; j++) {
        vis[i * p[j]] = 1;
        if (vis[i]) sfk[i * p[j]] = sfk[i];
        else
          sfk[i * p[j]] = pK[j];
        if (i % p[j] == 0) break;
        smu[i * p[j]] = -smu[i];
      }
    }
    while (p[tot] > m) --tot;
  }

  uint s[N];
  void init_s() {
    for (uint i = 1; i <= cnt; i++)
      s[i] = D[i] - 1;
    for (uint j = 1; j <= tot; j++)
      for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
        s[i] -= s[id(D[i] / p[j])] - (j - 1);
  }

  namespace fk {
    uint g[N];
    void init_g() {
      for (uint i = 1; i <= cnt; i++)
        g[i] = s[i];
      for (uint j = tot; j; --j)
        for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
          for (uint pk = p[j]; (ull)pk * p[j] <= D[i]; pk *= p[j])
            g[i] += g[id(D[i] / pk)] - s[id(D[i] / pk)]
                  + pK[j] * (s[id(D[i] / pk)] - j + 1);
    }
    inline uint get(uint x) {
      if (!x) return 0u;
      return g[id(x)];
    } 
  }
  
  namespace mu {
    uint g[N];
    void init_g() {
      for (uint i = 1; i <= cnt; i++)
        g[i] = -s[i];
      for (uint j = tot; j; --j)
        for (uint i = 1; i <= cnt && D[i] >= (ull)p[j] * p[j]; i++)
          g[i] -= g[id(D[i] / p[j])] + j;
    }
    inline uint get(uint x) {
      if (!x) return 0u;
      return g[id(x)] + 1;
    } 
  }

  void load(uint _n, uint _K) {
    n = _n, K = _K, m = sqrt(n);
    for (uint x = 1; x <= n; x = n / (n / x) + 1)
      D[++cnt] = n / x;
    sieve(std::max(gap, m));
    
    init_s();
    mu::init_g();
    fk::init_g();
  }
}

signed main() {
  scanf("%u%u", &n, &K);
  gap = pow(n, 2.0 / 3);
  min_25::load(n, K);

  for (uint i = 1; i <= gap; i++)
    for (uint j = 1; i * j <= gap; j++)
      sh[i * j] += sfk[i] * smu[j];
  for (uint i = 2; i <= gap; i++)
    sh[i] += sh[i - 1];
  
  auto calc_H = [&](uint x) {
    auto mu = min_25::mu::get;
    auto fk = min_25::fk::get;
    if (x <= gap) return sh[x];
    uint ret = 0;
    for (uint l = 1, r; l <= x; l = r + 1) {
      r = x / (x / l);
      ret += (mu(r) - mu(l - 1)) * fk(x / l);
    }
    return ret;
  };

  uint ans = 0, lst = 0;
  for (uint l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    uint cur = calc_H(r);
    ans += (cur - lst) * (n / l) * (n / l);
    lst = cur;
  }

  printf("%u\n", ans);
  return 0;
}

标签：11,right,return,Network,LibreOJ,sum,uint,ull,left
来源： https://www.cnblogs.com/-Wallace-/p/loj572.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。