[Nowcoder 2018ACM多校第十场E] Rikka with Equation

题目大意：
给出一组系数A， 和模数m
定义$f(A, m)$为模线性方程$\sum a_i * x_i = 0 \space (mod \space m)$ 的解数
给出一个长度为n的数组B, 对于B的所有非空子集A， 和一个属于[1,M]的m$w_m = \sum f(A, m) \space mod \space 998244353$， 求所有$w_m$的异或和。$(n, m, B_i \leq 10^5)$

题目思路：
先证明$f(A,m) = gcd(A, m) * m^{|A| - 1}$
考虑中国剩余定理， 先将m拆成质数幂的形式， 考虑对于某个单独质数幂的答案， 则原方程的解数就是各维度下解数的乘积。
考虑$p^t$， p是质数， $\sum a_i * x_i = 0 \space (mod \space p^t)$
由于p是质数， 对于$a_i$种与p互质的部分可以直接用逆元消去， 考虑$a_i$剩下的部分是$p^{q_i}$
原方程变为$\sum p^{q_i} * x_i = 0 \space (mod \space p^t)$
假设$p^{q_i}$种最小的是$p^{q_1}$， 将其余的部分看成一个整体$C$
$p^{q_1} * x_1 + C = 0 \space (mod \space p^t)$
仅考虑$x_1$的取值， 方程有解的条件是$C \space | \space gcd(p^{q_1}, p^t)$， 由于$p^{q_1}$是$p^{q_i}$中最小的那个， 所以条件成立， 方程有解， 解数就是$gcd(p^{q_1}, p^t)$即$p^{q_1}$, 然后剩下的x的值可以随便取有$p^t$种取值
所以在$p^t$下的解数是$p^{q_i} * (q^t)^{|A|-1}$, 其中$p^{q_i}$是所有$a_i$中拥有质数$p$的最小指数。
在用中国剩余定理合并一下答案， 可得$f(A,m) = gcd(A, m) * m^{|A| - 1}$
现在要求

$$\sum_m \sum_A gcd(A, m) * m^{|A| - 1}$$
然后使用套路1: 这种gcd求和， 先预处理$w_d$， 表示所有gcd是d的倍数情况下的和
要让gcd是d的倍数， 则gcd中的每一个元素都得是d的倍数， 预处理出数组$C_d$表示B数组中有多少个数是d的倍数， 则
$$w_d = \sum_{d | m} \sum_{i = 1} ^ {C_d} m ^ {i - 1}=\sum_{d|m}\frac{1}{m}(\sum_{i=0}^{C_d}m^i - 1)=\sum_{d | m}\frac{1}{m}*((m + 1)^{C_d}-1)$$
接着使用套路2: 容斥系数
然而求出了这个之后， 直接将$w_d$相加很明显不是答案， 一是$w_d$中对于每个gcd贡献都算作1， 实际上应该是gcd的值， 二是$w_d$考虑的是所有gcd是d的倍数的情况， 直接相加有很多重复。 然而我们还是想直接通过$w_d$的某种求和方式来表达我们相求的式子
考虑每个$w_d$前都乘上一个系数$a_d$
我们想要$\sum a_d * w_d$就是答案， 则根据定义， $a_d$需要满足
$$\sum_{d|n} a_d = n$$
然后我们发现这个性质和$\phi(d)$惊人的相似， 于是我就直接使用欧拉函数。
所以最终算的就是
$$\sum_d \sum_{d|m} \phi(d) * \frac{1}{m} * ((m + 1) ^ {C_d} - 1)$$
由于最后要求的是一个异或和， 把对每个m的贡献单独存起来， 最后求一遍异或和即可。

#include <map>
#include <cmath>
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>

#define ll long long
#define ull unsigned ll
#define db double
#define fi first
#define se second
#define pi pair<int, int >
#define ls (x << 1)
#define rs ((x << 1) | 1)
#define mid ((l + r) >> 1)
#define mp(x, y) make_pair((x), (y))
#define pb push_back
#define sqr(x) ((x) * (x))
#define eps 1e-8

using namespace std;

const int N = (int)1e5 + 10;
const int mo = 998244353;

int tot, phi[N], prim[N]; bool ok[N];
void getPhiAndPrime(){
    phi[1] = 1;
    for (int i = 2; i < N; i ++){
        if (!ok[i]){ok[i] = i; prim[++ tot] = i; phi[i] = i - 1;}
        for (int j = 1; j <= tot; j ++){
            if (1ll * i * prim[j] >= N) break;
            ok[i * prim[j]] = 1;
            if (i % prim[j]) phi[i * prim[j]] = phi[i] * (prim[j] - 1);
            else {phi[i * prim[j]] = phi[i] * prim[j]; break;}
        }
    }
}
ll pw(ll x, int k){
    ll ret = 1;
    while (k){
        if (k & 1) ret = ret * x % mo;
        x = x * x % mo;
        k >>= 1;
    }
    return ret;
}

int n, m; ll A[N], B[N], C[N], ans[N];

int main(){
    getPhiAndPrime();

    int T; scanf("%d", &T);
    while (T --){
        memset(A, 0, sizeof(A));
        memset(C, 0, sizeof(C));
        memset(ans, 0, sizeof(ans));

        scanf("%d %d", &n, &m);
        for (int i = 1; i <= n; i ++){
            scanf("%lld", B + i);
            A[B[i]] ++;
        }

        for (int i = 1; i <= m; i ++)
            for (int j = i; j < N; j += i)
                C[i] += A[j];

        for (int i = 1; i <= m; i ++)
            for (int j = i; j <= m; j += i)
                (ans[j] += phi[i] * pw(j, mo - 2) % mo * (pw(j + 1, C[i]) - 1) % mo) %= mo;

        ll ret = 0;
        for (int i = 1; i <= m; i ++) ret ^= ans[i];

        printf("%lld\n", ret);
    }

    return 0;
}