题解:
必知知识:
σ(n)是n的约数和。
若
n=∏pqii
则
σ(n)=∏∑qij=0pji
=∏pqi+1i−1pi−1
知道了这个就可以线性筛法了,如有不懂得见code。
必要结论:
σ(i∗j)=∑p|i∑q|jp∗jq(gcd(p,q)=1)
证明:
∑p|i∑q|jp∗jq(gcd(p,q)=1)
=∑p|i∑q|jp∗q(gcd(p,jq)=1)
对每一个质因子分类讨论:
设
i=∑puii
,
j=pvii
p=∑plii
,
q=∑prii
1.若
ui>=li>0
,则
ri=vi
,此时可以表示出
p(vi+1)−>(ui+vi)i
。
2.
li=0
,则
vi>=ri>=0
,此时可以表示出
p0−>vii
.
于是可以得出 p∗q 一定可以表示出每个质因子 p0−>(ui+vi)i ,合起来就包括了所有i*j的所有约数。
∴ σ(i∗j)=∑p|i∑q|jp∗jq(gcd(p,q)=1)
有了gcd就可以反演了。
Ans=∑ni=1∑nj=1σ(i∗j)∗max(i,j)
=∑ni=1∑nj=1∑p|i∑q|jp∗jq∗max(i,j)(gcd(p,q)=1)
=∑np=1∑nq=1∑⌊np⌋i=1∑⌊nq⌋j=1p∗j∗max(pi,qj)(gcd(p,q)=1)
=∑np=1∑nq=1∑⌊np⌋i=1∑⌊nq⌋j=1p∗j∗max(pi,qj)∑d|gcd(i,j)μ(d)
=∑nd=1μ(d)∗d∑⌊nd⌋p=1∑⌊nd⌋q=1∑⌊nd∗p⌋i=1∑⌊nd∗q⌋j=1p∗j∗max(d∗pi,d∗qj)
=∑nd=1μ(d)∗d2∗∑⌊nd⌋p=1∑⌊nd⌋q=1∑⌊nd∗p⌋i=1∑⌊nd∗q⌋j=1p∗j∗max(pi,qj)
注意到d后面是和n/d有关的一个函数。
设
g(n/d)=∑⌊nd⌋p=1∑⌊nd⌋q=1∑⌊nd∗p⌋i=1∑⌊nd∗q⌋j=1p∗j∗max(pi,qj)
g(n)=∑np=1∑nq=1∑⌊np⌋i=1∑⌊nq⌋j=1p∗j∗max(pi,qj)
仔细观察你发现这就是一个σ.
g(n)=∑ni=1∑nj=1σ(i)∗σ(j)∗max(i,j)
=(∑ni=1σ(i)∗i∗∑ij=1∗σ(j))∗2−(∑ni=1i∗σ(i)2)
用线性筛法预处理1-n的σ,就可以预处理g(1-n)了。
转回去,现在给出n,求
∑nd=1g(⌊nd⌋)∗μ(d)∗d2
设
M(d)=μ(d)∗d2
则求
∑nd=1g(⌊nd⌋)∗M(d)
1≤数据组数≤50000,1≤N≤1000000。
不要以为你的暴力分块能过。
用差分表O(n log n)预处理即可O(1)查询了。
Code:
#include<cstdio>
#include<cstring>
#define ll long long
#define fo(i, x, y) for(int i = x; i <= y; i ++)
#define min(a, b) ((a) < (b) ? (a) : (b))
#define max(a, b) ((a) > (b) ? (a) : (b))
using namespace std;
const ll N = 1000000;
const ll mo = 1e9 + 7, ni_2 = 5e8 + 4, ni_6 = 166666668;
bool bz[N + 5];
ll p[N];
ll u[N + 5], v[N + 5], s[N + 5];
ll mu[N + 5], s2[N + 5], s3[N + 5], s4[N + 5];
int T, n; ll ans;
int main() {
mu[1] = u[1] = v[1] = s[1] = 1;
fo(i, 2, N) {
if(!bz[i]) p[++ p[0]] = i, mu[i] = -1, u[i] = v[i] = i, s[i] = i + 1;
fo(j, 1, p[0]) {
ll k = i * p[j];
if(k > N) break;
bz[k] = 1;
if(i % p[j] == 0) {
mu[k] = 0;
u[k] = u[i]; v[k] = v[i] * u[k];
s[k] = s[k / v[k]] * ((v[k] * u[k] - 1) / (u[k] - 1));
break;
}
u[k] = v[k] = p[j];
s[k] = s[i] * s[p[j]];
mu[k] = -mu[i];
}
}
fo(i, 1, N) mu[i] = mu[i] * i % mo * i % mo;
fo(i, 1, N) s2[i] = (s2[i - 1] + s[i]) % mo;
fo(i, 1, N) s3[i] = (s3[i - 1] + s[i] * i % mo * s2[i] % mo) % mo;
fo(i, 1, N) s4[i] = (s4[i - 1] + s[i] * s[i] % mo * i % mo) % mo;
fo(i, 1, N) s[i] = (s3[i] * 2 - s4[i] + mo) % mo;
memset(s2, 0, sizeof(s2));
fo(i, 1, N) fo(j, 1, N / i) {
ll x = i * j, y = min(N, i * (j + 1) - 1);
s2[x] += mu[i] * s[j] % mo; s2[y + 1] -= mu[i] * s[j] % mo;
}
fo(i, 1, N) s3[i] = ((s3[i - 1] + s2[i]) % mo + mo) % mo;
scanf("%lld", &T);
fo(ii, 1, T) {
ans = 0;
scanf("%d", &n);
printf("Case #%d: %d\n", ii, s3[n]);
}
}