2018 ICPC 徐州 计蒜客 - Easy Math
题目给定
m<2∗109,n<1012
m
<
2
∗
10
9
,
n
<
10
12
计算:
显然: μ(n)=0 μ ( n ) = 0 时 ans=0 a n s = 0
当 mu(n)!=0 m u ( n ) ! = 0 时:
记:
显然,这里 F(k) F ( k ) 代表下指标遍历所有 gcd g c d 是 k k 的倍数的情况,那么必须是 n n 的约数才有意义:
反演有:
则:
记:
显然:
A(m,n)=μ(n)∑d|nμ(d)A(⌊md⌋,d)A(0,k)=0A(k,1)=M(k)A(1,k)=μ(k)
A
(
m
,
n
)
=
μ
(
n
)
∑
d
|
n
μ
(
d
)
A
(
⌊
m
d
⌋
,
d
)
A
(
0
,
k
)
=
0
A
(
k
,
1
)
=
M
(
k
)
A
(
1
,
k
)
=
μ
(
k
)
其中,
M(n)=∑ni=1μ(i)
M
(
n
)
=
∑
i
=
1
n
μ
(
i
)
计算方法参见:https://blog.csdn.net/ZLH_HHHH/article/details/77860249
在计算
M
M
的同事保留计算A所需的计算信息
计算A的时间可以表示为:
总时间复杂度不超过:
#include <algorithm>
#include <stdio.h>
#include <string.h>
#include <cmath>
#define MAXN 2000005
using namespace std;
typedef long long LL;
LL M[MAXN] = { 0,-1 };
LL tmp[MAXN];
LL mu[MAXN + 2];
void clat(LL n, int d)
{
LL ans = 0;
int m = (int)sqrt(n) + 1;
for (int L = 1;L<m;L++)
ans += M[L] * (n / L - n / (L + 1));
for (int i = (int)(n / m);i>1;i--)
{
LL u = n / i;
if (u<MAXN)
ans += M[u];
else
ans += tmp[i*d];
}
tmp[d] = 1 - ans;
}
LL slove(LL n)
{
if (n<MAXN)return M[n];
for (int i = (int)(n / (MAXN - 1));i;i--) clat(n / (LL)i, i);
return tmp[1];
}
LL DFS(LL m, LL c, LL n, LL mn)
{
if (n == 1)
{
if (m < MAXN)
return M[m];
else
return tmp[c];
}
if (m == 1)return mn;
if (m == 0)return 0;
LL ans = 0;
for (LL d = 1;d*d < n;d++)
{
if (n%d)continue;
ans += mu[d] * DFS(m / d, c*d, d, mu[d]);
int u_ = mn / mu[d];
LL d_ = n / (LL)d;
ans += u_*DFS(m / d_, c*d_, d_, u_);
}
return ans*mn;
}
int main()
{
for (int i = 1;i<MAXN;i++)
{
M[i] = -M[i];
mu[i] = M[i];
for (int j = i << 1;j<MAXN;j += i)M[j] += M[i];
M[i] += M[i - 1];
}
LL n, m;
scanf("%lld %lld", &m, &n);
LL mn = 1;
slove(m);
if (n >= MAXN)
{
LL a = n;
for (LL i = 2;i*i <= n;i++)
{
if (a%i == 0)
{
a /= i;
mn = -mn;
if (a%i == 0)
{
mn = 0;
break;
}
}
}
if (mn != 0 && a > 1)mn = -mn;
}
else
mn = mu[n];
if (mn == 0)
{
printf("0\n");
return 0;
}
printf("%lld\n", DFS(m, 1, n, mn));
return 0;
}