BZOJ4833 [Lydsy1704月赛] 最小公倍佩尔数

2018 年 12 月 26 日发布.

Description

令 $(1+\sqrt2)^n=e_n+f_n\sqrt2$，其中 $e_n, f_n$ 都是整数。

再令 $g_n=\mathop{\rm lcm}_{i=1}^n f_i$。

给定 $n, p$，求 $\left(\sum_{i=1}^n ig_i\right)\bmod p$。

数据组数 $T\leqslant210,n\leqslant10^6,2\leqslant p\leqslant10^9+7$，$p$ 是质数，$f_1\dots f_n$ 在 $\bmod p$ 意义下不为 $0$。

Solution

根据容斥我们可以得到

$$\mathop{\rm lcm}_{i\in S} f_i=\prod_{\emptyset\neq T\subseteq S}\left(\gcd_{i\in T}f_i\right)^{(-1)^{|T|-1}}$$

这实际上是 min-max 容斥

$$\max(S)=\sum_{\emptyset\neq T\subseteq S}(-1)^{|T|-1}\min(T)$$

的应用。

另外，可以发现 $\gcd(f_i,f_j)=f_{\gcd(i,j)}$ （实际上，所有 $f_i=Af_{i-1}+Bf_{i-2},f_0=0,f_1=1$ 且 $\gcd(A,B)=1$ 的序列 $f$ 都满足这个性质）。

于是

$$\begin{aligned} g_n=&\mathop{\rm lcm}_{i=1}^n f_i\\ =&\prod_{\emptyset\neq T\subseteq S}f_{\gcd(T)}^{(-1)^{|T|-1}}\\ =&\prod_{i=1}^nf_i^{\sum_{i\mid t\leqslant n}\mu(t/i)}\\ =&\prod_{i=1}^n\left(\prod_{d\mid i}f_d^{\mu(i/d)}\right) \end{aligned}$$

所以求出所有的 $f$ 之后莫比乌斯反演出 $\prod_{d\mid i}f_i^{\mu(i/d)}$，然后前缀积就可以得到 $g$ 了。

Code

#include <algorithm>
#include <cstdio>
#include <cstring>

typedef long long LL;

const int N = 1000050;

int n, mod, pr[N], cnt;
bool vis[N];
LL h[N], hi[N];

void Sieve() {
  for (int i = 2; i < N; ++i) {
    if (!vis[i]) pr[cnt++] = i;
    for (int j = 0; j < cnt && i * pr[j] < N; ++j) {
      vis[i * pr[j]] = 1;
      if (i % pr[j] == 0) break;
    }
  }
}

LL pow_mod(LL a, LL b) {
  if ((b %= mod - 1) < 0) b += mod - 1;
  LL ans = 1;
  for (a %= mod; b; b >>= 1, a = a * a % mod)
    if (b & 1) ans = ans * a % mod;
  return ans;
}

int main() {
  Sieve();
  int T;
  scanf("%d", &T);
  while (T--) {
    scanf("%d%d", &n, &mod);
    h[1] = hi[1] = 1;
    for (int i = 2; i <= n; ++i)
      hi[i] = pow_mod(h[i] = (h[i - 1] * 2 + h[i - 2]) % mod, mod - 2);
    for (int i = 0; i < cnt && pr[i] <= n; ++i)
      for (int j = n / pr[i]; j > 0; --j) {
        (h[j * pr[i]] *= hi[j]) %= mod;
        (hi[j * pr[i]] *= h[j]) %= mod;
      }
    LL ans = 0, g = 1;
    for (int i = 1; i <= n; ++i) {
      g = g * h[i] % mod;
      ans = (ans + i * g) % mod;
    }
    printf("%lld\n", (ans + mod) % mod);
  }
}