2019.10.18模拟赛T3

题目大意：

　　求$\sum\limits_{i=1}^n\sum\limits_{j=1}^n[lcm(i,j)>n](n\leq 10^{10})$的值。

题解：

　　这题貌似有n多种做法...

　　为了更好统计，把原式变为$n^2-\sum\limits_{i=1}^n\sum\limits_{j=1}^n[lcm(i,j)\leq n]$。

　　然后开始毒瘤...

　　首先，考虑枚举$lcm(i,j)$，设为$d$，计算有多少对$i.j$的最小公倍数为$d$。

　　设$i=p_1^{a_1}p_2^{a_2}\cdots p_k^{a_k}$，$tp(i)=k$

　　再枚举$gcd(i,j)$，设为$x$，又由于$\frac{i}{gcd(i,j)}$和$\frac{j}{gcd(i,j)}$互质，那么要统计的就是把$\frac{d}{x}$拆成两个互质数的方案数。

　　那么简单想一下，方案数就是$2^{tp(\frac{d}{x})}$，因为同一个质因子不能同时出现在两个数中。

　　于是答案变为：

　　　　$\sum\limits_{d=1}^n\sum\limits_{x|d}2^{tp(x)}$

　　枚举$x$，即$\sum\limits_{x=1}^n2^{tp(x)}\lfloor\frac{n}{x}\rfloor$。

　　然后我们发现$\lfloor\frac{n}{x}\rfloor$可以进行数论分块，所以需要求出$2^{tp(x)}$的前$n$项和。

　　但是...我不会啊！！

　　所以接下来我开始打表，然后我惊奇地发现：

　　　　$2^{tp(x)}=\sum\limits_{t|x}\mu^2(t)$！！！(别问我是怎么发现的)

　　带入原式，得：$\sum\limits_{x=1}^n\lfloor\frac{n}{x}\rfloor\sum\limits_{t|x}\mu^2(t)$

　　反手枚举$t$，即$\sum\limits_{t=1}^n\mu^2(t)\sum\limits_{x=1}^{n/t}\lfloor\frac{n}{tx}\rfloor$

　　另外，$\sum\limits_{i=1}^nd(i)=\sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor$，其中$d(i)$表示$i$的约数个数。

　　因此答案为：$\sum\limits_{t=1}^n\mu^2(t)\sum\limits_{x=1}^{n/t}d(x)$

　　式子结束了，需要调用两层数论分块，还有要预处理出$\mu^2(t)$以及$d(x)$的部分前缀和。

　　时间复杂度...我也不知道，大概在$O(n^{\frac{3}{4}})$左右吧...

　　

#include<bits/stdc++.h>
#define mod 1000000007
#define ll long long
#define N 20000010
using namespace std;
int pr[N],mu[N],t[N],d[N];
ll n,ans;
bitset<N>vis;
void make(int m)
{
    mu[1] = d[1] = 1;
    for(int i = 2;i<=m;i++)
    {
        if(!vis[i])pr[++pr[0]] = i,mu[i] = -1,d[i] = 2;
        for(int j = 1;i*pr[j]<=m;j++)
        {
            vis[i*pr[j]] = 1;
            d[i*pr[j]] = d[i]<<1;
            if(i%pr[j]==0){d[i*pr[j]]-=d[i/pr[j]];break;}
            mu[i*pr[j]] = -mu[i];
        }
    }for(int i = 1;i<=m;i++)t[i] = t[i-1]+mu[i]*mu[i],d[i] = (d[i]+d[i-1])%mod;
}
ll S(ll m)//mu^2[i]前缀和
{
    if(m<=20000000)return t[m];
    ll ret = 0,tt = (ll)sqrt(m)+1;
    for(ll i = 1;i<=tt;i++)ret+=mu[i]*(m/i/i);
    return ret;
}
ll SS(ll m)//d(i)前缀和
{
    if(m<=20000000)return d[m];
    ll ret = 0,nxt;
    for(ll i = 1;i<=m;i = nxt+1)
    {
        nxt = m/(m/i);
        ret = (ret+(m/i)*(nxt-i+1)%mod)%mod;
    }return ret;
}
int main()
{
    scanf("%lld",&n);
    make(20000000);
    ans = (n%mod)*(n%mod)%mod;
    ll nxt;
    for(ll i = 1;i<=n;i = nxt+1)
    {
        nxt = n/(n/i);
        ans = (ans-(S(nxt)-S(i-1))%mod*SS(n/i)%mod)%mod;
    }printf("%lld\n",(ans+mod)%mod);
    return 0;
}