标签:input 使用 求和 i+1 100% 就是 namespace cas 更换
去年的Lucas非常喜欢数论题,但是一年以后的Lucas却不那么喜欢了。
在整理以前的试题时,发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”,其中 表示i的约数个数。他现在长大了,题目也变难了。
求如下表达式的值:
其中 表示ij的约数个数。
他发现答案有点大,只需要输出模1000000007的值。
第一行一个整数n。
一行一个整数ans,表示答案模1000000007的值。
2
8
对于100%的数据n <= 10^9。
这个题的最大难点,在于对\(f(ij)\)的变形,只要这个变形正确了,后面的就是更换求和指标经典套路+数论分块经典套路了。
\(f(ij)=\sum_{x|i}\sum_{y|j}[(x,y)=1]\)
证明:ij的某个因子一定是i的某个因子*j的某个因子乘起来的,我们不妨设为i和\(\frac{j}{y}\)的某个因子,那么设\(p=(x,y)\),那么你在x中包括了p这个因子,又在\(\frac{j}{y}\)中把它消掉,就没意义了,也就会重复统计,所以只有\((x,y)=1\)的时候才会有合法的贡献。
然后有了\((x,y)=1\)这个条件,直接上莫比乌斯函数:
\(\sum^{n}_{i=1}\sum_{j=1}^{n}f(i,j)=\sum^{n}_{i=1}\sum_{j=1}^{n}\sum_{x|i}\sum_{y|j}\sum_{d|x,d|y}\mu(d)\)
后面一步就是更换求和指标啦,把d提到最前面,x,y其次,i,j最后面,因为有两部分完全相同,所以这个式子就变成了:
$\sum\limits_{d=1}^n \mu(d)(\sum\limits_{i=1}^{\lfloor {n\over i}\rfloor} {\lfloor {n\over i*d}\rfloor})^2 $
先对于n/i分块,然后对于n/(i/d)分块,前面的莫比乌斯函数根据n/i的分块范围在线使用杜教筛计算。
时间复杂度$O(n^{3/4}+n^{2/3}logn) $
#include <bits/stdc++.h>
#define ll long long
using namespace std;
map<ll,ll>mmp;int sum[1000005],pri[1000005],vis[1000005],tot,mu[1000005],n,ls,ans,P=1000000007;
ll djs(int x)
{
if(x<=1e6) return sum[x];
if(mmp.count(x)) return mmp[x];
int ans=1,ls;
for(int i=2;i<=x;i=ls+1) ls=x/(x/i),ans=(ans-1ll*(ls-i+1)*djs(x/i)%P+P)%P;
return mmp[x]=ans;
}
ll cal(ll x)
{
ll ans=0,ls;
for(int i=1;i<=x;i=ls+1) ls=x/(x/i),ans=(ans+1ll*(ls-i+1)*(x/i)%P)%P;
return 1ll*ans*ans%P;
}
int main()
{
mu[1]=sum[1]=1;
for(int i=2;i<=1000000;i++)
{
if(!vis[i]) pri[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot&&i*pri[j]<=1000000;j++)
{
vis[i*pri[j]]=1;
if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
mu[i*pri[j]]=-mu[i];
}
sum[i]=(sum[i-1]+mu[i]+P)%P;
}
scanf("%d",&n);
for(int i=1;i<=n;i=ls+1) ls=n/(n/i),ans=(ans+1ll*(djs(ls)-djs(i-1)+P)%P*cal(n/i)%P)%P;
printf("%d\n",ans);
}
标签:input 使用 求和 i+1 100% 就是 namespace cas 更换
原文地址:https://www.cnblogs.com/CK6100LGEV2/p/9398608.html