推公式看这里
https://www.cnblogs.com/acjiumeng/p/10743919.html
推出来公式为
∑
t
=
1
(
t
)
3
∗
⌊
n
t
⌋
3
∗
(
⌊
n
t
⌋
+
1
)
2
∗
(
⌊
n
t
⌋
∗
2
+
1
)
12
∑
d
∣
t
ϕ
(
d
)
∗
μ
(
t
d
)
\sum_{t=1}(t)^3*\frac{{\left \lfloor \frac{n}{t} \right \rfloor}^3*({\left \lfloor \frac{n}{t} \right \rfloor}+1)^2*({\left \lfloor \frac{n}{t} \right \rfloor}*2+1)}{12}\sum_{d|t}\phi(d)*\mu(\frac{t}{d})
t=1∑(t)3∗12⌊tn⌋3∗(⌊tn⌋+1)2∗(⌊tn⌋∗2+1)d∣t∑ϕ(d)∗μ(dt)
但是这位博主的代码实现感觉不太好,我综合了群里大佬的各种写法,写出了一下代码
- 求和公式最后出现了12,而模数是1<<30 ,这样不好求,我们发现直接求 s u m ( i 2 ) , s u m ( i 3 ) sum(i^2),sum(i^3) sum(i2),sum(i3)是可以接受的,这样就避免了除法
- 关于杜教筛求解积性函数
bzoj 4804
f ( n ) = ∑ d ∣ t ϕ ( d ) ∗ μ ( t d ) f(n) = \sum_{d|t}\phi(d)*\mu(\frac{t}{d}) f(n)=d∣t∑ϕ(d)∗μ(dt)
f ( p k ) = ϕ ( p k ) ∗ u ( 1 ) + ϕ ( p ( k − 1 ) ) ∗ u ( p ) = ϕ ( p k ) − ϕ ( p ( k − 1 ) ) f(p^k) = \phi(p^k)*u(1)+\phi(p^{(k-1)})*u(p) = \phi(p^k)-\phi(p^{(k-1)}) f(pk)=ϕ(pk)∗u(1)+ϕ(p(k−1))∗u(p)=ϕ(pk)−ϕ(p(k−1))
当 k = = 1 k==1 k==1
f ( p ) = p − 2 f(p) = p-2 f(p)=p−2
当 k > = 2 k >= 2 k>=2
f ( p k ) = p ( k − 2 ) ∗ ( p − 1 ) 2 f(p^k) = p^{(k-2)}*(p-1)^2 f(pk)=p(k−2)∗(p−1)2
假设 n = n ′ ∗ p k n = n'*p^k n=n′∗pk p 为k的最小的素因子k = 1时
f [ n ] = f [ n ] ∗ f [ p ] f[n] = f[n]*f[p] f[n]=f[n]∗f[p]
当 k = = 2 k == 2 k==2
f [ n ] = f [ n / p / p ] ∗ ( p − 1 ) ∗ ( p − 1 ) f[n] = f[n/p/p]*(p-1)*(p-1) f[n]=f[n/p/p]∗(p−1)∗(p−1)
当 k > = 3 k >= 3 k>=3
f [ n ] = f [ n / p ] ∗ p f[n] = f[n/p]*p f[n]=f[n/p]∗p
根据线性筛的性质,一个数一定是被它最小的素因子筛到,于是我们便可得到代码
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = 1<<30;
const int maxn = 1e7+10;
int pow1[maxn],pow2[maxn];//,pow3[maxn];//,pow4[maxn];
// \sum{d|n} phi*u
LL prime[maxn],cnt,f[maxn];
bool vis[maxn];
void init(int n){
// pow2[i] =
for(LL i = 1;i <= n; ++i){
pow1[i] = pow1[i-1]+i,pow1[i] %= mod;
pow2[i] = pow2[i-1]+i*i%mod,pow2[i] %= mod;
}
f[1] = 1;
for(int i = 2;i <= n; ++i){
if(!vis[i]) prime[cnt++] = i,f[i] = i-2;
for(int j = 0;j < cnt; ++j){
if(prime[j] * i > n)break;
vis[prime[j]*i] = 1;
if(i%prime[j] == 0){
if(i/prime[j]%prime[j] == 0)
f[i*prime[j]] = f[i]*prime[j]%mod;
else
f[i*prime[j]] = f[i/prime[j]]*(prime[j]-1)%mod*(prime[j]-1)%mod;
// if(i*prime[j]== 4)
// cout<<f[i*prime[j]]<<" "<<prime[j]<<endl;
break;
}
else
f[i*prime[j]] = f[i]*f[prime[j]]%mod;
}
}
// cout<<f[4]<<endl;
for(int i = 1;i <= n; ++i){
f[i] = 1ll*i*i%mod*i%mod*f[i]%mod;
f[i] = (f[i-1]+f[i])%mod;
}
}
LL solve(LL n){
LL ans = 0;
for(int i = 1;i <= n; ){
LL t = n/i;
LL t2= min(n/t,n);
// n/t*pow2[n/t]*pow3[n/t]
ans = (ans + (f[t2]-f[i-1])%mod*t%mod*pow1[t]%mod*pow2[t]%mod)%mod;
i = t2+1;
}
return ans;
}
int main(void){
init(maxn-1);
// cout<<f[4]<<endl;
int t;cin>>t;
while(t--){
LL n;scanf("%lld",&n);
printf("%lld\n",(solve(n)%mod+mod)%mod);
}
return 0;
}