LOJ#6229. 這是一道簡單的數學題(莫比烏斯反演+杜教篩)

題目連接函數


\(Description\)
\[\sum_{i=1}^n\sum_{j=1}^i\frac{lcm(i,j)}{gcd(i,j)}\]
答案對\(10^9+7\)取模。spa

\(n<=10^9\)code


\(Solution\)
之前作的反演題都是\(j\)枚舉到\(n\),可是如今\(j\)只枚舉到\(i\)就很是難受,考慮怎麼求\(\sum_{i=1}^n\sum_{j=1}^n\frac{lcm(i,j)}{gcd(i,j)}\)
能夠把它當作是一個\(n*n\)的網格,第\(i\)行第\(j\)列上的數是\(\frac{lcm(i,j)}{gcd(i,j)}\),須要咱們求的是包括對角線在內的下三角矩陣的權值和。
因此答案爲(全部網格權值之和+對角線上的權值和)/2。
\[\sum_{i=1}^n\sum_{j=1}^n\frac{lcm(i,j)}{gcd(i,j)}\]
\[=\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^n\frac{ij}{d^2}[gcd(i,j)==d]\]
\[=\sum_{d=1}^n\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij[gcd(i,j)==1]\]
考慮怎麼求後半部分
\[\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)==1]\]
\[=\sum_{i=1}^n\sum_{j=1}^nij\sum_{d|gcd(i,j)}\mu_d\]
枚舉\(d\)
\[=\sum_{d=1}^n\mu_d\sum_{d|i}^n\sum_{d|j}^nij\]
\[=\sum_{d=1}^n\mu_dd^2\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij\]
\(sum(n)=\sum_{i=1}^ni\)
因此原式
\[=\sum_{i=1}^n\mu_ii^2sum(n/i)^2\]
帶回到一開始的式子裏去
\[\sum_{d=1}^n\sum_{i=1}^{n/d}\mu_ii^2sum(\frac{n}{id})^2\]
按照套路令\(T=id\)
\[=\sum_{T=1}^nsum(n/T)^2\sum_{d|T}\mu_dd^2\]
\(f(x)=\sum_{d|x}\mu_dd^2\),如今若是咱們能夠快速的求出\(f(x)\)的前綴和,那麼就能夠數論分塊算答案了。
但是\(f(x)\)並非一個熟悉的數論函數,怎麼才能用杜教篩呢?
能夠把\(f(x)\)寫成幾個函數的卷積的形式。
\(g(x)=\mu_xx^2\)。那麼\(f=g*1\)。如今要找一個函數\(h\)使得\(f*h=g*1*h\)好算。咱們知道\(\sum_{d|x}\mu_d=e\),因此令\(h(x)=x^2\)來把\(g(x)中的乘x^2\)消掉。
因此就構造出了\(s=f*h=g*1*h=e*1=1\),不難發現\(f\)是個積性函數,能夠線篩。ip

#include<complex>
#include<cstdio>
#include<map>
using namespace std;
const int mod=1e9+7;
const int N=2e6+7;
int n,tot,inv2=mod+1>>1,inv6=166666668;
int prime[N],mu[N],f[N];
bool check[N];
map<int,int>mp;
int qread()
{
    int x=0;
    char ch=getchar();
    while(ch<'0' || ch>'9')ch=getchar();
    while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x;
}
void Init()
{
    int nn=min(n,N-1);
    check[1]=f[1]=1;
    for(int i=2;i<=nn;i++)
    {
        if(!check[i])prime[++tot]=i,f[i]=1-1ll*i*i%mod;
        for(int j=1;j<=tot && i*prime[j]<=nn;j++)
        {
            check[i*prime[j]]=1;
            if(i%prime[j])f[i*prime[j]]=1ll*f[i]*f[prime[j]]%mod;
            else
            {
                f[i*prime[j]]=f[i];
                break;
            }
        }
    }
    for(int i=1;i<=nn;i++)
        f[i]=(f[i-1]+f[i])%mod;
}
int Calc1(int x)
{
    long long res=1ll*x*(x+1)/2%mod;
    return res*res%mod;
}
int Calc2(int x)
{
    return 1ll*x*(x+1)%mod*(x+x+1)%mod*inv6%mod;
}
int Sum(int x)
{
    if(x<N)return f[x];
    if(mp[x])return mp[x];
    long long res=x;
    for(int l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        res=(res-1ll*(Calc2(r)-Calc2(l-1)+mod)*Sum(x/l))%mod;
    }
    return mp[x]=(res+mod)%mod;
}
int main()
{
    scanf("%d",&n);
    Init();
    long long ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=(ans+1ll*Calc1(n/l)*(Sum(r)-Sum(l-1)))%mod;
    }
    printf("%d\n",1ll*(ans+n+mod)*inv2%mod);
    return 0;
}
相關文章
相關標籤/搜索