题意
给定 (n,c,d) 以及 (b) 数组，解方程组：
[egin{equation}
sum_{j = 1}^{n} gcd(i, j)^c cdot lcm(i, j)^d cdot x_j equiv b_i pmod{p}
end{equation}
]
(n le 10^5)
题解
推式子：
[sum_{j=1}^n (i,j)^{c-d}j^dx_j=frac{b_i}{i^d}\
sum_{j=1}^n(i,j)^{c-d}X_j=B_i
]
发现这个 ((i,j)^{c-d}) 很烦人，如果 (c-d=1)，我们或许可以用 (id=varphi * 1) 之类的东西来反演，可是 (c-d) 不一定为 (1)。不过我们可以直接莫比乌斯反演来搞到个什么类似 (d|i,d|j) 的东西来摆脱掉那个 (gcd)：
[egin{aligned}
f(x)&=x^{c-d}\
f(x)&=sum_{d|x}g(d)\
g(x)&=sum_{d|x}f(d)mu(x/d)\
end{aligned}
]
(g) 可以直接算。那么继续推:
[egin{aligned}
sum_{j=1}^nsum_{d|i,d|j}g(d)X_j&=B_i\
sum_{d|i}g(d)sum_{j=1}^{lfloor frac{n}{d}
floor}X_{jd}&=B_i\
sum_{d|i}g(d)S(d)&=B_i
end{aligned}
]
其中 (S(d)) 表示 (sum_{j=1}^{lfloor frac{n}{d} floor}X_{jd})。于是我们可以类似杜教筛那样搞出 (S(d))：
[g(i)S(i)=B_i-sum_{d|i,d
ot= i}g(d)S(d)
]
现在我们能够搞出 (S)，但是没搞出 (X)。不过 (X) 也可以被类似的方法推出来：
[X_d=S(d)-sum_{d|k,d
ot= k}X_k
]
只不过这回是从大往小推。