洛谷传送门
SPOJ传送门
BZOJ传送门
注:
B
Z
O
J
BZOJ
BZOJ数据范围稍有不同, 博主程序以
B
Z
O
J
BZOJ
BZOJ为标准实际上是在SPOJ上面卡不过QAQ
Description
Input
第一行一个自然数 T T T 表示数据组数
下接 T T T 行,每行一个正整数 n n n
T ≤ 10 , n ≤ 1 0 9 T ≤ 10, n ≤ 10^9 T≤10,n≤109
Output
输出 T T T 行,每行一个数字表示答案
Sample Input
4
233333
2333333
23333333
233333333
Sample Output
14145459
191698371
2494643347
31475021381
解题分析
上次在这道题中我们分析出了一个
σ
0
(
i
)
\sigma_0(i)
σ0(i)的性质, 这里还需要用到一个:
σ
0
(
n
2
)
=
∑
i
∣
n
2
ω
(
i
)
\sigma_0(n^2)=\sum_{i|n}2^{\omega(i)}
σ0(n2)=i∣n∑2ω(i)
其中
ω
(
n
)
\omega(n)
ω(n)表示
n
n
n的不同质因数的个数。
仍然考虑将 n n n分解为 n = p 1 a 1 ∗ p 2 a 2 ∗ p 3 a 3 ∗ . . . . . . ∗ p k a k n=p_1^{a_1}*p_2^{a_2}*p_3^{a_3}*......*p_k^{a_k} n=p1a1∗p2a2∗p3a3∗......∗pkak, 那么对于每个原有的约数 m m m, 那么我们可以让任何一个质数次幂 + a i +a_i +ai得到一个新的组合, 而这样枚举恰好是覆盖了所有情况的。
注意到 2 ω ( n ) 2^{\omega(n)} 2ω(n)还有一个意义: n n n的不含平方因子的约数个数。
因此实际上 2 w ( n ) = ∑ d ∣ n μ 2 ( d ) 2^{w(n)}=\sum_{d|n}\mu^2(d) 2w(n)=∑d∣nμ2(d)。
所以我们可以设
f
(
x
)
=
σ
0
(
x
2
)
f(x)=\sigma_0(x^2)
f(x)=σ0(x2),
g
(
x
)
=
2
ω
(
x
)
g(x)=2^{\omega(x)}
g(x)=2ω(x),
h
(
x
)
=
μ
2
(
x
)
h(x)=\mu^2(x)
h(x)=μ2(x), 那么有
f
(
x
)
=
g
(
x
)
∗
I
(
x
)
g
(
x
)
=
h
(
x
)
∗
I
(
x
)
f
(
x
)
=
h
(
x
)
∗
I
(
x
)
∗
I
(
x
)
=
h
(
x
)
∗
(
I
(
x
)
∗
I
(
x
)
)
=
h
(
x
)
∗
σ
0
(
x
)
f(x)=g(x)*I(x) \\ g(x)=h(x)*I(x) \\ f(x)=h(x)*I(x)*I(x) \\ =h(x)*(I(x)*I(x)) \\ =h(x)*\sigma_0(x)
f(x)=g(x)∗I(x)g(x)=h(x)∗I(x)f(x)=h(x)∗I(x)∗I(x)=h(x)∗(I(x)∗I(x))=h(x)∗σ0(x)
所以就有
∑
i
=
1
n
f
(
x
)
=
∑
i
=
1
n
μ
2
(
i
)
∑
j
=
1
⌊
n
i
⌋
σ
0
(
j
)
\sum_{i=1}^nf(x) \\ =\sum_{i=1}^n\mu^2(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\sigma_0(j)
i=1∑nf(x)=i=1∑nμ2(i)j=1∑⌊in⌋σ0(j)
而这两个函数也可以在
N
\sqrt N
N的时间内筛出来:
∑
i
=
1
n
μ
2
(
i
)
=
∑
i
=
1
n
μ
(
i
)
⌊
n
i
2
⌋
∑
i
=
1
n
σ
0
(
i
)
=
∑
i
=
1
n
⌊
n
i
⌋
\sum_{i=1}^n\mu^2(i)=\sum_{i=1}^{\sqrt n}\mu(i)\lfloor\frac{n}{i^2}\rfloor \\ \sum_{i=1}^n\sigma_0(i)=\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor
i=1∑nμ2(i)=i=1∑nμ(i)⌊i2n⌋i=1∑nσ0(i)=i=1∑n⌊in⌋
然后可以类似杜教筛, 处理出前 N 2 3 N^\frac{2}{3} N32项, 使得总复杂度在下底分块后达到 O ( N 2 3 ) O(N^\frac{2}{3}) O(N32)。
代码如下:
#include <cstdio>
#include <cctype>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <bitset>
#include <cstdlib>
#include <tr1/unordered_map>
#include <map>
#define R register
#define IN inline
#define W while
#define gc getchar()
#define MX 1000050
#define ll long long
template <class T>
IN void in(T &x)
{
x = 0; R char c = gc;
for (; !isdigit(c); c = gc);
for (; isdigit(c); c = gc)
x = (x << 1) + (x << 3) + c - 48;
}
int pcnt;
short miu[MX], cnt[MX];
int sgm0[MX], sum1[MX], pri[MX / 10], sum2[MX];
bool npr[MX];
std::tr1::unordered_map <ll, ll> mp1, mp2;
IN void get()
{
miu[1] = sgm0[1] = 1;
R int i, j, tar;
for (i = 2; i <= 1e6; ++i)
{
if (!npr[i]) miu[i] = -1, sgm0[i] = 2, pri[++pcnt] = i, cnt[i] = 1;
for (j = 1; j <= pcnt; ++j)
{
tar = i * pri[j];
if (tar > 1e6) break;
npr[tar] = true;
if (!(i % pri[j]))
{
miu[tar] = 0;
sgm0[tar] = sgm0[i] / (cnt[i] + 1) * (cnt[i] + 2);
cnt[tar] = cnt[i] + 1;
break;
}
miu[tar] = -miu[i];
sgm0[tar] = sgm0[i] << 1;
cnt[tar] = 1;
}
}
for (R int i = 1; i <= 1e6; ++i)
{
sum1[i] = sum1[i - 1] + (miu[i] != 0);
sum2[i] = sum2[i - 1] + sgm0[i];
}
}
IN ll solve1(R ll pos)//for sigma mu(i)^2
{
if (pos <= 1e6) return sum1[pos];
if (mp1.count(pos)) return mp1[pos];
ll bd = std::sqrt(pos);
ll ret = 0;
for (R ll i = 1; i <= bd; ++i) ret += miu[i] * (pos / (i * i));
return mp1[pos] = ret;
}
IN ll solve2(R ll pos)//for sigma sig0(i)
{
if (pos <= 1e6) return sum2[pos];
if (mp2.count(pos)) return mp2[pos];
ll ret = 0;
for (R ll lef = 1, rig; lef <= pos; lef = rig + 1)
{
rig = pos / (pos / lef);
ret += (pos / lef) * (rig - lef + 1);
}
return mp2[pos] = ret;
}
int main(void)
{
int T;
ll n, ans;
in(T); get();
W (T--)
{
in(n); ans = 0;
for (ll lef = 1, rig; lef <= n; lef = rig + 1)
{
rig = n / (n / lef);
ans += (solve1(rig) - solve1(lef - 1)) * solve2(n / lef);
}
printf("%lld\n", ans);
}
}