洛谷3768：简单的数学题

洛谷3768：简单的数学题

题意描述：

• 求([sum_{i=1}^nsum_{j=1}^nijgcd(i,j)] mod p)。
• 数据范围(nleq10^{10},p)为质数（long long）。

思路：

首先这个(n)的数据范围大概率是要杜教筛的节奏啊。

先化式子。

枚举(gcd(i,j)=d)。

• (sum_{d=1}^ndsum_{i=1}^nsum_{j=1}^nij[gcd(i,j)=d]).

提取(d)。

• (sum_{d=1}^nd^3sum_{i=1}^frac{n}{d}sum_{j=1}^frac{n}{d}ij[gcd(i,j)=1]).

我们知道，从(1)到(n)中与(n)互质的数的和为：

• (sum_{i=1}^ni[gcd(i,n)=1]=frac{n*varphi(n)}{2}).

提一下：

• (sum_{d=1}^nd^3sum_{i=1}^frac{n}{d}isum_{j=1}^frac{n}{d}j[gcd(i,j)=1]).

原式就等于：

• (sum_{d=1}^nd^3sum_{i=1}^frac{n}{d}i*i*varphi(i)).

前面(frac{n}{d})有(sqrt{n})种选择，所以可以数论分块一下，后面的算是很明显的杜教筛了。

设(f=varphi idid,S(n)=sum_{i=1}^nf(i))。设(g(i)=idid)。

那么有：

• ((f*g)(n)=sum_{d|n}(varphi(d)d^2)frac{n^2}{d^2}=n^3)。

这玩意前缀和是(frac{n^2(n+1)^2}{4})。

根据杜教筛有：

• (g(1)S(n)=sum_{i=1}^n(f*g)(i)-sum_{i=2}^ng(i)S(frac{n}{i})).

也就是(S(n)=frac{n^2(n+1)^2}{4}-sum_{i=2}^ni^2S(frac{n}{i}))。

原式就是：

• (sum_{d=1}^nd^3S(frac{n}{d})).

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int inv[15];
const int maxn = 5e6;
ll n, mod;

void get_inv(int n){
    inv[1] = 1;
    for(int i = 2; i <= n; i++)
    inv[i] = ((ll)(mod-mod/i)*inv[mod%i])%mod;
}

int primes[maxn+10], cnt;
bool vis[maxn+10];
ll phi[maxn+10];
void init(int n)
{
    phi[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!vis[i])
        {
            primes[++cnt] = i;
            phi[i] = i-1;
        }
        for(int j = 1; primes[j] <= n/i; j++)
        {
            vis[primes[j]*i] = 1;
            if(i % primes[j] == 0)
            {
                phi[primes[j]*i] = phi[i]*primes[j];
                break;
            }
            else phi[i*primes[j]] = phi[i]*(primes[j]-1);
        }
    }
    for(ll i = 1; i <= n; i++)
    phi[i] = (phi[i-1]%mod+i%mod*i%mod*phi[i]%mod)%mod;
}

inline ll get_cube(ll n){
    n %= mod;
    return n%mod*(n+1)%mod*(2*n+1)%mod*inv[6]%mod;
}

unordered_map<ll, ll> Sphi;

ll getSphi(ll n)
{
    if(n <= maxn) return phi[n];
    if(Sphi[n]) return Sphi[n];
    ll res = n%mod*n%mod*(n+1)%mod*(n+1)%mod*inv[4]%mod;
    for(ll l = 2, r; l <= n; l = r+1)
    {
        r = n/(n/l);
        res -= getSphi(n/l)%mod*(((get_cube(r)-get_cube(l-1))%mod)+mod)%mod;
        res = (res%mod+mod)%mod;
    } return Sphi[n] = res;
}

inline ll get_n3(ll n){
    n %= mod;
    return n%mod*n%mod*(n+1)%mod*(n+1)%mod*inv[4]%mod;
}

int main()
{
    cin >> mod >> n; get_inv(10);
    init(maxn);
    ll ans = 0;
    for(ll l = 1, r; l <= n; l = r+1)
    {
        r = n/(n/l);
        ans = (ans+(((get_n3(r)-get_n3(l-1))%mod)+mod)%mod*getSphi(n/l)%mod) % mod;
    }
    cout << ans << endl;
    return 0;
}