UVA 11426 GCD-Extreme(II) ★ (欧拉函数)

题意

求Σ_{1<=i<N} Σ_{i<j<=N} GCD(i, j)     (N<=4000000)

分析

原始思路

暴力求明显是不行的，我们把式子简化形式一下发现它可以写成Σ_{2<=j<=N} GCD(1~j-1, j) 这个形式就给我们一种思路：可以只枚举j，然后快速算出GCD(1~j-1, j) 我们当然不能枚举1~j-1那么算，那么再换种思路，枚举可能的答案k，即j的所有约数。分别计算GCD(1~j-1, j) = k的方案数（HDU 1695），然后加起来就能求出和了。 【求GCD(x, j) = k，x ∈ (1, j-1)的个数】我们知道j必须是k的倍数，所以可以在等式两边同时除以k变成：GCD(x, j/k) = 1，x ∈ (1, j/k-1)。那么显然答案等于phi(j/k)。

进一步优化

上面的方法还是超时，我们还需要优化。在求j的约数时会有很多无用状态，它的过程很像是暴力判断素数一样，联想到筛法求素数，我们也可以想到类似的思路：枚举k，那么k的倍数就是j，然后再算GCD(1~j-1, j) = k。  

#include 
#include 
#include 
#include 
#include 
#include 
#define MID(x,y) ((x+y)/2)
#define mem(a,b) memset(a,b,sizeof(a))
using namespace std;
const int MAX = 4000002;
int phi[MAX];
bool noprime[MAX];
vector  prime;
void Euler(int n){
    phi[1] = 0;
    for (int i = 2; i <= n; i ++){
        if (!noprime[i]){
            prime.push_back(i);
            phi[i] = i - 1;
        }
        for (int j = 0; j < (int)prime.size() && prime[j] * i <= n; j ++){
            noprime[prime[j]*i] = 1;
            if (i % prime[j] == 0){
                phi[prime[j]*i] = phi[i] * prime[j];
            }
            else{
                phi[prime[j]*i] = phi[i] * phi[prime[j]];
            }
        }
    }
}
long long f[MAX];
int main(){
	//freopen("test.in", "r", stdin);
	//freopen("test.out", "w", stdout);
	Euler(MAX);
	mem(f, 0);
	for (int i = 1; i < MAX; i ++){
        for (int n = i*2; n < MAX; n += i){
            f[n] += i * phi[n/i];
        }
    }
	int n;
	while(scanf("%d", &n), n){
        long long res = 0;
        for (int i = 2; i <= n; i ++){
            res += f[i];
        }
        printf("%lld
", res);
	}
	return 0;
}