拉格朗日插值

基础操作with 一个人的数论

首先推一下式子。

我们考虑多项式科技。构造一个多项式 \(F(x)\)，形式如下：

\[\large F(x) = \sum_{i = 1}^{x}i^d [\gcd(i,n)=1] \]

我们进行一些基础的反演：

\[\large\begin{aligned} F(x) = & \sum_{i = 1}^{x}i^d [\gcd(i,x)=1] \\= & \sum_{i = 1}^{x}i^d \sum_{e|\gcd(i,x)} \mu(e) \\= & \sum_{e|x} \mu(e) \sum_{e|i}^{x}i^d \\= & \sum_{e|x} \mu(e) \sum_{i = 1}^{ \frac x e}(ie)^d \\= & \sum_{e|x} \mu(e) e^d \sum_{i = 1}^{\frac x e }i^d \end{aligned} \]

可以发现，我们如果用另一个多项式 \(G(x) = \sum\limits_{i = 1}^x i^d\) 表示后面的部分，那么式子化成

\[\large F(x) = \sum_{e|n} \mu(e) e^d G(\frac x e) \]

可以发现，\(G(x)\) 为一个关于 \(x\) 的 \(d+1\) 次多项式，于是可以换一下表示方法。

\[\large G(x) = \sum_{i = 0}^{d+1} f_i x^i \]

将它带入 \(F(x)\) 的式子，并化简，有如下过程：

\[\large\begin{aligned} F(x) = & \sum_{e|x} \mu(e) e^d \sum_{i = 1}^{\frac x e }i^d \\= & \ \sum_{e|x} \mu(e) e^d G(\frac x e) \\= & \ \sum_{e|x} \mu(e) e^d \sum_{i = 0}^{d+1} f_i (\frac x e)^i \\= & \ \sum_{e|x} \mu(e) e^{d-i} \sum_{i = 0}^{d+1} f_i x^i \\= & \ \sum_{i = 0}^{d+1} f_i x^i \sum_{e|x} \mu(e) e^{d-i} \end{aligned} \]

我们把 \(e\) 给整数唯一分解掉，用以下式子表示：$\large \prod_{p_i | e} p_i^{\alpha_i} $
于是有

\[\large\begin{aligned} F(x) = & \sum_{i = 0}^{d+1} f_i x^i \sum_{e|x} \mu(e) e^{d-i} \\= & \ \sum_{i = 0}^{d+1} f_i x^i \prod_{p_j | e} \sum_{t = 0}^{\alpha_i} \mu(p_j ^ t) p_j ^ {t (d-i)} \end{aligned} \]

看到 \(\mu\) 函数里出现了平方，我们考虑它的取值。容易发现以下取值规律：

\[\large \mu(p^\alpha)=\left\{ \begin{aligned} & 1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \alpha=0 \\ & -1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \alpha=1 \\ & 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \alpha>1 \\ \end{aligned} \right. \]

回到原来的式子。我们可以发现 \(t\) 只在取 0 与 1 时有意义，因此可以拆开：

\[\large\begin{aligned} F(x) = & \sum_{i = 0}^{d+1} f_i x^i \prod_{p_j | e} \sum_{t = 0}^{\alpha_i} \mu(p_j ^ t) p_j ^ {t (d-i)} \\= & \ \sum_{i = 0}^{d+1} f_i x^i \prod_{p_j | e} (1 - p_j ^ {d-i}) \end{aligned} \]

如果我们能把 \(f_i\) 搞出来，那我们就能每次带入一个质因子计算最终结果了。考虑拉格朗日插值。

我们考虑带入 \([1, k+3]\) 这 \(k+3\) 个数值为 \(x\)，求得共 \(k+2\) 个 \(f_i\)。同时，由于带入值的特殊性，我们有$ x_i = i$。代码如下：

inline int qp(int a, int b = mod-2) {
    int ans = 1, bse = a;
    while (b) {
        if (b & 1) ans = 1ll * ans * bse % mod;
        bse = 1ll * bse * bse % mod;
        b >>= 1;
    } return ans;
}

inline void Poly_Mul(int a, int b) {
	// 多项式乘法，每次向A里面乘一个 (ax + b)
	++len;
	pre(i,len,1) {
		A[i] = (1ll * A[i] * b + A[i-1] * a) % mod;
	} A[0] = 1ll * A[0] * b % mod;
}

inline void init(int k) {
	// 用f的定义插出fi
	A[0] = 1, len = 0;
	int s = 0, mul;
	rep(i,1,k+3) {
		s = (s + qp(i, k)) % mod; // 这是 f 的前缀和,现在循环求出的是 F(i)
		mul = qp(s); // 1 / f(n)，作分母
		// 我们带入的是 x = 1~k+3
		// 所以 xi = i
		rep(j,1,k+3) {
			if (i == j) continue;
			Poly_Mul(1, mod-j); // 分子，乘进去 (x - x_j) = (x - j)
			mul = 1ll * mul * (i+mod-j) % mod; // 分母乘上 (x_i - x_j) = (i - j)
		} mul = qp(mul);
		rep(p,0,len) A[p] = 1ll * A[p] * mul % mod;// 再把分母乘进去，我们就直接插出了系数
		// 每次系数相加作答案
		rep(p,0,len) f[p] = (f[p] + A[p]) % mod, A[p] = 0;
		len = 0, A[0] = 1;
	}
}