HDU 5528 Count a * b 欧拉函数

题意：

定义函数(f(n))为(i cdot j otequiv 0 ; (mod ; n))的数对((i,j))的个数((0 leq i,j leq n))
(g(n)=sum_{d|n}f(d))，求(g(n) ; mod ; 2^{64})，其中(1 leq n leq 10^9)

仰慕叉老师手推公式

分析：

计算(f(n))

(f(n))不容易计算，但可以计算它的反面，也就是计算满足(ij equiv 0 ; (mod ; n))的数对((i,j))的个数。
为了方便计算，对(i,j)的范围稍作改动，(1 leq i,j leq n)，这不会影响最终计算结果。
(i cdot j)是(n)的整数倍，也就是(n)被分到了这里面。
不妨设(i=ax,j=by)，其中(ab=n)。我们可以枚举(a)来统计答案，但是这样有重复计算。
比如((3 cdot 2) cdot (2 cdot 1) = (6 cdot 1) cdot (2 cdot 1) = 6 cdot 2)，这样数对((6, 2))在(a=3)和(a=6)时被计算了两次。
所以计数时我们规定，数对((i,j))只会在(a)最大的时候被计算，所以有(gcd(b,x)=1)成立。
答案为(sumlimits_{a|n} sumlimits_{1 leq x leq frac{n}{a}}[gcd(x,b)=1] cdot frac{n}{b}=sumlimits_{a|n}a phi(frac{n}{a}))

计算(g(n))

根据题中给的关系，得到(g(n)=sumlimits_{m|n}(m^2-f(m))=sumlimits_{m|n}(m^2-sumlimits_{d|m}d phi(frac{m}{d})))
继续化简后面那部分，(sumlimits_{m|n}sumlimits_{d|m}d phi(frac{m}{d}))交换求和顺序，得到(sumlimits_{d|n}dsumlimits_{frac{m}{d}|frac{n}{d}}phi(frac{m}{d}))
由于(sumlimits_{d|n} phi(d)=n)，所以得到(sumlimits_{d|n}d cdot frac{n}{d}=nsumlimits_{d|n}1=n sigma(n))。其中(sigma(n))为(n)的约数的个数。

(sigma(n))可以这样计算：
首先将(n)质因数分解为(prodlimits_{i=1}^{r}p_{i}^{e_i})，这样(sigma(n)=prodlimits_{i=1}^{r}(e_i+1))。
因为我们已经将(n)分解，所以计算前面那部分的时候不要在(sqrt{n})内枚举(n)的约数，会超时的。要根据分解出来的结果直接枚举。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;

const int maxn = 32000 + 10;
int pcnt, prime[maxn];
bool vis[maxn];

void preprocess() {
	for(int i = 2; i < maxn; i++) {
		if(!vis[i]) prime[pcnt++] = i;
		for(int j = 0; j < pcnt && i * prime[j] < maxn; j++) {
			vis[i * prime[j]] = true;
			if(i % prime[j] == 0) break;
		}
	}
}

int tot, p[maxn], e[maxn];

unsigned long long ans;
void dfs(int d, unsigned long long product) {
	if(d == tot) {
		ans += product * product;
		return ;
	}

	dfs(d + 1, product);
	for(int i = 1; i <= e[d]; i++) {
		product *= p[d];
		dfs(d + 1, product);
	}
}

int main()
{
	preprocess();
	int T; scanf("%d", &T);
	while(T--) {
		int n;
		scanf("%d", &n);

		int t = n;
		tot = 0;
		for(int i = 0; i < pcnt && prime[i] * prime[i] <= n; i++) {
			if(t % prime[i] == 0) {
				p[tot] = prime[i];
				e[tot] = 0;
				while(t % prime[i] == 0) {
					t /= prime[i];
					e[tot]++;
				}
				tot++;
			}
		}
		if(t > 1) { p[tot] = t; e[tot++] = 1; }

		unsigned long long tmp = 1;
		for(int i = 0; i < tot; i++) tmp = tmp * (e[i] + 1);
		tmp *= n;

		ans = 0;
		dfs(0, 1);

		ans -= tmp;

		printf("%llu
", ans);
	}

	return 0;
}