UOJ62【UR#5】怎样跑得更快 <莫比乌斯反演>

Problem

【UR #5】怎样跑得更快

时间限制： $\mathrm{1\;Sec}$
空间限制： $\mathrm{256\;MB}$
大力水手问禅师：“大师，我觉得我光有力气是不够的。比如我吃菠菜可以让力气更大，但是却没有提升跑步的速度。请问怎样才能跑得更快？我试过吃白菜，没有效果。”
禅师浅笑，答：“方法很简单，不过若想我教你，你先看看这道$\mathrm{UOJ\;Round}$的$\mathrm{C}$题。”
令 $p=998244353$（$7\times 17\times 223+1$，一个质数）。
给你整数 $n,c,d$。现在有整数 $x_1,\cdots,x_n$ 和 $b_1,\cdots,b_n$ 满足 $0\le x_1,\cdots,x_n,b_1,\cdots,b_n<p$，且对于 $1\le i\le n$ 满足：

$$\sum_{j=1}^{n}\gcd(i,j)^c\cdot\mathrm{lcm}(i,j)^d\cdot x_j\equiv b_i\mod{p}$$

其中 $v\equiv u\mod{p}$ 表示 $v$ 和 $u$ 除以 $p$ 的余数相等。$\gcd(i,j)$ 表示 $i$ 和 $j$ 的最大公约数，$\mathrm{lcm}(i,j)$ 表示 $i$ 和 $j$ 的最小公倍数。
有 $q$ 个询问，每次给出 $b_1,\cdots,b_n$，请你解出 $x_1,\cdots,x_n$ 的值。

输入格式

第一行四个整数 $n,c,d,q$。保证 $n,q\ge1$。
接下来 $q$ 行，每行 $n$ 个整数依次表示 $b_1,\cdots,b_n$。保证 $0\le b_1,\cdots,b_n<p$。

输出格式

共 $q$ 行，每行对给出的 $b_1,\cdots,b_n$，输出对应的 $x_1,\cdots,x_n$。
如果有多组解输出任意一组即可。
如果无解那么这一行只用输出一个整数 $-1$。

样例输入输出

样例一

Input

3 1 0 2
1 0 0
1 2 3

Output

499122179 998244352 499122176
998244352 1 1

Explanation
对于第一个询问，要满足的等式为：

$$\begin{cases} x_1+x_2+x_3\equiv1\mod{p}\\ x_1+2x_2+x_3\equiv0\mod{p}\\ x_1+x_2+3x_3\equiv0\mod{p}\\ \end{cases}$$

样例二

见样例数据下载。

限制与约定

对于所有数据，$n\cdot q\le3\times10^5$，$0\le c,d\le10^9$。

测试点编号	$n$	$q$	其他
$1$	$\le3$	$=1$	$c,d\le1000$
$2$	$\le100$	$=1$	$c=d$
$3\sim4$	$\le100$	$\le10$	保证有唯一解
$5$	$\le100$	$\le1000$	$c=1,d=0$
$6\sim8$	$\le100$	$\le1000$	保证有唯一解
$9\sim11$	$\le1000$	$=10$	保证有唯一解
$12$	$\le1000$	$\le10$	$c=d$
$13\sim14$	$\le1000$	$\le10$	保证有唯一解
$15\sim16$	$\le10^5$	$\le3$	$c=1,d=0$
$17$	$\le10^5$	$\le3$	$c=d$
$18\sim20$	$\le10^5$	$\le3$	无

后记

还没听完题，大力水手就嘶吼着：“太难了我不会我不会！”，飞快地跑掉了。
禅师看着大力水手消失的背影，叹了口气说：“你们这些人啊，每天就想做些大水题，一碰到难题，跑得不知道比谁都快。”
后来大力水手把$\mathrm{UOJ\;Round}$的$\mathrm{C}$题题面贴在了汽车的后挡风玻璃上，人类从此掌握了光速旅行的正确方式。

下载

样例数据下载

标签：莫比乌斯反演

Solution

先膜一发vfk的题解。

首先将题目中的式子转换一下：

$$\sum_{j=1}^{n}\gcd(i,j)^c\cdot\mathrm{lcm}(i,j)^d\cdot x_j\equiv b_i\mod{p}\\ \sum_{j=1}^{n}\gcd(i,j)^{c-d}\cdot i^d\cdot j^d\cdot x_j\equiv b_i\mod{p}\\$$

设$f(x)=x^{c-d},\;g(x)=x^d$，那么原式化为

$$\sum_{j=1}^{n}f(\gcd(i,j))\cdot g(i)\cdot g(j)\cdot x_j\equiv b_i\mod{p}$$

此时可以强行”说一句废话“来加入反演。存在数论函数$f_r(n)$使得$f(n)=\sum_{d|n}f_r(d)$，那么$f_r(n)=f(n)-\sum_{d|n且d\ne n}f_r(d)$，这样若知道$f$，即可枚举约数求出$f_r$，即

for (int i = 1; i <= n; i++) fr[i] = f[i];
for (int i = 1; i <= n; i++)
    for (int j = i+i; j <= n; j += i)
        fr[j] -= fr[i];

将$f(n)=\sum_{d|n}f_r(d)$带入原式即可得到

$$\sum_{j=1}^{n}\sum_{d|i,j}f_r(d)\cdot g(i)\cdot g(j)\cdot x_j\equiv b_i\mod{p}\\ \sum_{d|i}f_r(d)\sum_{j为d的倍数且j\le n}g(j)\cdot x_j\equiv b_i\cdot g(i)^{-1}\mod{p}\\$$

令$z(d)=\sum_{j为d的倍数且j\le n}g(j)\cdot x_j$，$f_z(d)=f_r(d)\cdot z(d)$，那么

$$\sum_{d|i}f_z(d)\equiv b_i\cdot g(i)^{-1}\mod{p}\\$$

我们知道右边，要求左边，于是再次做反演，即$f_z(n)=b_n\cdot g(n)^{-1}-\sum_{d|n且d\ne n}f_z(d)\mod{p}$，即

for (int i = 1; i <= n; i++) fz[i] = b[i]*Pow(g(i), P-2);
for (int i = 1; i <= n; i++)
    for (int j = i+i; j <= n; j += i)
        fz[j] -= fz[i];

这样就得到了$f_z$的值，即得到$f_r\cdot z$的值，于是$z(n)\equiv f_z(n)\cdot f_r(n)^{-1}$即可算出$z$的值。
然后再尝试用$z$推出$x$。由于$z(d)=\sum_{j为d的倍数且j\le n}g(j)\cdot x_j$，设$hx(n)=g(n)\cdot x_n$，那么$\sum_{j为d的倍数且j\le n}hx(j)=z(d)$。
仍然是反演的形式，再次反演得到$hx(d)=z(d)-\sum_{j=1}^{n}[j>d][d|j]\cdot hx(j)$，即

for (int i = 1; i <= n; i++) hx[i] = z[d];
for (int i = n; i >= 1; i--)
    for (int j = i+i; j <= n; j += i)
        hx[i] -= hx[j];

知道了$hx$，即可由$x_i=hx(i)\cdot g(i)^{-1}$求出$x$的值。

总结一下，分为三次反演：

• 预处理$f,g$的值（线性筛），反演求$f_r$
• 读入$b$，反演求$f_z$，乘$f_r$的逆元即可得到$z$
• 通过$z$反演求$hx$，乘$g$的逆元即可得到$x$

Code

#include <bits/stdc++.h>
#define MAX_N 100000
#define P 998244353
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
int n, c, d, q, b[MAX_N+5], x[MAX_N+5];
lnt g[MAX_N+5], fr[MAX_N+5], fz[MAX_N+5], hx[MAX_N+5];
int cnt, pri[MAX_N+5]; bool NotPri[MAX_N+5];
inline lnt Pow(lnt x, int k) {
    lnt ret = 1; bool flag = (k < 0);
    for (k = abs(k); k; k >>= 1, (x *= x) %= P)
        if (k&1) (ret *= x) %= P;
    return flag ? Pow(ret, P-2) : ret;
}
void init() {
    fr[1] = g[1] = NotPri[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!NotPri[i])
            pri[cnt++] = i, fr[i] = Pow(i, c-d), g[i] = Pow(i, d);
        for (int j = 0; j < cnt; j++) {
            if (i*pri[j] > n) break;
            NotPri[i*pri[j]] = true;
            fr[i*pri[j]] = fr[i]*fr[pri[j]]%P;
            g[i*pri[j]] = g[i]*g[pri[j]]%P;
            if (i%pri[j] == 0) break;
        }
    }
    for (int i = 1; i <= n; i++)
        for (int j = i+i; j <= n; j += i)
            (fr[j] += (P-fr[i])) %= P;
}
int main() {
    read(n), read(c), read(d), read(q), init();
    for (bool flag = false; q--; flag = false) {
        for (int i = 1; i <= n; i++) read(b[i]);
        for (int i = 1; i <= n; i++)
            fz[i] = 1LL*b[i]*Pow(g[i], P-2)%P;
        for (int i = 1; i <= n; i++)
            for (int j = i+i; j <= n; j += i)
                (fz[j] += (P-fz[i])) %= P;
        for (int i = 1; i <= n; i++)
            if (!fr[i] && fz[i]) flag = true;
        if (flag) {puts("-1"); continue;}
        for (int i = 1; i <= n; i++)
            hx[i] = fz[i]*Pow(fr[i], P-2)%P;
        for (int i = n; i >= 1; i--)
            for (int j = i+i; j <= n; j += i)
                (hx[i] += (P-hx[j])) %= P;
        for (int i = 1; i <= n; i++)
            x[i] = (int)(hx[i]*Pow(g[i], P-2)%P);
        for (int i = 1; i <= n; i++) printf("%d ", x[i]); puts("");
    }
    return 0;
}