Product（2019年西安邀请赛B+杜教筛）

题目链接

传送门

题面

思路：

我们先将公式进行化简：
$$\begin{array}{rlr}& \prod _{i=1}^n\prod _{i=1}^n\prod _{i=1}^n{m}^{gcd(i,j)[k|gcd(i,j)]}& \\ =& m^{\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^{n}gcd(i,j)[k|gcd(i,j)]}& \end{array}$$
现在我们单独考虑指数部分：
$$\begin{array}{rlr}& \sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^{n}gcd(i,j)[k|gcd(i,j)]& \\ =& \sum _{d=1}^{n}df(d)\sum _{i=1}^n\sum _{k=1}^n[gcd(i,j)=d]& \\ =& \sum _{d=1}^{n}df(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i,j)=1]& \\ =& \sum _{d=1}^{n}df(d)(2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\varphi (i)-1)& \\ =& 2\sum _{d=1}^{n}df(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\varphi (i)-\sum _{d=1}^{n}df(d),\text{其中}f(d)\text{为}d的因子个数& \end{array}$$
对于$\sum _{d=1}^{n}df(d)$我们发现可以进行如下化简求解:
$$\begin{array}{rlr}& \sum _{i=1}^{n}if(i)& \\ =& \sum _{i=1}^{n}i\sum _{j|i}1& \\ =& \sum _{j=1}^n\sum _{j|i}i& \\ =& \sum _{j=1}^n\frac{j\times (1+\lfloor \frac{n}{j}\rfloor )\lfloor \frac{n}{j}\rfloor }{2}& \end{array}$$
对于$\sum _{d=1}^{n}df(d)$我们可以直接通过数论分块来求，对于$\sum _{d=1}^{n}df(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\varphi (i)$我们可以对$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\varphi (i)$进行数论分块用杜教筛求出来，然后在同一块中对$\sum _{d=1}^{n}df(d)$进行数论分块求出来。

代码实现如下

#include <set>
#include <map>
#include <deque>
#include <queue>
#include <stack>
#include <cmath>
#include <ctime>
#include <bitset>
#include <cstdio>
#include <string>
#include <vector>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <unordered_map>
using namespace std;

typedef long long LL;
typedef pair<LL, LL> pLL;
typedef pair<LL, int> pLi;
typedef pair<int, LL> pil;;
typedef pair<int, int> pii;
typedef unsigned long long uLL;

#define lson rt<<1
#define rson rt<<1|1
#define lowbit(x) x&(-x)
#define name2str(name) (#name)
#define bug printf("*********\n")
#define debug(x) cout<<#x"=["<<x<<"]" <<endl
#define FIN freopen("D://code//in.txt","r",stdin)
#define IO ios::sync_with_stdio(false),cin.tie(0)

const double eps = 1e-8;
const int mod = 1000000007;
const int maxn = 5e6 + 3;
const double pi = acos(-1);
const int inf = 0x3f3f3f3f;
const LL INF = 0x3f3f3f3f3f3f3f3fLL;

bool v[maxn];
int n, m, p, cnt, MOD;
int isp[maxn/10], phi[maxn];
int sum[maxn], ssum[maxn];
unordered_map<int, int> dp1, dp2;

int add(int x, int y) {
    if(y < 0) x += y;
    else x = x - MOD + y;
    if(x < 0) x += MOD;
    return x;
}

void init() {
    phi[1] = ssum[1] = 1;
    for(int i = 2; i < maxn; ++i) {
        if(!v[i]) {
            ssum[i] = 2;
            phi[i] = i - 1;
            isp[cnt++] = i;
        }
        for(int j = 0; j < cnt && i * isp[j] < maxn; ++j) {
            v[i*isp[j]] = 1;
            if(i % isp[j] == 0) {
                phi[i*isp[j]] = phi[i] * isp[j];
                int pp = 1, x = i;
                while(x % isp[j] == 0) x /= isp[j], pp++;
                ssum[i*isp[j]] = ssum[i] / pp * (pp + 1);
                break;
            }
            phi[i*isp[j]] = phi[i] * (isp[j] - 1);
            ssum[i*isp[j]] = ssum[i] * 2;
        }
    }
    for(int i = 1; i < maxn; ++i) {
        sum[i] = add(sum[i-1], phi[i]);
        ssum[i] = add(1LL * ssum[i] * i % MOD, ssum[i-1]);
    }
}

int getphi(int x) {
    if(x < maxn) return sum[x];
    if(dp1[x]) return dp1[x];
    int ans = (1LL * x * (x + 1) / 2) % MOD;
    for(int l = 2, r; l <= x; l = r + 1) {
        r = x / (x / l);
        int tmp = 1LL * (r - l + 1) * getphi(x/l) % MOD;
        ans = add(ans, -tmp);
    }
    return dp1[x] = ans;
}

int qpow(int x, int n) {
    int res = 1;
    while(n) {
        if(n & 1) res = 1LL * res * x % p;
        x = 1LL * x * x % p;
        n >>= 1;
    }
    return res;
}

int getnum(int x) {
    if(x < maxn) return ssum[x];
    if(dp2[x]) return dp2[x];
    int ans = 0;
    for(int l = 1, r; l <= x; l = r + 1) {
        r = x / (x / l);
        LL tmp1 = 1LL * ((x / l) + 1) * (x / l);
        LL tmp2 = 1LL * (l + r) * (r - l + 1) / 2;
        ans = add(ans, tmp1 * tmp2 / 2 % MOD);
    }
    return dp2[x] = ans;
}

int main(){
    scanf("%d%d%d", &n, &m, &p);
    MOD = p - 1;
    init();
    int ans2 = getnum(n);
    int ans1 = 0;
    for(int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        int tmp1 = getphi(n / l);
        int tmp2 = add(getnum(r), -getnum(l-1));
        ans1 = add(ans1, 1LL * tmp1 * tmp2 % MOD);
    }
    ans1 = 2LL * ans1 % MOD;
    int ans = add(ans1, -ans2);
    printf("%d\n", qpow(m, ans));
    return 0;
}