洛谷p3768 简单的数学题

题目大意

求

$$\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)$$ 对质数p取膜的结果

n<=10^10 5*10^8<=p<=1.1*10^9

传送门https://www.luogu.org/problem/show?pid=3768

分析

看到这种两层求和带gcd的式子 可以考虑直接用莫比乌斯反演求解。
但注意到i j的上界都是n，所以用$\phi$推导公式也可以解决这题。(感觉可能会简便一些)

$$\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)$$
$$=\sum_{d=1}^{n}d^3\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}ij[i⊥j]$$
$$=\sum_{d=1}^{n}d^3G(\lfloor \frac{n}{d}\rfloor)$$
右边的式子是关于 $\lfloor \frac{n}{d}\rfloor$的 可以看成一个函数
如果能够快速计算出G函数 那么这个式子可以分块快速求出
$$G(n)=\sum_{i=1}^{n}\sum_{j=1}^{n}ij[i⊥j]$$
$$=\sum_{i=1}^{n}i\sum_{j=1}^{n}j[i⊥j]$$
$$=2\sum_{i=1}^{n}i\sum_{j=1}^{i}j[i⊥j]-1$$
$$=2\sum_{i=1}^{n}i (\frac{i\phi(i)+[i==1]}{2})-1$$
$$=\sum_{i=1}^{n}(i^2\phi(i)+[i==1])-1$$
$$=\sum_{i=1}^{n}i^2\phi(i)$$
(嗯 看起来很简单吧..)

G(n)显然是一个积性函数的前缀和
如果预处理出G函数 那么可以设计一个O(n)的算法 可以通过60%的数据。

记$G(n)=\sum_{i=1}^{n}f(i)$
考虑将$\phi$通过狄利克雷卷积消去
已知$\phi*Id0=Id1$

假设 $f*h=T$
$T(n)=\sum_{d|n}f(d)h(\frac{n}{d})=\sum_{d|n}d^2\phi(d)h(\frac{n}{d})=\sum_{d|n}(d^2h(\frac{n}{d}))\phi(d)$
当$h=Id2$时 $d^2可以与h(\frac{n}{d})配对成n^2 后面的\sum phi(d)也可以合并成一个n$
卷积即可得到$T=Id3$

那么就可以用杜教筛快速计算出G(n)了

$$G(n)=n^3-\sum_{i=2}^{n}i^2G(\lfloor \frac{n}{i}\rfloor)$$
总时间复杂度 $O(n^\frac{2}{3})$

#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
#include<map>
#include<set>
#include<algorithm>
#include<queue>
#include<time.h>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define fe(x) for(int ii=be[x];ii;ii=e[ii].ne)
using namespace std;
const int mn=20000000;
int p,pri[mn/10],pt,mm;
long long n,p2,p6,ans,g[mn+20],mg[100000];
bool vis[100000],isp[mn+20];
inline long long qpow(int a,int n){
    long long x=a,ans=1;
    while (n){
        if (n&1) ans=ans*x%p;
        x=x*x%p;
        n>>=1;
    }
    return ans;
}
inline long long sum2(long long n){
    n%=p;
    return (n+n+1)%p*(n+1)%p*n%p*p6%p;
}
inline long long sum3(long long n){
    n%=p;
    n=n*(n+1)%p*p2%p;
    return n*n%p;
}
inline long long G(long long x){
    if (x<=mm) return g[x];
    int ms=n/x;
    if (vis[ms]) return mg[ms];
    vis[ms]=1;
    long long i=2,j,te;
    long long &ans=mg[ms];
    ans=sum3(x);
    while (i<=x){
        te=x/i;
        j=x/te;
        ans-=(sum2(j)-sum2(i-1))*G(te)%p;
        i=j+1;
    }
    return ans%=p;
}
int main(){
    scanf("%d%lld",&p,&n);
    mm=mn;if (n<mm) mm=n;
    g[1]=1;
    fo(i,2,mm){
        if (!isp[i]){
            g[i]=i-1;
            pri[++pt]=i;
        }
        fo(j,1,pt){
            long long te=i*(long long)pri[j];
            if (te>mn) break;
            isp[te]=1;
            if (i%pri[j]==0){
                g[te]=g[i]*pri[j];
                break;
            }
            g[te]=g[i]*(pri[j]-1);
        }
    }
    fo(i,2,mm) g[i]=(g[i-1]+g[i]*i%p*i)%p;
    p2=qpow(2,p-2);p6=qpow(6,p-2);
    long long i=1,j,te;
    while (i<=n){
        te=n/i;
        j=n/te;
        ans+=(sum3(j)-sum3(i-1))*G(te)%p;
        i=j+1;
    }
    ans=(ans%p+p)%p;
    printf("%lld\n",ans);
    return 0;
}