Description:
1<=n<=10^10,k = 1、40
题解:
一看就是反演题,立马进行简单变形:
原式=
∑kx=1∑nd=1fx(d)∗∑ni=1∑nj=1[gcd(i,j)=d]
=∑kx=1∑nd=1fx(d)∗∑⌊ni⌋i=1∑⌊ni⌋j=1[gcd(i,j)=1]
=∑kx=1∑nd=1fx(d)∗((∑⌊ni⌋i=1φ(i)∗2)−1)
当k=1时,
f1(x)=μ(x)
,直接杜教筛μ和φ的前缀和即可。
当k=40时,就不好做了,首先需要注意到2^40>10^10。
设
λ(n)=fmax(n)
若
n=∏pqii
,则
λ(n)=(−1)∑qi
。
设
Fd(n)=∑ni=1fd(i)
。
实际上可以用λ来容斥出F_d,思路是这样的:
用总和,减去指数至少为d+1的,就可以枚举一个i,那么
id+1∗j(j<=n/id+1)
一定是指数大于i+1的,这样求会有重复,直接乘上一个μ即可。
所以:
Fd(n)=∑ni=1μ(i)∗∑nid+1j=1λ(id+1∗j)
λ显然是一个完全积性函数,所以:
Fd(n)=∑ni=1μ(i)∗λ(id+1)∗∑nid+1j=1λ(j)
缩小范围:
Fd(n)=∑n1d+1i=1μ(i)∗λ(id+1)∗∑nid+1j=1λ(j)
注意式子后面是λ的前缀和,杜教筛求出这个,对每个 Fd(n) 就能快速求了,复杂度证明不会,但是感性认识一下,这个是指数级减小的,所以应该是很小很小的。
最后一个问题是如何杜教筛λ。
∑d|nλ(d)
=∑d|nλ(d)∗I(nd)
这是个裸的狄利克雷卷积式,假设
n=pq
,随便算一下,易得:
∑d|nλ(d)=[n是完全平方数]
∑ni=1∑j|iλ(j)
=n√
=∑ni=1∑⌊nj⌋j=1λ(j)
设
S(n)=∑ni=1λ(i)
则
S(n)=n√−∑ni=2S(⌊ni⌋)
#include<cmath>
#include<cstdio>
#define ll long long
#define ld long double
#define fo(i, x, y) for(ll i = x; i <= y; i ++)
#define max(a, b) ((a) > (b) ? (a) : (b))
using namespace std;
const ll N = 5000000, M = 16516513;
const ll mo = 1 << 30, ni_2 = 5e8 + 4;
ll n, m;
bool bz[N + 5];
ll p[N], s[N + 5], mu[N + 5], v[N + 5], c[N + 5], phi[N + 5];
ll la[N + 5], cc[N + 5], sm[N + 5], ss[N + 5];
ll ans;
ll h[M][2];
ll hash(ll x) {
ll y = x % M;
while(h[y][0] != 0 && h[y][0] != x) y = (y + 1) % M;
return y;
}
ll dg(ll n) {
if(n <= N) return s[n];
ll y = hash(n);
if(h[y][0] == n) return h[y][1];
h[y][0] = n;
ll ans = sqrt(n * 1.0);
fo(i, 2, n) {
ll j = n / (n / i);
ans -= (j - i + 1) * dg(n / i) % mo;
i = j;
}
ans = (ans % mo + mo) % mo;
h[y][1] = ans;
return ans;
}
ll h2[M][2];
ll hash2(ll x) {
ll y = x % M;
while(h2[y][0] != 0 && h2[y][0] != x) y = (y + 1) % M;
return y;
}
ll dg2(ll n) {
if(n <= N) return phi[n];
ll y = hash2(n);
if(h2[y][0] == n) return h2[y][1];
h2[y][0] = n;
ll n_ = n % mo;
ll ans = n_ * (n_ + 1) / 2 % mo;
fo(i, 2, n) {
ll j = n / (n / i);
ans -= (j - i + 1) * dg2(n / i); ans %= mo;
i = j;
}
ans = (ans % mo + mo) % mo;
h2[y][1] = ans;
return ans;
}
void Build() {
mu[1] = 1; c[1] = 1; phi[1] = 1;
fo(i, 2, N) {
if(!bz[i]) p[++ p[0]] = i, v[i] = 1, c[i] = 1, s[i] = 1, mu[i] = -1, phi[i] = i - 1;
fo(j, 1, p[0]) {
int k = i * p[j];
if(k > N) break;
bz[k] = 1;
s[k] = s[i] ^ 1;
if(i % p[j] == 0) {
v[k] = v[i] + 1;
c[k] = max(c[i], v[k]);
mu[k] = 0;
phi[k] = phi[i] * p[j];
break;
}
v[k] = 1;
c[k] = c[i];
mu[k] = -mu[i];
phi[k] = phi[i] * phi[p[j]];
}
}
}
ll SS(ll p) {
if(p <= N) return ss[p];
ll gg = sqrt(p * 1.0);
fo(k, 1, gg) cc[k] = k;
ll sum = 0;
fo(d, 1, m) {
ll mm = pow(p, (ld) 1 / (d + 1));
fo(i, 1, mm) {
cc[i] *= i;
ll lb = (((d + 1) % 2) ? la[i] : 1) * mu[i];
sum += lb * dg(p / cc[i]); sum %= mo;
}
}
return sum;
}
int main() {
Build();
scanf("%lld %lld", &n, &m);
fo(i, 1, N) phi[i] = (phi[i - 1] + phi[i]) % mo;
fo(i, 1, N) {
ss[i] = ss[i - 1];
if(m == 1)
ss[i] = (ss[i] + (c[i] <= 1) * (s[i] ? -1 : 1)) % mo; else
ss[i] = (ss[i] + (s[i] ? -1 : 1) * (m - c[i] + 1)) % mo;
}
fo(i, 1, N) la[i] = s[i] ? -1 : 1, s[i] = s[i - 1] + la[i];
fo(i, 1, N) sm[i] = sm[i - 1] + mu[i];
fo(i, 1, n) {
ll c = n / i;
ll j = n / c;
ans += (SS(j) - SS(i - 1)) * (2 * dg2(n / i) - 1); ans %= mo;
i = j;
}
ans = (ans % mo + mo) % mo;
printf("%lld", ans);
}