【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，那么最后答案就是N2−所求
那么我们可以首先写出比较暴力的东西

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

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

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

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

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

现在考虑分块。
针对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);}