P3768 简单的数学题

题意描述

洛谷

求 (displaystyle(sum_{i=1}^{n}sum_{j=1}^{n} ij gcd(i,j) )mod p)

数据范围： (nleq 10^{10}, p 为质数)。

solution

推了好久才推出来的一道题。

先枚举一下 (gcd(i,j)) 可得：

(displaystylesum_{d=1}^{n} dsum_{i=1}^{n}sum_{j=1}^{n} ij[gcd(i,j) == d])

后面柿子可以提出来一个 (d) 变为：

(displaystyle sum_{d=1}^{n}d^3 sum_{i=1}^{{nover d}}sum_{j=1}^{nover d} ij[gcd(i,j) == 1])

莫比乌斯反演一下可以变为：

(displaystylesum_{ d=1}^{n} d^3 sum_{i=1}^{nover d}sum_{j=1}^{nover d}sum_{pmid gcd(i,j)} mu(p) ij)

(displaystylesum_{d=1}^{n} d^3sum_{p=1}^{nover d} mu(p)sum_{i=1}^{nover d}sum_{j=1}^{nover d} ij[pmid i,pmid j])

我们观察一下后面的 (displaystylesum_{i=1}^{nover d}sum_{j=1}^{nover d} ij[pmid i,pmid j]) 这个柿子，你会发现：

(i) 的取值可能为： ( p, 2p, 3p, 4p....{nover dp})

(j) 的取值可能为：( p,2p, 3p, 4p...{nover dp})

设 (S(n) = 1+2+3+4....+n), 那么 (displaystylesum_{i=1}^{nover d}sum_{j=1}^{nover d} ij[pmid i,pmid j] = p^2 S({nover dp})^2)。

带入原式可得：

(displaystylesum_{d=1}^{n} d^3 sum_{p=1}^{nover d} mu(p)p^2 S({nover dp})^2)

(dp) 化 (Q) 可得：

(displaystylesum_{Q=1}^{n} S({nover Q})^2 sum_{dmid Q} d^3 ({Qover d})^2 mu({Qover d}))

(= displaystyle sum_{Q=1}^{n} S({nover Q})^2 sum_{dmid Q} Q^2 d mu({Qover d}))

(= displaystylesum_{Q=1}^{n}S({nover Q})^2 Q^2 sum_{dmid Q} d mu({Qover d}))

然后你会发现 (displaystyle sum_{dmid Q} dmu({Qover d})) 他是个卷积的形式，由 (phi = mu * id) 可得： (displaystyle sum_{dmid Q} dmu({Qover d}) = phi(Q)), 代入原式可得：

(displaystylesum_{Q=1}^{n} S({nover Q}) ^2 Q^2 phi(Q))

前面那个可以拿整除分块来解决，设 (f(i) = i^2 phi(i)) ,现在关键是怎么求出来 (f(i)) 的前缀和。

推到这里自己就不会了，看了看题解总算懂了。

(n leq 10^{10}), 所以线性筛又萎了，在这个数据范围下，能做的好像只剩下杜教筛了。

套一下积性函数的柿子 ： (displaystyle h(n) = sum_{dmid n} f(d) g({nover d}))

因为 (f(i)) 里面有 (i^2) ，所以考虑构造 (g({nover i})) 把这一项消掉，那么就 设 (g(i) = i^2) 。

根据小学知识可得：(displaystyle sum_{i=1}^{n} g(i) = {n(n+1)(2n+1)over 2}), (g(i)) 的前缀和就可以 (O(1)) 的算出来了。

把 (g(i)) 代入原式可得：

(displaystyle h(n) = sum_{dmid n} d^2 phi(d) ({nover d})^2)

(displaystyle h(n) = n^2 sum_{dmid n} phi(d))

又因为 (displaystyle sum_{dmid n}phi(d) = n), 代入可得： (h(n) = n^3)

根据小学知识可得： $displaystylesum_{i=1}^n h(i) = ({n imes (n+1)over 2})^2 $, 这可也可以 (O(1)) 的来求。

然后使用杜教筛就可以在非线性的时间内把 (f(i)) 的前缀和算出来了。

最后的答案用数论分块计算一下即可，然后这道题就做完了QWQ。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<map>
using namespace std;
#define int long long
using namespace std;
const int N = 1e7+10;
int p,n,tot,inv2,inv6,ans;
int prime[N],phi[N],f[N];
bool check[N];
map<int,int> vis;
inline int read()
{
    int s = 0,w = 1; char ch = getchar();
    while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
    return s * w; 
}
int ksm(int a,int b)
{
    int res = 1;
    for(; b; b >>= 1)
    {
        if(b & 1) res = res * a % p;
        a = a * a % p;
    }
    return res;
}
void YYCH()
{
    phi[1] = 1;
    for(int i = 2; i <= N-5; i++)
    {
        if(!check[i])
        {
            prime[++tot] = i;
            phi[i] = i-1;
        }
        for(int j = 1; j <= tot && i * prime[j] <= N-5; j++)
        {
            check[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                phi[i * prime[j]] = phi[i] * prime[j] % p;
                break;
            }
            else phi[i * prime[j]] = phi[i] * phi[prime[j]] % p;
        }
    }
    for(int i = 1; i <= N-5; i++) f[i] = i * i % p * phi[i] % p;
    for(int i = 1; i <= N-5; i++) f[i] = (f[i] + f[i-1]) % p;
}
int sqr(int x)
{
	x %= p; 
	return x * (x+1) % p * (2 * x % p + 1) % p * inv6 % p;
}
int Sum(int n)
{
    if(n == 0) return 0;
    n %= p;//这里 n*n 直接算的话会炸long long 所以要预先模一下
    int tmp = n * (n + 1) % p * inv2 % p;
    return tmp * tmp % p;
}
int get(int n)
{
    if(n <= N-5) return f[n];
    if(vis[n]) return vis[n];
    int res = Sum(n);
    for(int l = 2, r; l <= n; l = r+1)
    {
        r = min(n,n/(n/l));
        res -= get(n/l) * ((sqr(r)-sqr(l-1) % p + p) % p) % p;
	res = (res + p) % p; 
    }
    vis[n] = (res % p + p) % p;
    return vis[n];
}
signed main()
{
    p = read(); n = read(); YYCH(); inv2 = ksm(2,p-2); inv6 = ksm(6,p-2);
    for(int l = 1, r; l <= n; l = r+1)
    {
        r = min(n,n/(n/l));
        ans = (ans + Sum(n/l) * ((get(r)-get(l-1) + p) % p) % p) % p;
    }
    printf("%lld
",ans);
    return 0;
}