P3768 简单的数学题

· · 题解

这是一道给人启示的一道数学题

但是我太菜了，只会推一下版版的式子

原式：

\sum_{i=1}^n\sum_{j=1}^ni*j*gcd(i,j)

化简：

\sum_{d=1}^nd^3\sum_{i=1}^{[n/d]}\sum_{j=1}^{[n/d]}i*j*[gcd(i,j)==1]

接着：

\sum_{d=1}^nd^3\sum_{i=1}^{[n/d]}\sum_{j=1}^{[n/d]}i*j\sum_{k|gcd(i,j)}\mu(k)

把\mu(k)提前

然后：

\sum_{d=1}d^3\sum_{k=1}^{[n/d]}\mu(k)*k^2\sum_{i=1}^{[n/kd]}\sum_{i=1}^{[n/kd]}i*j

我们记S(k)为(\sum_{i=1}^ki)^2

我们令T=k*d

\sum_{T=1}^nT^2*\sum_{k|T}\mu(k)*id(T/d)*S(T) \sum_{T=1}^nT^2\phi(T)*S(T)

我们知道f(x)=x^2和\phi(x)是积性函数

那么我们考虑狄利克雷卷积

令h(x)=T^2*\phi(T),我们现在需要一个适合的函数g(x)

\sum_{i=1}^n(h*g)(i) \sum_{i=1}^n\sum_{d|i}h(i)*g(i/d)

于是：

\sum_{i=1}^n\sum_{T|i}^nT^2*\phi(T)*g(i/T)

那么我们需要找到一个合适的g

于是猛然发现函数f(x)=x^2是个不错的选择

\sum_{i=1}\sum_{T|i}\phi(i)*i^2

然而：

\sum_{i|n}\phi(i)=n

∴

\sum_{i=1}^n(f*g)(i)=\sum_{i=1}i^3

又有不知道从哪里知道的：

\sum_{i=1}^ni^3=S(n)

于是我们达到了目的，在O(1)时间内求出卷积

\sum_{i=1}(f*g)(i)

的值，并且快速求出函数

g(x)=x^2

的前缀和

还是不知道是在哪离学的:

\sum_{i=1}^ni^2=n*(n+1)*(2*n+1)/6

其中S(n)是前面提到的

于是我们便可以数论分块再套一个杜教筛解决此题

而我还是太菜了，还是只会推一些死的式子

贴一波蒟蒻的代码，好让大佬嘲讽

模数叫remoon，不怕不A题，

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<map>
using namespace std;
#define int long long
const int N = 5000000 ;
int vis[N+1000] , p[N+1000] , phi[N+1000] , sumphi[N+1000] ;
int num , n , remoon , inv ;
map<int ,int> pphi ;
int pow(int a , int b)
{
    int Ans = 1;
    while(b) {
        if(b & 1) Ans = Ans * a % remoon;
        a = a * a % remoon , b >>= 1;
    }
    return Ans;
}
void find(int x)
{
    sumphi[1] = 1 ;
    phi[1] = 1 ;
    vis[1] = 1 ;
    for(int i = 2 ; i <=  x + 10; i++ ){
        if( !vis[i] ) p[++num] = i , phi[i] = i - 1 ;
        for(int j = 1 ; j <=num ; j++ ){
            if( p[j] * i > x + 10 ) break ;
            vis[ i * p[j] ] = 1 ;
            if( i % p[j] == 0 )
            {
                phi[ p[j] * i ] = phi[i] * p[j] ;
                break ;
            }
            phi[ i * p[j] ] = ( p[j] - 1 ) * phi[i] ;
        }
    }
    for(int i = 2 ; i <= x ; i++ ) ( sumphi[i] = ( sumphi[i - 1] + i % remoon * i % remoon * phi[i] % remoon )) %= remoon ;
    return ;
}
int h(int x)
{
    x %= remoon ;
    return x * (x + 1) % remoon * (x + x + 1) % remoon * inv % remoon ;
}
int s(int x)
{
    x %= remoon ;
    return (x * (x + 1) / 2) % remoon * ((x * (x + 1) / 2) % remoon) % remoon; 
}
int get_phisum(int x)
{
    if(x < N) return sumphi[x];
    if(pphi[x]) return pphi[x];
    int Ans = s(x);
    for(int l =  2 , r ; l <= x ; l = r + 1) {
        r = x / (x / l);
        (Ans -= (h(r) - h(l - 1) + remoon) % remoon * get_phisum(x / l) % remoon) %= remoon ;
        //printf("%lld\n" , Sum_Phi(x / l));
    }
    (Ans += remoon) %= remoon;
    return pphi[x] = Ans;
}
signed main()
{
   // freopen("toutput.out","w",stdout) ;
    scanf("%lld%lld" , &remoon , &n ) ;
    inv = pow( 6 , remoon - 2 ) ;
    find(N + 3) ;
    int Ans = 0 ;
    for(int l = 1 , r ; l <= n ; l = r + 1) {
        r = n / (n / l);
        //printf("%lld %lld %lld\n" , Sum_Phi(r) , Sum_Phi(l - 1) , Sum3(n / l));
        (Ans += ( get_phisum(r) - get_phisum(l - 1) + remoon) % remoon * s(n / l) % remoon) %= remoon ;
    }

    printf("%lld\n" , ( Ans += remoon ) %= remoon ) ;
    return 0 ;
}

最后膜一发zyx：

zyx：我曾经在比赛的时候打了个暴力就过来，后来有人告诉我：这个

其实是可持久化zkw线段树套可持久化平衡树以维护上下凸包并以类似

CDQ分治的方法进行不强制在线的动态快速排序来解决斜率优化加四边

形不等式优化动态dp问题以实现可持久化并查集在对偶图上跑

dijkstra代替原图中用HLPP优化上下界的网络最大流将模数用唯一分

解定理分解后答案用拓展卢卡斯定理求出并用中国剩余定理合并答案

最后使用套欧拉函数和莫比乌斯反演进行一次树上莫队分块再将最后

的答案使用矩阵优化斐波那契数列套快速幂求出还要记得用线性时间

求出逆元。