2019ICPC西安邀请赛 B. Product

题目链接:计蒜客

对于给定的 n, m, p,我们要计算

\prod_{i=1}^n \prod_{j=1}^n \prod_{k=1}^n m^{\gcd(i,j) [k|\gcd(i,j)]} \mod p

首先,我们可以把乘积扔到指数上,变成求和;然后由于 p 是素数,我们求其指数模 p-1

\sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n \gcd(i,j) [k|\gcd(i,j)] \mod p-1

考虑先对 k 求和,则此式可以变成

f(n) = \sum_{k=1}^n k \left(\sum_{i=1}^{\lfloor n/k\rfloor} \sum_{j=1}^{\lfloor n/k\rfloor} \gcd(i,j) \right)

考虑这样的一个函数,并使用莫比乌斯反演、适当的求和变换和狄利克雷卷积的结论

\begin{aligned} g(N) & = \sum_{i=1}^N \sum_{j=1}^N \gcd(i,j) \\ & = \sum_{s=1}^N s \sum_{i=1}^N \sum_{j=1}^N [\gcd(i,j) = s] \\ & = \sum_{s=1}^N s \sum_{t=1}^{\lfloor N/s \rfloor} \left\lfloor\frac{N}{st}\right\rfloor^2 \mu(t) \\ & = \sum_{d=1}^N \left\lfloor\frac{N}{d}\right\rfloor^2 \sum_{st=d} \mu(t) s \\ & = \sum_{d=1}^N \left\lfloor\frac{N}{d}\right\rfloor^2 \varphi(d) \end{aligned}

回到原式,再使用求和变换

\begin{aligned} f(n) & = \sum_{k=1}^n k \sum_{t=1}^{\lfloor n/k \rfloor} \left\lfloor\frac{n}{kt}\right\rfloor^2 \varphi(k) \\ & = \sum_{d=1}^n \left\lfloor\frac{n}{d}\right\rfloor^2 \sum_{kt=d} \varphi(k)t \end{aligned}

如果后面那个 \varphi * id 可以用杜教筛计算出来,那么再预处理部分前缀和,我们就能 O(n^{2/3}) 计算出结果啦!

既然杜教筛要配一个函数使得其狄利克雷卷积好计算,那么…我们好像只能配上 1 了。

(1 * \varphi * id)(n) = (1 * \mu * id * id)(n) = (id * id)(n) = n \times d(n)

好吧这个式子好像不太好求前缀和?不过考虑一下 d(n) 求前缀和的方法,我们好像还是可以求的。

\sum_{t=1}^n t d(t) = \sum_{t=1}^n \sum_{sk=t} sk = \sum_{s=1}^n s \sum_{k=1}^{\lfloor n/s \rfloor}k

行吧,三个数论分块,写吧…馍馍片。

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;
const size_t MAXN = 1e6+5;
bool isnp[MAXN];
int pcn, pri[MAXN];
lld g[MAXN], q[MAXN], s[MAXN];
lld _P1[MAXN];
map<lld,lld> _P2;
int mod;

void init_prime()
{
    _P1[1] = g[1] = q[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i])
        {
            g[i] = s[i] = 2*i-1;
            q[i] = i-1;
            pri[pcn++] = i;
        }

        for (int j = 0; j < pcn; j++)
        {
            int k = i * pri[j];
            if (k >= MAXN) break;
            isnp[k] = true;

            if (i % pri[j] == 0)
            {
                s[k] = (s[i] + q[i] + 0ll) * pri[j];
                q[k] = q[i] * pri[j];
                g[k] = g[i] / s[i] * s[k];
                break;
            }
            else
            {
                q[k] = q[pri[j]];
                s[k] = s[pri[j]];
                g[k] = g[i] * s[k];
            }
        }

        _P1[i] = _P1[i-1] + g[i];
    }
}

inline lld sumn(lld n)
{
    return n*(n+1)/2%mod;
}

lld sum_xdx(lld n)
{
    lld ans = 0;

    for (lld l = 1, r; l <= n; l = r+1)
    {
        r = n / (n / l);
        (ans += (sumn(r) - sumn(l-1)) * sumn(n/l) % mod) %= mod;
    }

    return ans;
}

lld G(lld n)
{
    if (n < MAXN) return _P1[n] % mod;
    else if (_P2.count(n)) return _P2[n];

    lld ans = sum_xdx(n);
    for (lld k = 2, fk; k <= n; k = fk+1)
    {
        fk = n / (n/k);
        ans -= (fk-k+1) * G(n/k) % mod;
    }

    ans = (ans % mod + mod) % mod;
    return _P2[n] = ans;
}

lld fpow(lld b, lld n, lld MOD)
{
    lld a = 1;

    while (n)
    {
        if (n & 1) a = a * b % MOD;
        b = b * b % MOD;
        n >>= 1;
    }

    return a;
}

inline lld sumn2(lld n)
{
    __int128 s = n*2+1;
    s *= n * (n+1);
    s /= 6;
    return lld(s % mod);
}

int main()
{
    lld n, m, p;
    scanf("%lld %lld %lld", &n, &m, &p);
    mod = p-1;
    init_prime();
    G(n);

    lld fff = 0;
    for (lld l = 1, r; l <= n; l = r+1)
    {
        r = n / (n / l);
        (fff += (G(r) - G(l-1)) * ((n/l) * (n/l) % mod)) %= mod;
    }

    if (fff < 0) fff += mod;
    printf("%lld\n", fpow(m, fff, p));
    return 0;
}

上面式子中,为了求某函数前缀和,写了一堆数组,在这里解释一下:

首先考虑 g(p^n),则有

\begin{aligned} g(p^n) & = \sum_{ij=p^n} \varphi(i) j = \sum_{i=0}^n \varphi(p^i) p^{n-i} \\ & = p^n + n\varphi(p^n) = (n+1)p^n – np^{n-1} \end{aligned}

于是不停的利用最小因子和它的积性去线性筛就好啦~

  • g[n] 表示 g(n)
  • q[n] 表示 n 最小因子 p 以及其对应指数 e 所对应的 (e+1)p^e – e p^{e-1}
  • s[n] 表示 n 最小因子 p 以及其对应指数 e 所对应的 g(p^e)

好麻烦啊QAQ

发表评论

电子邮件地址不会被公开。 必填项已用*标注