优美函数 题解

题意简述

给你函数 $f$：

$$f(x,y,u)=\{\begin{array}{rc}u-y,& x=1\\ u,& 1<x\le y,gcd(x,y)=1\\ -x\cdot y,& x\ne 1,gcd(x,y)=x\\ 0,& \text{otherwise}.\end{array}$$

对于一个长度为 $n$ 的序列 $a$，定义函数 $g$：

$$g(l,r,u)=\sum _{i=1}^{{10}^6}\sum _{j=l}^{r}f(i,a_j,u)$$

给你 $q$ 个询问，告诉你 $(l,u)$，求 $\underset{r=l}{\overset{n}{max}}g(l,r,u)$，以及对应的 $r$，多个最值输出 $r$ 较小的。

$n,q\le 5\times {10}^5$，$1\le a_i\le {10}^6$，$1\le u\le 1.8\times {10}^7$。

题目分析

显然要先观察 $f$，我们把 $\sum _{i=1}^{{10}^6}f(i,a_j,u)$ 的代码写出来：

lint res = 0;
for (int i = 1; i <= 1000000; ++i) {
    if (i == 1) res += u - a[j];
    else if (i <= a[j] && gcd(i, a[j]) == 1) res += u;
    else if (gcd(i, a[j]) == i) res += -i * a[j];
}

显然要把 $i=1$ 提出来，然后后面两个 if 是不交的，可以分别考虑。发现此时 $i$ 的上界可以下调至 $a_j$。$gcd(i,a_j)=i$ 等价于 $i\mid a_j$。再把不变的东西提出来，可以得到如下代码：

lint res = u - a[j];
lint cnt = 0, sum = 0;
for (int i = 2; i <= a[j]; ++i) {
    if (gcd(i, a[j]) == 1) ++cnt;
    if (a[j] % i == 0) sum += i;
}
res += cnt * u, res -= sum;

我们惊奇地发现，$u-a_j$ 满足我们 for 里面两个 if 的条件，可以放回去：

lint cnt = 0, sum = 0;
for (int i = 1; i <= a[j]; ++i) {
    if (gcd(i, a[j]) == 1) ++cnt;
    if (a[j] % i == 0) sum += i;
}
lint res = cnt * u - sum;

这时候，$cnt$ 和 $sum$ 就有意义了，分别表示：「$1\sim a_j$ 中和 $a_j$ 互质的数的个数」、「$a_j$ 的所有约数之和」，也就是 $\phi (a_j)$ 和 $\sigma (a_j)$，这两个都可以线性筛地搞。

可是我不会

考场上忘了，只好现推一遍。其实以前根本不知道啥是 $\sigma$，但是实力地推出式子来了！我们只用考虑 $imod{p}_j\ne 0$ 和 $imod{p}_j=0$ 两种情况。

1. $\phi (x)$。

   1. $imod{p}_j\ne 0$。
      经典积性函数曾说：如果你给我 $a\perp b$，我能告诉你，$\phi (ab)=\phi (a)\phi (b)$。所以 $\phi (i\cdot p_j)=\phi (i)\phi (p_j)$。
   2. $imod{p}_j=0$。
      我们记 $i=p_j^t\prod _{p_k\mid i}{p}_k^{t_k}$。考虑欧拉函数计算式 $\phi (x)=x\prod _{p_j\mid x}{\displaystyle \frac{p_j-1}{p_j}}$，就有 $\phi (i)=\phi \left({\displaystyle \frac{i}{p_j^t}}\right)p_j^t{\displaystyle \frac{p_j-1}{p_j}}$。而 ${\displaystyle \frac{i}{p_j^t}}\perp p_j^{t+1}$，所以：
   $$\begin{array}{rl}& \phi (i\cdot p_j)\\ & =\phi \left({\displaystyle \frac{i}{p_j^t}}{p}_j^{t+1}\right)\\ & =\phi \left({\displaystyle \frac{i}{p_j^t}}\right)\phi (p_j^{t+1})\\ & ={\displaystyle \frac{\phi (i)}{p_j^t}}{\displaystyle \frac{p_j}{p_j-1}}(p_j^{t+1}-p_t^t)\\ & ={\displaystyle \frac{\phi (i)}{\cancel{p_j^t}}}{\displaystyle \frac{p_j}{\cancel{p_j-1}}}\cancel{(p_j-1)}\cancel{p_j^t}\\ & =\phi (i)p_j\end{array}$$

2. $\sigma (x)$。
   考虑 $x=\prod _{p_k\mid x}{p}_k^{t_k}$，那么 $\sigma (x)=\mathrm{sum}\{\prod p_k^{t_k^{\prime }}\mid t_k^{\prime }\le t_k\}$，由于 $t_k^{\prime }$ 之间互不影响，考虑乘法原理，$\sigma (x)=\prod (1+p_k+\dots +p_k^{t_k})$，这样一定能凑出所有 $x$ 的因数。

   1. $imod{p}_j\ne 0$。
      $\sigma (i\cdot p_j)=\sigma (i)(1+p_j)=\sigma (i)\sigma (p_j)$。这里没有一步到位的原因是，我之前不知道 $\sigma$，也就没想它是不是积性函数了。
   2. $imod{p}_j=0$。
      同样设 $i$ 中有 $p_j^t$，那么 $\sigma (i\cdot p_j)=\sigma \left({\displaystyle \frac{i}{p_j^t}}\right)(1+p_j+\cdots +p_j^{t+1})$，这很难受，考虑变换一下 $\sigma \left({\displaystyle \frac{i}{p_j^t}}\right)(1+p_j+\cdots +p_j^{t+1})=\sigma (i)+\sigma \left({\displaystyle \frac{i}{p_j^t}}\right)p_j^{t+1}$。接下来，我们只需要对一个 $x$ 维护 $c(x)$ 表示 $p^t$，其中 $p$ 是 $x$ 的最小质因数，$t$ 是幂次。这个很好维护，当 $imod{p}_j\ne 0$，$c(i\cdot p_j)=p_j$；否则，$c(i\cdot p_j)=c(i)p_j$。如此，$\sigma (i\cdot p_j)=\sigma (i)+\sigma \left({\displaystyle \frac{i}{c(i)}}\right)c(i)p_j$

对于一次询问，求的式子可以简化成：

$$\underset{r=l}{\overset{n}{max}}\{\sum _{i=l}^r\phi (a_i)\cdot u+\sigma (a_i)\}$$

对于一个 $i$ 来说，$\phi (a_i)$ 和 $\sigma (a_i)$ 都是定值，我们可以把它们分别做个前缀和。由于很像一次函数，我们将这两个前缀和记作 $K,B$。

$$\begin{array}{rl}& \underset{r=l}{\overset{n}{max}}\{K_r-K_{l-1}\cdot u+B_r-B_{l-1}\}\\ & =\underset{r=l}{\overset{n}{max}}\{K_r\cdot u+B_r\}-K_{l-1}\cdot u-B_{l-1}\end{array}$$

一次函数最值，且 $K$ 有单调性，我们考虑凸包上二分。我们考虑把查询离线到 $l$ 上，从右往左扫，维护一个下凸包（一次函数的最值），查询最值就去二分出 $x$ 在凸包哪条直线上。

时间复杂度：$\mathrm{\Theta }(V+n+q\log n)$。

代码

#include <cstdio>
#include <iostream>
#include <vector>
using namespace std;
 
const int N = 500010, V = 1000010;
 
using lint = long long;
 
int n, q, val[N];
int phi[V], prime[V], pcnt, cj[V];
bool nprime[V];
lint subsum[V], K[N], B[N];
vector<pair<int, int>> qry[N];
lint mx[N];
int mxR[N];
 
int stack[N], top;
 
inline bool cmptail(int a, int b, int c) {
    return (__int128_t)(B[c] - B[a]) * (K[a] - K[b]) > (__int128_t)(B[b] - B[a]) * (K[a] - K[c]);
}
 
inline void insert(int idx) {
    while (top - 1 >= 1 && cmptail(stack[top - 1], stack[top], idx)) --top;
    stack[++top] = idx;
}
 
inline int query(int x) {
    int l = 1, r = top - 1, res = stack[top];
    while (l <= r) {
        int mid = (l + r) >> 1;
        if (1ll * x * (K[stack[mid]] - K[stack[mid + 1]]) >= B[stack[mid + 1]] - B[stack[mid]])
            res = stack[mid], r = mid - 1;
        else
            l = mid + 1;
    }
    return res;
}
 
signed main() {
    #ifndef XuYueming
    freopen("function.in", "r", stdin);
    freopen("function.out", "w", stdout);
    #endif
    scanf("%d%d", &n, &q);
    for (int i = 1; i <= n; ++i)
        scanf("%d", &val[i]);
    phi[1] = 1, subsum[1] = 1;
    for (int i = 2; i <= 1000000; ++i) {
        if (!nprime[i]) prime[++pcnt] = i, phi[i] = i - 1, subsum[i] = 1 + i, cj[i] = i;
        for (int j = 1; j <= pcnt && i * prime[j] <= 1000000; ++j) {
            nprime[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                subsum[i * prime[j]] = subsum[i] + subsum[i / cj[i]] * cj[i] * prime[j];
                cj[i * prime[j]] = cj[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * phi[prime[j]];
            subsum[i * prime[j]] = subsum[i] * subsum[prime[j]];
            cj[i * prime[j]] = prime[j];
        }
    }
    for (int i = 1; i <= n; ++i) {
        K[i] = phi[val[i]];
        B[i] = -1ll * subsum[val[i]] * val[i];
        K[i] += K[i - 1];
        B[i] += B[i - 1];
    }
    for (int i = 1, u, l; i <= q; ++i) {
        scanf("%d%d", &u, &l);
        qry[l].emplace_back(u, i);
    }
    for (int i = n; i >= 1; --i) {
        insert(i);
        for (auto [x, idx]: qry[i]) {
            int ps = query(x);
            mx[idx] = (1ll * x * K[ps] + B[ps]) - (1ll * x * K[i - 1] + B[i - 1]);
            mxR[idx] = ps;
        }
    }
    for (int i = 1; i <= q; ++i)
        printf("%lld %d\n", mx[i], mxR[i]);
    return 0;
}