题意
设 (f(n)=sum_{d|n}{|mu(d)|}),求 (sum_{i=1}^{m}{f(ni)} mod 1e9+7)。
(1leq n,m leq 10^7)
分析
设 (omega(n)) 表示 (n) 的不同质因子个数,根据 (mu) 的定义,有:
[f(n)=sum_{k=0}^{omega(n)}{(omega(n),k)·|(-1)^k|}=2^{omega(n)}
]
显然 (f(n)) 是一个积性函数,且有(f(ab)=frac{f(a)f(b)}{f(gcd(a,b))})。
公式推导
[egin{align}
sum_{i=1}^{m}{f(ni)} & =f(n)sum_{i=1}^{m}{frac{f(i)}{f(gcd(i,n))}}\
& = f(n)sum_{d|n}{frac{1}{f(d)}}sum_{i=1}^{m}{f(i)·[gcd(n,i)==d]}\
& = f(n)sum_{d|n}{frac{1}{f(d)}}sum_{i=1}^{lfloor m/d
floor}{f(id)·[gcd(frac{n}{d},i)==1]}\
& = f(n)sum_{d|n}{frac{1}{f(d)}}sum_{i=1}^{lfloor m/d
floor}{f(id)·sum_{p|i,pd|n}{mu(p)}}\
& = f(n)sum_{d|n}{frac{1}{f(d)}sum_{dp|n}{mu(p)sum_{i=1}^{lfloor m/(dp)
floor}{f(idp)}}}\
& = f(n)sum_{d|n}{sum_{dp|n}{frac{mu(p)}{f(d)}sum_{i=1}^{lfloor m/(dp)
floor}{f(idp)}}}\
end{align}
]
令 (T=dp),有:
[egin{align}
& = f(n)sum_{T|n}{sum_{d|T}{frac{mu(T/d)}{f(d)} sum_{i=1}^{lfloor m/T
floor}{f(iT)}}}\
& = f(n)sum_{T|n}{g(T)} sum_{i=1}^{lfloor m/T
floor}{f(iT)}\
& = sum_{T|n}{f(n)g(T)·sum_{i=1}^{lfloor m/T
floor}{f(iT)}}
end{align}
]
其中 (g(T)=sum_{d|T}{frac{mu(T/d)}{f(d)}}),即 (g=frac{1}{f}mu) ,也是一个积性函数,用线性筛的方法推导可得:
[g(T)=
egin{cases}
left(-frac{1}{2}
ight)^{omega(T)}&, T=prod_{i=1}^{omega(T)}{p_i}\
\
0&, ext{others}\
end{cases}
]
由于只有当 (T) 为 (n) 的各质因子最多出现一次时有效,所以,可以用二进制枚举 (n) 的质因子,(O(8·2^8)) 时间内处理出来。
计算内存限制,可开的 ( ext{int}) 型的数组的大小为 (6e7) 。无法预处理出 (sum_{i=1}^{lfloor m/T floor}{f(iT)}) 。因此,采用离线处理,把所有 (sum_{i=1}^{lfloor m/T floor}{f(iT)}) 的查询,按 (T) 为第一关键字,(lfloor frac{m}{T} floor ·T) 为第二关键字排序,然后直接计算。代码实现上,要注意优化,否则容易超时。
时间复杂度:(O(nlogn+T(sqrt{n}+8*2^8)))
代码
#include <bits/stdc++.h>
#define pb push_back
using namespace std;
const int mod=1e9+7;
const int N=1e7+7;
const int M=1e4+5;
const int inv2=500000004;
typedef long long ll;
struct node
{
int T,u,id;
ll w;
bool operator <(const node b) const
{
if(T==b.T) return u<b.u;
else return T<b.T;
}
};
vector<node>num;
ll ans[M];
int prime[N],cnt;
bool vis[N];
int factor[10],f[N],g[N];
void init()
{//线性筛
int maxn=1e7;
f[1]=g[1]=1;//注意
for(int i=2;i<=maxn;i++)
{
if(!vis[i])
{
prime[++cnt]=i;
f[i]=2;
g[i]=-inv2+mod;
}
for(int j=1;j<=cnt&&prime[j]*i<=maxn;j++)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
f[i*prime[j]]=f[i];
g[i*prime[j]]=0;
break;
}
else
{
f[i*prime[j]]=1LL*f[i]*2%mod;
g[i*prime[j]]=1LL*g[i]*g[prime[j]]%mod;//积性函数
}
}
}
}
void divide(int x,int &cot)
{
cot=0;
for(int i=2;i*i<=x;i++)
{
if(x%i==0)
{
factor[++cot]=i;
while(x%i==0)
x/=i;
}
}
if(x>1) factor[++cot]=x;
}
int main()
{
int t,n,m,tol=0,cot;
scanf("%d",&t);
init();
while(t--)
{
tol++;
scanf("%d%d",&n,&m);
divide(n,cot);
for(int i=0;i<(1<<cot);i++)//二进制枚举质因子,求T
{
int d=1;
for(int j=0;j<cot;j++)
{
if((i>>j)&1) d*=factor[j+1];
}
if(d>m||g[d]==0) continue;//优化
num.pb(node{d,(m/d)*d,tol,1LL*f[n]*g[d]%mod});
}
}
sort(num.begin(),num.end());
ll w=0;
int d=0;
for(int i=0,j=0;i<num.size();i++)
{//注意优化,否则超时
if(num[i].T!=d)
d=num[i].T,w=0,j=num[i].T;
for(;j<=num[i].u;j+=d)
w=(w+f[j])%mod;
ans[num[i].id]=(ans[num[i].id]+num[i].w*w%mod)%mod;
}
for(int i=1;i<=tol;i++)
printf("%lld
",ans[i]);
return 0;
}
参考博客:https://blog.csdn.net/BeNoble_/article/details/108088288