【GDOI2017模拟11.7】太阳神

题目描述

太阳神拉很喜欢最小公倍数，有一天他想到了一个关于最小公倍数的题目。
求满足如下条件的数对(a,b)对数:a,b均为正整数且a,b<=n而lcm(a,b)>n。其中的lcm当然表示最小公倍数。答案对1,000,000,007取模
对于20%的数据n<=2000；
对于40%的数据n<=10000000；
对于60%的数据n<=100000000；
对于80%的数据n<=1000000000；
对于100%的数据n<=10000000000。

分析

（中途有些上下界写得不太准确）
20分暴力。
然后我们来看，先把问题反过来，求多少个LCM(A,B)<=N，那么最后答案就是$N^2-所求$
那么我们可以首先写出比较暴力的东西

其中d表示最大公约数，ij就是枚举两个数看。这里ij都是在原本基础上除了d的。
怎么优化呢？
其实中间这一块可以反演，那么可以消除一个条件。

这里为什么n/D^2呢，因为ij枚举的时候已经除了D嘛。
好了可以发现d、D可以交换一下，不影响的。

然后可以把j化掉，毕竟除一下就出来了。

这个东西分块来算的话只能过80，怎么办呢。
我们把字母弄少一点观察一下。

哦那个x就是M
也就是说确定了D,d之后，后面那个东西就是那种形式。
仔细观察发现，它就是叫你求1~M的约数个数和。
这个在$10^7$内的可以预处理。因为n/D是衰减的特别快的，那么当M大于$10^7$时暴力做是没有问题的。
那么最后整理一下上界。

现在考虑分块。
针对D是不行的了，毕竟本身只有根号个。
针对d：让相等的成为一块
针对超过预处理范围的i：相等的成为一块
那么完美解决了。
当然还有别的方法，别的比较好的思路是令d

代码

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
#define fo(i,j,k) for(i=j;i<=k;i++)
const ll mo=1000000007;
const int N=10000000;
ll D,di,dj,m,x,i,j,itmp,y,sqn,dtmp,n,ans,is,dtt;
int miu[N+5],pri[N/5],cur[N+5],t;
ll divs[N+5];
bool pd[N+5];
void work()
{
    fo(D,1,sqn)
    {
        m=n/D/D;
        di=1;
        dtt=dtmp=0;
        while (di<=m)
        {
            x=m/di;
            dj=m/x;
            if (x<=N)
                dtmp=(dtmp+divs[x]*(dj-di+1)%mo)%mo;
            else
            {
                i=1;
                itmp=0;
                while (i<x)
                {
                    y=x/i;
                    j=x/y;
                    itmp=(itmp+(j-i+1)*y%mo)%mo;
                    i=j+1;
                }
                dtmp=(dtmp+itmp*(dj-di+1)%mo)%mo;
            }
            di=dj+1;
        }
        ans=((ll)ans+miu[D]*dtmp+mo)%mo;
    }
}
void predo()
{
    miu[1]=1;
    divs[1]=1;
    int i;
    fo(i,2,N)
    {
        if (!pd[i])
        {
            pri[++pri[0]]=i;
            miu[i]=-1;
            divs[i]=1;
            cur[i]=1;
        }
        fo(j,1,pri[0])
        {
            if (i>N/pri[j]) break;
            t=i*pri[j];
            pd[t]=1;
            miu[t]=-miu[i];
            divs[t]=divs[i];
            if (i%pri[j]==0) 
            {
                cur[t]=cur[i]+1;
                miu[t]=0;
                break;
            }
            divs[t]=(ll)divs[t]*(cur[i]+1)%mo;
            cur[t]=1;
        }
    }
    fo(i,2,N) 
    {
        divs[i]=((ll)divs[i]*(cur[i]+1)%mo+divs[i-1])%mo;
    }
}
int main()
{
    freopen("ra.in","r",stdin);
    freopen("ra.out","w",stdout);
    scanf("%lld",&n);
    sqn=trunc(sqrt(n));
    predo();
    work();
    n%=mo;
    printf("%lld",(n*n%mo-ans+mo)%mo);
}