原题题面
Given you n,x,k , find the value of the following formula:
∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)f(gcd(a1,a2,...ax))∗gcd(a1,a2,...ax)\sum_{a_1=1}^{n}\sum_{a_2=1}^{n}...\sum_{a_x=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})f(gcd(a_1,a_2,...a_x))*gcd(a_1,a_2,...a_x)∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)f(gcd(a1,a2,...ax))∗gcd(a1,a2,...ax)
gcd(a1,a2,…,an)gcd(a_1,a_2,…,a_n)gcd(a1,a2,…,an) is the greatest common divisor of a1,a2,...,a−na_1,a_2,...,a-na1,a2,...,a−n.
The function f(x)f(x)f(x) is defined as follows:
If there exists an ingeter kkk (k>1)(k>1)(k>1) , and k2k^2k2 is a divisor of xxx,
then f(x)=0f(x)=0f(x)=0, else f(x)=1f(x)=1f(x)=1.
输入格式
The first line contains three integers t,k,x(1≤t≤104,1≤k≤109,1≤x≤109)t,k,x (1≤t≤10^4,1≤k≤10^9,1≤x≤10^9)t,k,x(1≤t≤104,1≤k≤109,1≤x≤109)
Then ttt test cases follow. Each test case contains an integer nnn (1≤n≤2×105)(1≤n≤2×10^5)(1≤n≤2×105)
输出格式
For each test case, print one integer — the value of the formula.
Because the answer may be very large, please output the answer modulo 109+710^9+7109+7.
输入样例
3 1 3
56
5
20
输出样例
139615686
4017
11554723
题面分析
设gcd(a1,a2,...ax)=dgcd(a_1,a_2,...a_x)=dgcd(a1,a2,...ax)=d,得
∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)f(d)∗d[gcd(a1,a2,...ax)=d]\sum_{a_1=1}^{n}\sum_{a_2=1}^{n}...\sum_{a_x=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})f(d)*d[gcd(a_1,a_2,...a_x)=d]∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)f(d)∗d[gcd(a1,a2,...ax)=d]
枚举ddd,得到
∑d=1nd∗f(d)∗∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)[gcd(a1,a2,...ax)=d]\sum_{d=1}^{n}d*f(d)*\sum_{a_1=1}^{n}\sum_{a_2=1}^{n}...\sum_{a_x=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})[gcd(a_1,a_2,...a_x)=d]∑d=1nd∗f(d)∗∑a1=1n∑a2=1n...∑ax=1n(∏j=1xajk)[gcd(a1,a2,...ax)=d]
∑d=1ndkx+1∗f(d)∗∑a1=1⌊nd⌋∑a2=1⌊nd⌋...∑ax=1⌊nd⌋(∏j=1xajk)[gcd(a1,a2,...ax)=1]\sum_{d=1}^{n}d^{kx+1}*f(d)*\sum_{a_1=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{a_2=1}^{\lfloor \frac{n}{d}\rfloor}...\sum_{a_x=1}^{\lfloor \frac{n}{d}\rfloor}(\prod_{j=1}^{x}a_{j}^{k})[gcd(a_1,a_2,...a_x)=1]∑d=1ndkx+1∗f(d)∗∑a1=1⌊dn⌋∑a2=1⌊dn⌋...∑ax=1⌊dn⌋(∏j=1xajk)[gcd(a1,a2,...ax)=1]
将[gcd(a1,a2,...ax)=1][gcd(a_1,a_2,...a_x)=1][gcd(a1,a2,...ax)=1]化作μ\muμ,得到
∑d=1ndkx+1∗f(d)∗∑a1=1⌊nd⌋∑a2=1⌊nd⌋...∑ax=1⌊nd⌋(∏j=1xajk)∑t∣a1,t∣a2,...t∣axμ(t)\sum_{d=1}^{n}d^{kx+1}*f(d)*\sum_{a_1=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{a_2=1}^{\lfloor \frac{n}{d}\rfloor}...\sum_{a_x=1}^{\lfloor \frac{n}{d}\rfloor}(\prod_{j=1}^{x}a_{j}^{k})\sum_{t|a_1,t|a_2,...t|a_x}\mu(t)∑d=1ndkx+1∗f(d)∗∑a1=1⌊dn⌋∑a2=1⌊dn⌋...∑ax=1⌊dn⌋(∏j=1xajk)∑t∣a1,t∣a2,...t∣axμ(t)
∑d=1ndkx+1∗f(d)∗(∑i=1⌊nd⌋ik)x∗∑t∣a1,t∣a2,...t∣axμ(t)\sum_{d=1}^{n}d^{kx+1}*f(d)*(\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}i^k)^x*\sum_{t|a_1,t|a_2,...t|a_x}\mu(t)∑d=1ndkx+1∗f(d)∗(∑i=1⌊dn⌋ik)x∗∑t∣a1,t∣a2,...t∣axμ(t)
枚举ttt得
∑d=1ndkx+1∗tkx∗f(d)∗∑t=1⌊nd⌋μ(t)(∑i=1⌊ndt⌋ik)x\sum_{d=1}^{n}d^{kx+1}*t^{kx}*f(d)*\sum_{t=1}^{\lfloor \frac{n}{d}\rfloor}\mu(t)(\sum_{i=1}^{\lfloor \frac{n}{dt}\rfloor}i^k)^x∑d=1ndkx+1∗tkx∗f(d)∗∑t=1⌊dn⌋μ(t)(∑i=1⌊dtn⌋ik)x
设T=tdT=tdT=td,得到
∑T=1nTkxdf(d)∗(∑i=1⌊nT⌋ik)x∗∑t∣Tμ(t)\sum_{T=1}^{n}T^{kx}df(d)*(\sum_{i=1}^{\lfloor \frac{n}{T}\rfloor}i^k)^x*\sum_{t|T}\mu(t)∑T=1nTkxdf(d)∗(∑i=1⌊Tn⌋ik)x∗∑t∣Tμ(t)
为了方便计算,我们把右边的t∣Tt|Tt∣T改成d∣Td|Td∣T,得到
∑T=1nTkx(∑i=1⌊nT⌋ik)x∗∑d∣Tμ(d)df(Td)\sum_{T=1}^{n}T^{kx}(\sum_{i=1}^{\lfloor \frac{n}{T}\rfloor}i^k)^x*\sum_{d|T}\mu(d)df(\frac{T}{d})∑T=1nTkx(∑i=1⌊Tn⌋ik)x∗∑d∣Tμ(d)df(dT)
设G(x)=∑d∣Tμ(d)df(Td)G(x)=\sum_{d|T}\mu(d)df(\frac{T}{d})G(x)=∑d∣Tμ(d)df(dT)
首先μ\muμ可以O(n)O(n)O(n)预处理,fff也可以O(nlogn)O(nlogn)O(nlogn)预处理
那么G可以在O(n)O(n)O(n)下预处理
对于(∑i=1⌊nT⌋ik)x(\sum_{i=1}^{\lfloor \frac{n}{T}\rfloor}i^k)^x(∑i=1⌊Tn⌋ik)x,也可以在O(nlogn)O(nlogn)O(nlogn)下处理完
对于每一次查询,利用分块和预处理前缀和,只需要O(n)O(\sqrt{n})O(n)。
故总复杂度为O(nlogn+Tn)O(nlogn+T\sqrt{n})O(nlogn+Tn)
AC代码(2199ms)
#include <bits/stdc++.h>
using namespace std;
const long long maxn=2e5;
const long long mod=1e9+7;
bool check[maxn+10];//访问标记
int prime[maxn+10];//质数
int mu[maxn+10];//mu函数
int f[maxn+10];//f函数
long long powk[maxn+10];
long long G[maxn+10];
long long sumG[maxn+10];
long long n, k, x;
inline long long quick_pow(long long a, long long b)//快速幂
{
long long ans=1, base=a;
while(b!=0)
{
if (b&1)
ans=(long long) ans*base%mod;
base=(long long) base*base%mod;
b>>=1;
}
return ans;
}
void init()
{
f[1]=mu[1]=1;
int tot=0;
for(int i=2; i<=maxn; i++)
{
f[i]=1;
if (!check[i])
{
prime[tot++]=i;
mu[i]=-1;
}
for(int j=0; j<tot; j++)
{
if (i*prime[j]>maxn) break;
check[i*prime[j]]=true;
if (i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else
{
mu[i*prime[j]]=-mu[i];
}
}
}
for(long long d=2; d<=maxn; d++)//处理f
{
for(long long j=d*d; j<=maxn; j+=d*d)
{
f[j]=0;
}
}
G[0]=0;
for(long long d=1; d<=maxn; d++)//处理G
{
for(long long i=d; i<=maxn; i+=d)
{
G[i]=(G[i]+d*f[d]%mod*mu[i/d]%mod+mod)%mod;
}
}
powk[0]=0;
sumG[0]=0;
for(long long i=1; i<=maxn; i++)//处理累加
{
long long z=quick_pow(i, k);
powk[i]=(powk[i-1]+z)%mod;
sumG[i]=(sumG[i-1]+quick_pow(z, x)*G[i]%mod)%mod;
}
}
long long answer()
{
long long ans=0;
for(long long l=1, r; l<=n; l=r+1)
{
r=n/(n/l);
ans=(ans+(sumG[r]-sumG[l-1]+mod)%mod*quick_pow(powk[n/l], x)%mod)%mod;
}
return ans;
}
void solve()
{
int t;
scanf("%d%lld%lld", &t, &k, &x);
init();
while(t--)
{
scanf("%lld", &n);
printf("%lld\n", answer());
}
}
int main()
{
// ios_base::sync_with_stdio(false);
// cin.tie(0);
// cout.tie(0);
#ifdef ACM_LOCAL
freopen("in.txt", "r", stdin);
freopen("out.txt", "w", stdout);
long long test_index_for_debug=1;
char acm_local_for_debug;
while(cin>>acm_local_for_debug)
{
cin.putback(acm_local_for_debug);
if (test_index_for_debug>100)
{
throw runtime_error("Check the stdin!!!");
}
auto start_clock_for_debug=clock();
solve();
auto end_clock_for_debug=clock();
cout<<"\nTest "<<test_index_for_debug<<" successful"<<endl;
cerr<<"Test "<<test_index_for_debug++<<" Run Time: "
<<double(end_clock_for_debug-start_clock_for_debug)/CLOCKS_PER_SEC<<"s"<<endl;
cout<<"--------------------------------------------------"<<endl;
}
#else
solve();
#endif
return 0;
}
DrGilbert 2020.8.7