[bzoj2671] Calc

Description

给出N，统计满足下面条件的数对(a,b)的个数：

1. 1<=a<b<=N
2. a+b整除a*b　

Input

　一行一个数N

Output

　一行一个数表示答案

Sample Input

15

Sample Output

4

Solution

题目求的是：

[sum_{i=1}^nsum_{j=i+1}^n[i+j|ij] ]

不妨设(d=(i,j),xcdot d=i,ycdot d =j)，那么后面的式子就是：

[d(x+y)|xyd^2\ x+y|xyd ]

因为(x perp y)，可得：

[x+y|d ]

设(d=k(x+y))，得(j=yk(x+y)leqslant n)，所以对于互质的数对(x,y)，有(lfloorfrac{n}{y(x+y)} floor)的贡献。

考虑到(y(x+y)leqslant n)，有(yleqslantsqrt{n})，所以答案就是：

[ans=sum_{y=1}^sqrt{n}sum_{x=1}^{min(sqrt{n},y-1)}[(x,y)=1]lfloorfrac{n}{y(x+y)} floor ]

然后莫反下可得：

[ans=sum_{d=1}^sqrt{n}mu(d)sum_{y=1}^{lfloorfrac{sqrt{n}}{d} floor}sum_{x=1}^{y-1}lfloorfrac{n}{d^2y(x+y)} floor ]

然后后面数论分块下然后全程暴力，(凭借信仰过掉此题)

考虑下时间复杂度，总枚举次数大概是：

[sum_{i=1}^{sqrt{n}}frac{sqrt{n}}{i}sqrt{frac{sqrt{n}}{i}}=n^{frac{3}{4}}cdot sum_{i=1}^{sqrt{n}}i^{-frac{3}{2}} ]

(因为后面不会算)带入数据可得后面那个(sum)大概是2，所以总复杂度就是(O(n^frac{3}{4}))

update: 后面那个是黎曼函数(zeta(x)=sum_{i=1}^infty i^{-x})，可以证明当(x=1.5)时这个函数值收敛，且趋近于(2.6)左右，可以看作是(O(1))，具体证明可以看(wikipedia)，虽然我也不会。

#include<bits/stdc++.h>
using namespace std;

#define int long long

#define ONLINE_JUDGE

#ifdef ONLINE_JUDGE
#define getchar() ((p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin)),p1==p2)?EOF:*p1++)
#endif

namespace fast_IO {
	char buf[1<<21],*p1=buf,*p2=buf;

	template <typename T> inline void read(T &x) {
		x=0;T f=1;char ch=getchar();
		for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
		for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
	}
	template <typename T,typename... Args> inline void read(T& x,Args& ...args) {
		read(x),read(args...);
	}

	char buf2[1<<21],a[80];int p,p3=-1;

	inline void flush() {fwrite(buf2,1,p3+1,stdout),p3=-1;}
	template <typename T> inline void write(T x) {
		if(p3>(1<<20)) flush();
		if(x<0) buf2[++p3]='-',x=-x;
		do {a[++p]=x%10+48;} while(x/=10);
		do {buf2[++p3]=a[p];} while(--p);
		buf2[++p3]='
';
	}
	template <typename T,typename... Args> inline void write(T x,Args ...args) {
		write(x),write(args...);
	}
}

using fast_IO :: read;
using fast_IO :: write;
using fast_IO :: flush;

const int maxn = 1e5+10;

int vis[maxn],mu[maxn],pri[maxn],tot;

void sieve() {
	mu[1]=1;
	for(int i=2;i<maxn;i++) {
		if(!vis[i]) pri[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pri[j]<maxn;j++) {
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
			mu[i*pri[j]]=-mu[i];
		}
	}
}

int calc(int n,int m) {
	int ans=0;
	for(int i=1;i<=m;i++) {
		int t=n/i,T=i+1;
		while(T<=t&&T<i*2) {
			int pre=T;T=min(i*2-1,t/(t/T));
			ans+=(T-pre+1)*(t/T);T++;
		}
	}return ans;
}

signed main() {
	sieve();
	int n,ans=0,m;read(n);m=sqrt(n);
	for(int i=1;i<=m;i++) if(mu[i]) ans+=calc(n/(i*i),m/i)*mu[i];
	write(ans);
	flush();
	return 0;
}