标签归档:筛法

[2018 Multi-University Training Contest 10] C. Calculate

Link: http://acm.hdu.edu.cn/showproblem.php?pid=6428

Statement:

Given \(1 \leq A, B, C \leq 10^7\), compute

\[ \sum_{i=1}^A \sum_{j=1}^B \sum_{k=1}^C \phi(\gcd(i, j^2, k^3)) \]

no more than 10 test cases.

Solution:

Consider the number of triples (i, j, k) such that \( \gcd(i, j^2, k^3) = d\), denoted as \(f(d)\). By Mobius inversion,

\[ f(d) = \sum_{i} \mu(i) g(di) \]

where \(g(d)\) is the number of triples satisfying \( d | \gcd(i, j^2, k^3) \), and the answer should be

\[ \sum_{d=1}^A f(d) \phi(d) \]

Now we try to compute \(g(d)\):

\[ g(d) = \sum_{i = 1}^A [d | i] \sum_{j = 1}^B [d | j^2] \sum_{k = 1}^C [d | k^3] \]

Let \(d = \prod {p_i}^{a_i} \) denote the unique factorization of d. Note that \(d | k^n \) whenever \( \prod {p_i}^{\lceil a_i / n \rceil} | k \).

Let \( h_n(d) = \prod {p_i}^{\lceil a_i / n \rceil} \), which denotes the minimum positive integer whose nth power is divisible by d. Now \(g(d)\) could be rewritten as

\[ g(d) = \lfloor \frac{A}{d} \rfloor \lfloor \frac{B}{h_2(d)} \rfloor \lfloor \frac{C}{h_3(d)} \rfloor \]

Plug in answer:

\[ \sum_{d=1}^A f(d) \phi(d) = \sum_{d=1}^A \phi(d) \sum_{i} \mu(i) g(di) = \sum_{q = 1}^A g(q) \sum_{ij = q} \mu(i) \phi(j) \]

where the Dirichlet convolution \( (\mu * \phi)(q) = \sum_{ij = q} \mu(i) \phi(j) \) is multiplicative and thus could be effectively precomputed (along with \(h_2\), \(h_3\)) by the sieve of Euler in linear time. Hence each test case could be done in \(O(A)\) time.

Source Code:

#include <bits/stdc++.h>
using namespace std;

#ifdef __LOCAL_DEBUG__
# define _debug(fmt, ...) fprintf(stderr, "\033[94m%s: " fmt "\n\033[0m", \
    __func__, ##__VA_ARGS__)
#else 
# define _debug(...) ((void) 0)
#endif
#define rep(i, n) for (int i=0; i<(n); i++)
#define Rep(i, n) for (int i=1; i<=(n); i++)
#define range(x) (x).begin(), (x).end()
typedef long long LL;
typedef unsigned long long ULL;

namespace sieve {
  constexpr int MAXN = 10000007;
  bool p[MAXN];
  int prime[MAXN], sz;
  int pval[MAXN], pcnt[MAXN];
  unsigned h2[MAXN], h3[MAXN], phi[MAXN], muphi[MAXN];

  void exec(int N = 10000007) {
    p[0] = p[1] = 1;
    
    pval[1] = 1;
    pcnt[1] = 0;
    h2[1] = h3[1] = phi[1] = muphi[1] = 1;

    for (int i = 2; i < N; i++) {
      if (!p[i]) {
        prime[sz++] = i;
        for (LL j = i; j < N; j *= i) {
          int b = j / i;
          pval[j] = i * pval[b];
          pcnt[j] = pcnt[b] + 1;
          h2[j] = h2[b] * (pcnt[j] % 2 == 1 ? i : 1);
          h3[j] = h3[b] * (pcnt[j] % 3 == 1 ? i : 1);
          phi[j] = (i == j ? i-1 : phi[b] * i);
          muphi[j] = phi[j] - phi[b];
        }
      }
      for (int j = 0; i * prime[j] < N; j++) {
        int x = i * prime[j]; p[x] = 1;
        if (i % prime[j] == 0) {
          pval[x] = pval[i] * prime[j];
          pcnt[x] = pcnt[i] + 1;
        } else {
          pval[x] = prime[j];
          pcnt[x] = 1;
        }
        if (x != pval[x]) {
          int b = x / pval[x];
          h2[x] = h2[b] * h2[pval[x]];
          h3[x] = h3[b] * h3[pval[x]];
          muphi[x] = muphi[b] * muphi[pval[x]];
        }
        if (i % prime[j] == 0) break;
      }
    }
  }
}

int main() {
  using namespace sieve;
  exec();
  int T; scanf("%d", &T);
  while (T--) {
    unsigned a, b, c; scanf("%d%d%d", &a, &b, &c);
    unsigned res = 0;
    for (unsigned i = 1; i <= a; i++) 
      res += muphi[i] * (a / i) * (b / h2[i]) * (c / h3[i]); 
    printf("%d\n", res & 0x3fffffff);
  }
  return 0;
}

 

k次幂前缀和的线性做法

给定正整数 \(x, k\) 和质数 \(p\),求
\[ \sum_{i=1}^{x} i^k \bmod p \]

这是一道很经典的题目,在杜瑜皓的《多项式及求和》中有详细讨论。比较暴力的做法是对于每个 \(i^k\) 都用快速幂算出其值,然后相加,时间复杂度为 \(O(x \log k) \),当 \(x\) 很大时,就超时了。然而利用拉格朗日插值法,这题最快可以做到 \(O(k)\) 的复杂度,其做法主要由以下几部分构成。

线性预处理逆元

求 1..n 之间每个数的逆元,如果都用费马小定理或者扩展欧几里得算,那么复杂度将会达到 \(O(n \log p)\) 或 \(O(n \log n)\)。利用一些递推式,可以线性地求出 1..n 中每个数的逆元,从而复杂度可以减少一个 log。常用的一个递推公式是

\[ i^{-1} = -\left\lfloor \frac{p}{i} \right\rfloor \cdot (p \bmod i)^{-1} \bmod p \]

这一公式的正确性也很好证明,等式右端可以写成

\[ -\frac{\left\lfloor \frac{p}{i} \right\rfloor}{p-\left\lfloor \frac{p}{i} \right\rfloor i} = \frac{1}{i} \]

利用此公式,可以在 \(O(n)\) 时间内求出 1..n 中每一个数的逆元。

线性筛求积性函数的值

注意到 \(i^k\) 是一个积性函数,这样,利用线性筛,我们只需要用快速幂计算素数 \(i\) 的值,其他点的函数值就可以直接得到。根据素数定理,前 \(n\) 个数中素数的个数约为 \(O(\frac{n}{\log n})\),从而利用快速幂计算出素数处的值,可以做到 \(O(\frac{n}{\log n}\log{k}\) ,再加上筛法本身 \(O(n)\) 的复杂度,我们可以在 \(O(n + \frac{n}{\log n}\log{k})\) 的时间复杂度内计算出 \(i^k\) 的前 \(n\) 项值。

拉格朗日插值法

注意到答案一定是关于 \(x\) 的 \(k+1\) 次多项式。而根据线性代数的知识,我们只要知道这个多项式的任意 \(k+2\) 个点处的值,我们就能确定这个多项式。拉格朗日插值公式就给出了这个多项式

\[ f(x) = \sum_{i=0}^{k+1} y_i l_i(x) \]

其中

\[ l_i(x) = \prod_{0 \leq m \leq k+1, m \neq i} \frac{x – x_m}{x_j – x_m} \]

具体到这道题上,插值公式就是

\[ f(n) = \sum_{i = 0}^{k+1} (-1)^{k+i+1} f(i) \frac{n(n-1) \cdots (n-i+1)(n-i-1) \cdots (n-d)}{(d-i)!i!} \]

其中 \(f(i)\) 就是 \(i^k\) 的前缀和,上式右端求和的每一项,分子都是 n 到 n-d 的前缀积乘后缀积,分母可以用递推线性求出,这样就可以在 \(O(k)\) 时间内解决整个问题。

#include <bits/stdc++.h>
using namespace std;

#define rep(i, n) for (int i=0; i<(n); i++)
#define Rep(i, n) for (int i=1; i<=(n); i++)
#define range(x) (x).begin(), (x).end()
typedef long long LL;
typedef unsigned long long ULL;

#define pow owahgrauhgso

ULL n, k, m;

ULL powmod(ULL b, ULL e) {
  ULL r = 1;
  while (e) {
    if (e & 1) r = r * b % m;
    b = b * b % m;
    e >>= 1;
  }
  return r;
}

const int MAXN = 1000006;
ULL inv[MAXN];

void init_inv() {
  inv[1] = 1;
  for (int i = 2; i <= k+1; i++) { 
    inv[i] = (m - m / i) * inv[m % i] % m; 
    assert(inv[i] * i % m == 1);
  }
}

ULL pow[MAXN];
ULL prime[MAXN], cnt;

void sieve() {
  pow[1] = 1;
  for (int i = 2; i <= k+1; i++) {
    if (!pow[i]) {
      pow[i] = powmod(i, k);
      prime[cnt++] = i;
      for (int j = 0; j < cnt && i*prime[j] <= k+1; j++) {
        pow[i * prime[j]] = pow[i] * pow[prime[j]] % m;
        if (i % prime[j] == 0) break;
      }
    }
  }
}

ULL sum[MAXN];
ULL ans[MAXN];

auto addmod = [](ULL a, ULL b) -> ULL {return (a+b)%m;};

ULL lagrange() {
  ULL p;
  p = 1;
  for (int i=0; i<=k+1; i++) {
    if (i) p = p * inv[i] % m;
    ans[i] = (k+1-i)&1 ? m-sum[i] : sum[i];
    ans[i] = ans[i] * p % m;
    p = p * (m + n - i) % m;
  }
  p = 1;
  for (int i=k+1; i>=0; i--) {
    if (k+1-i) p = p * inv[k+1-i] % m;
    ans[i] = ans[i] * p % m;
    p = p * (m + n - i) % m;
  }
  return accumulate(ans, ans+k+2, 0, addmod); 
}

 

int main() { cin >> n >> k >> m; init_inv(); sieve(); partial_sum(pow, pow+k+2, sum, addmod); cout << lagrange() << endl; return 0; }

2017-2018年问题求解(三)期末上机考试 Day1 题解

A. 智能机器

题意:求 \(\sum_{i=l}^r d(i^k)\) 模 16290047,其中 \(d(n)\) 表示 \(n\) 的约数个数。

数据范围:\(r – l \leq 10^6, l \leq r \leq 10^{12}, k \leq 10^6\)。

题解:首先,对于任意正整数 \(n\),如果它的唯一分解为

$$ n = q_1^{a_1} q_2^{a_2} \cdots q_m^{a_m} $$

那么它的约数个数就是

$$ d(n) = (1+a_1)(1+a_2)\cdots(1+a_m) $$

从而

$$ d(n^k) = (1+ka_1)(1+ka_2)\cdots(1+ka_m) $$

这样,我们只要对每个要求的 \(i\),找出所有形如 \(p^x | i\) 中 \(x\) 最大的因子。注意到这样的因子不可能超过 \(\sqrt{i}\),除非 \(x = 1\)。从而我们对素数 \(p\),可以枚举 l~r 区间范围内的 \(rp^x\) 形式的数(其中 \(r\) 不是 \(p\) 的倍数),并乘上 \(kx+1\) 的贡献。如果某个数除以所有枚举过的 \(p^x\) 形式的因子仍然不为 1,说明剩下的数一定是一个素数,那么还要乘上 \(k+1\) 的贡献。

经过仔细的分析,上述三重循环枚举的时间复杂度实际上不会超过 \(O(\frac{\sqrt{r}}{\log r} \log (r – l))\),完全可以接受。

B.  Lunatic

题意:一共有 \(n\) 个动作,每个动作可以在 \([l_i, r_i]\) 之间的某个时刻完成,完成后可得到 \(v_i\) 分数。每一时刻只能做一样动作,计算最多能拿多少分。

数据范围:\(1 \leq n \leq 1000, 1 \leq s_i \leq t_i \leq 30000, 1 \leq v_i \leq 1000\)

题解:显然是带权二分图的最大权匹配。考场上没带 Kuhn-Munkres 的板子,只好敲费用流,最后 2 分钟卡过去的…(不过这题用裸的 KM 算法应该也过不了 233)其实这题标答是贪心匹配,将边按收益排序,优先匹配收益高的就行了…

说一下费用流的做法吧。每个动作连一条容量为 1,费用为 \(-v_i\) 的边到 \([l_i, r_i]\) 中的每个时间点,源点到每个动作连费用为 0 容量为 1 的边,每个时间点连费用为 0 容量为 1 的边到汇点,然后跑最小费用流取负就好了。不过这样是过不了全部数据的,要过全部数据,还要加两个优化:一是把时间区间离散化;二是如果某个动作有 \(r_i – l_i + 1 \geq n\),那么这个动作一定可以在某个时刻完成,就不需要加到图里面了。