loj2085题解

挺优美的一道数论题. 当然也有可能是我太菜了.

首先显然是考虑只计算既约分数来避免算重,即只计算$x\perp y$. 考虑$k=10$的时候,$\frac{x}{y}$符合题目条件当且仅当$y$没有因数$2,5$. 容易证明对于所有情况,$\frac{x}{y}$符合题目条件等价于$y\perp k$. 所以我们要计算的其实就是

$$\sum _{x=1}^n\sum _{y=1}^m[x\perp y][y\perp k]$$

一开始我很傻地把两个中括号都拆开了,其实应该先只拆一个观望一下. 因为$[x\perp y]$这一项复杂一些,所以考虑先拆这个(其实还是在瞎猜),得到

$$\begin{aligned}&\sum _{y=1}^m[y\perp k]\sum _{x=1}^n\sum _{d|x,~d|y}\mu(d)\\=&\sum _{y=1}^m[y\perp k]\sum _{d|y}^n\mu(d)\left\lfloor\frac{n}{d}\right\rfloor\\=&\sum _{d=1}^n\mu(d)\left\lfloor\frac{n}{d}\right\rfloor\sum _{y=1}^{\lfloor m/d\rfloor}[dy\perp k]\\=&\sum _{d=1}^n[d\perp k]\mu(d)\left\lfloor\frac{n}{d}\right\rfloor\sum _{y=1}^{\lfloor m/d\rfloor}[y\perp k]\end{aligned}$$

定义$f(n)=\sum _{y=1}^m[y\perp k]$,那么有$f(km+r)=mf(k)+f(r)~(0\le r\lt k)$,所以预处理$f(1),f(2),\ldots,f(k)$就可以$\mathcal O(1)$计算了. 之前的式子化为

$$\sum _{d=1}^n[d\perp k]\mu(d)\left\lfloor\frac{n}{d}\right\rfloor f(\left\lfloor\frac{m}{d}\right\rfloor)$$

对这个式子整除分块,于是我们需要快速计算

$$f(n,k)=\sum _{d=1}^n[d\perp k]\mu(d)$$

如果没有$[d\perp k]$,那就是一个简单的杜教筛,现在加上这一项,它还是积性函数没关系就可以考虑从总和里面去掉不互质的,设$k=p^rq$,其中$p$为素数,$p\perp q$,那么有

$$\begin{aligned}f(n,k)=&f(n,q)-\sum _{y=1}^{\lfloor n/p\rfloor}[py\perp q]\mu(py)\\=&f(n,q)-\sum _{y=1}^{\lfloor n/p\rfloor}[y\perp q]\mu(py)\\=&f(n,q)-\sum _{y=1}^{\lfloor n/p\rfloor}[y\perp q][y\perp p]\mu(p)\mu(y)\\=&f(n,q)+\sum _{y=1}^{\lfloor n/p\rfloor}[y\perp k]\mu(y)\\=&f(n,q)+f(\left\lfloor\frac{n}{p}\right\rfloor,k)\end{aligned}$$

然后注意到$f(n,1)$可以直接杜教筛,这道题就可以做了.

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define ele long long
using namespace std;
#define maxn 1000010
#define maxk 2010
template<class T1,class T2>struct hashmap{
static const ele R=8000009;
T1 a[R];
T2 b[R];
bool c[R];
hashmap(){ memset(c,0,sizeof(c)); }
inline ele find(T1 x){
ele i=x%R;
for (; c[i] && a[i]!=x; (i+=1)%=R);
return i;
}
inline bool test(T1 x){ return c[find(x)]; }
inline T2& operator[](T1 x){
ele i=find(x);
!c[i] && (c[i]=true,a[i]=x,b[i]=0);
return b[i];
}
};
ele n,m,k,S,f[maxk],pcnt,plst[maxn],mu[maxn],qcnt,q[maxk];
bool flag[maxn];
hashmap<ele,ele> g;
ele gcd(ele a,ele b){
return b?gcd(b,a%b):a;
}
ele calc(ele n){
if (n<=S) return mu[n];
if (g.test(n)) return g[n];
g[n]=1;
for (ele d=n; d>1; ){
ele d1=max(1ll,n/(n/d+1));
g[n]-=calc(n/d)*(d-d1);
d=d1;
}
return g[n];
}
inline ele calcf(ele n){
return f[k]*(n/k)+f[n%k];
}
ele calcg(ele n,ele i){
if (!i) return calc(n);
if (!n) return 0;
if (n==1) return 1;
return calcg(n,i-1)+calcg(n/q[i-1],i);
}
int main(){
scanf("%lld%lld%lld",&n,&m,&k);
S=exp(log(n)/3*2);
f[0]=0;
for (int i=1; i<=k; ++i) f[i]=f[i-1]+(gcd(i,k)==1);
mu[1]=1; pcnt=0; memset(flag,0,sizeof(flag));
for (int i=2; i<=S; ++i){
if (!flag[i]) plst[pcnt++]=i,mu[i]=-1;
for (int j=0; j<pcnt && i*plst[j]<=S; ++j){
flag[i*plst[j]]=true;
if (i%plst[j]) mu[i*plst[j]]=-mu[i];
else{
mu[i*plst[j]]=0;
break;
}
}
}
mu[0]=0;
for (int i=2; i<=S; ++i) mu[i]+=mu[i-1];
ele tmp=k; qcnt=0;
for (int i=2; i*i<=k; ++i)
if (tmp%i==0){
q[qcnt++]=i;
while (tmp%i==0) tmp/=i;
}
tmp>1 && (q[qcnt++]=tmp);
ele ans=0;
for (ele d=min(n,m); d; ){
ele d1=n/(n/d+1),d2=m/(m/d+1);
d1=max(d1,d2);
ans+=calcf(m/d)*(n/d)*(calcg(d,qcnt)-calcg(d1,qcnt));
d=d1;
}
printf("%lld\n",ans);
return 0;
}