【codeforces 235E】 Number Challenge

http://codeforces.com/problemset/problem/235/E (题目链接)

题意

　　给出${a,b,c}$，求${sum_{i=1}^asum_{j=1}^bsum_{k=1}^cd(ijk)}$

extra

　　有这样一个公式，就是约数个数和那道题的推广吧。$${sum_{i=1}^asum_{j=1}^bsum_{k=1}^cd(ijk)=sum_{i=1}^asum_{j=1}^bsum_{k=1}^c[gcd(i,j)=gcd(i,k)=gcd(j,k)=1]lfloorfrac{a}{i} floorlfloorfrac{b}{j} floorlfloorfrac{c}{k} floor}$$

　　然并卵，这个式子我不会证也不会用。。还是直接推吧。

Solution

　　莫比乌斯反演。

egin{aligned}  sum_{i=1}^asum_{j=1}^bsum_{k=1}^cd(ijk)=sum_{i=1}^{ab}f(i)sum_{j=1}^cd(ij)   end{aligned}

　　其中${f(n)=sum_{i=1}^asum_{j=1}^b[ab=n]}$。

egin{aligned}     & sum_{i=1}^{ab}f(i)sum_{j=1}^cd(ij)  \  =&sum_{i=1}^{ab}f(i)sum_{j=1}^csum_{u|i}sum_{v|j}[gcd(u,v)=1]  \  =&sum_{u=1}^{ab}sum_{v=1}^c[gcd(u,v)=1]sum_{i=1}^{lfloor{ab/u} floor}f(iu)lfloorfrac{c}{v} floor     end{aligned}

　　我们令${S(n)=sum_{i=1}^{lfloor{ab/n} floor}f(in)}$，并把变量${u,v}$换成${i,j}$，因为${u,v}$看起来太丑了→_→。

egin{aligned}     & sum_{i=1}^{ab}sum_{j=1}^c[gcd(i,j)=1]S(i)lfloorfrac{c}{j} floor  \  =&sum_{i=1}^{ab}sum_{j=1}^csum_{t|i,t|j}μ(t)S(i)lfloorfrac{c}{j} floor  \  =&sum_{t=1}^cμ(t)sum_{i=1}^{lfloor{ab/t} floor}S(it)sum_{j=1}^{lfloor{c/t} floor}lfloorfrac{c}{jt} floor     end{aligned}

　　看到这个式子是不是感到了满满的套路，我们故技重施，令${P(n)=sum_{i=1}^{lfloor{ab/n} floor}S(in)}$，${Q(n)=sum_{i=1}^nlfloorfrac{n}{i} floor}$

egin{aligned}     sum_{t=1}^cμ(t)P(t)Q(lfloorfrac{c}{t} floor)     end{aligned}

　　这样，我们${O(ab)}$的求出${f}$，之后按顺序${O(ablog ab)}$的求出${S}$和${P}$，然后${O(c^2)}$的求出${Q}$，当然如果想分段也可以分段。最后只需要从${1}$枚举到${c}$就可以得出答案了。

细节

　　时限卡得有点紧，你需要常数优化→_→

代码

#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<queue>
#define LL long long
#define inf 2147483640
#define MOD (1ll<<30)
#define Pi acos(-1.0)
#define free(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout);
using namespace std;

const int maxn=4000010,maxm=2010;
LL s[maxn],f[maxn],P[maxm],Q[maxm];
int a,b,c,p[maxm],vis[maxm],mu[maxm];

void inc(LL &a,LL b) {
    a+=b;if (a>MOD) a-=MOD;
}
int main() {
	scanf("%d%d%d",&a,&b,&c);
	mu[1]=1;
	for (int i=2;i<=c;i++) {
		if (!vis[i]) p[++p[0]]=i,mu[i]=-1;
		for (int j=1;j<=p[0] && i*p[j]<=c;j++) {
			vis[i*p[j]]=1;
			if (i%p[j]==0) {mu[i*p[j]]=0;break;}
			else mu[i*p[j]]=-mu[i];
		}
	}
	for (int i=1;i<=a;i++)
		for (int j=1;j<=b;j++) f[i*j]++;
	for (int i=1;i<=a*b;i++)
		for (int j=1;j<=a*b/i;j++) inc(s[i],f[i*j]);
	for (int i=1;i<=c;i++)
		for (int j=1;j<=a*b/i;j++) inc(P[i],s[i*j]);
	for (int i=1;i<=c;i++)
		for (int j=1;j<=i;j++) inc(Q[i],i/j);
	LL ans=0;
	for (int i=1;i<=c;i++)
		ans=(ans+mu[i]*P[i]*Q[c/i]%MOD)%MOD;
	printf("%lld",(ans+MOD)%MOD);
	return 0;
}