问题为:$$ans=\sum^{n}_{i=1}\sum^{n}_{j=1}\sum^{n}_{k=1} \frac{\phi(i)\phi(j^2)\phi(k^3)}{\phi(i)\phi(j)\phi(k)}\phi(gcd(i,j,k))\ \ \ 1\leq n\leq10^7$$
对上式作以下推导
到现在为止关于表达式的推导已经结束了,接下来就是怎么快速的计算出答案。
显然对于表达式的前半截$\sum^{n}_{t=1}\frac{\lfloor\frac{n}{t}\rfloor^3(\lfloor\frac{n}{t}\rfloor+1)^2(2\lfloor\frac{n}{t}\rfloor+1)}{12}$可以在每次询问的时候使用整数分块来回答。这儿还有一个不好处理的点,我们没办法求在模$2^{30}$的意义下12的逆元。于是我们观察上式$\lfloor\frac{n}{t}\rfloor和\lfloor\frac{n}{t}\rfloor+1$中必定有一个能够被2整除,于是我们在整数分块的时候先判断上式哪个能够被2整除并且除于两个2。然后再取模并对于最后的答案乘上3的逆元。计算3在模$2^{30}$的意义下的逆元不能用费马小定理,只能使用拓展欧几里得算法或者是欧拉定理。$\phi(2^{30})=\phi(1073741824)=536870912$最后求出来3的逆元为715827883。
对于后半截$设g(n)=d^3\sum_{d|n}\mu(\frac{n}{d})\phi(d)$,本来对于一般的题目这儿用$O(NlogN)$的方法预处理是完全足够的,但是$O(NlogN)$对于$10^7$还是太慢了啊。于是我们我们考虑到$\mu(x)和\phi(x)$为两个积性函数。
根据积性函数的性质,我们知道积性函数的狄利克雷卷积还是积性函数。所以$g(n)$为一个积性函数。我们可以使用线性筛来推得$g(n)$,要线性筛推$g(i*j)$对于i与j互质的情况$g(i*j)=g*(i)*g(j)$,我们可以把i与j互质的情况简化为需要推导出$g(p^k)$
这样的话就可以用线性筛递推求得g(n)
#include <bits/stdc++.h>
const int MAXN = 10000000;
const long long MOD = 1073741824;
using namespace std;
int T,n,cnt,inv;
int prime[MAXN+5],num[MAXN+5];
long long g[MAXN+5],ans;
bool vis[MAXN+5];
int pow_mod(long long a,int b){
long long ans = 1;
while(b){
if(b & 1) ans = ans * a % MOD;
a = a * a % MOD;
b >>= 1;
}
ans %= MOD;
return ans;
}
inline void init()
{
g[1] = 1;
for(int i = 2;i <= MAXN;i++){
if(!vis[i]){
prime[++cnt] = i;
g[i] = i - 2;
num[i] = 1;
}
for(int j = 1;j <= cnt && prime[j] * i <= MAXN;j++){
vis[prime[j]*i] = true;
if(i%prime[j] == 0){
num[i*prime[j]] = num[i] + 1;
if(num[i] == 1){
if(prime[j] == 2) g[i*prime[j]] = g[i/prime[j]];
else g[i*prime[j]] = g[i]/(prime[j]-2)*(prime[j]-1)*(prime[j]-1);
}
else g[i*prime[j]] = g[i]*prime[j];
break;
}
else {
g[prime[j]*i] = g[prime[j]] * g[i];
num[prime[j]*i] = 1;
}
}
}
for(int i = 1;i <= MAXN;i++) g[i] = (g[i-1]+g[i]*i%MOD*i%MOD*i%MOD)%MOD;
}
int main()
{
init();
inv = 715827883;
scanf("%d",&T);
while(T--){
scanf("%d",&n);
ans = 0;
int R;
long long m = 0,ret = 0;
for(int i = 1;i <= n;i = R + 1)
{
R = n/(n/i);
m = n/i;ret = (g[R]-g[i-1]+MOD)%MOD;
if(m%2 == 0)
{
ans += ((m/2)*(m/2)%MOD*m%MOD*(m+1)%MOD*(m+1)%MOD*(m*2+1)*inv%MOD*ret%MOD)%MOD;
}
else {
ans += (m*m%MOD*m%MOD*((m+1)/2)%MOD*((m+1)/2)%MOD*(m*2+1)*inv%MOD*ret%MOD)%MOD;
}
ans %= MOD;
//printf("i=%d R=%d m=%lld ret=%lld ans=%lld\n",i,R,m,ret,ans);
}
printf("%lld\n",ans);
}
return 0;
}
写的真好,懂了,谢谢学长