Count a*b 数论推公式+求因子平方和

设$$f(x)=\sum _{a=0}^{x-1}\sum _{b=0}^{x-1}[x\nmid ab]$$
$$g(x)=\sum _{d\mid x}f(d)$$
求$g(n)\%2^{64}$，$T$组询问，$T\le 20000$, $n\le 2\ast 10^9$

解法：

$$f(x)=\sum _{a=1}^x\sum _{b=1}^x[x\nmid ab]$$
容斥得到，
$$=x^2-\sum _{a=1}^x\sum _{b=1}^x[x\mid ab]$$
a,x同时除去最大公因子得，
$$=x^2-\sum _{a=1}^x\sum _{b=1}^x[\frac{x}{gcd(a,x)}\mid \frac{ab}{gcd(a,x)}]$$
$\frac{a}{gcd(a,x)}$肯定和$\frac{x}{gcd(a,x)}$互质，所以，
$$=x^2-\sum _{a=1}^x\sum _{b=1}^x[\frac{x}{gcd(a,x)}\mid b]$$
${\sum }_{b=1}^x[\frac{x}{gcd(a,x)}\mid b]$其实求1到x有多少个数被$\frac{x}{gcd(a,x)}$整除，所以，
$$=x^2-\sum _{a=1}^x\frac{x}{\frac{x}{gcd(a,x)}}$$
$$=x^2-\sum _{a=1}^{x}gcd(a,x)$$
换枚举方式，枚举x的所有因子d
$$=x^2-\sum _{d\mid x}d\sum _{a=1}^x[gcd(a,x)==d]$$
$$=x^2-\sum _{d\mid x}d\sum _{a=1}^{\frac{x}{d}}[gcd(a,\frac{x}{d})==1]$$
引用欧拉函数的定义$\phi (x)$，
$$f(x)=x^2-\sum _{d\mid x}d\ast \phi (\frac{x}{d})$$
带入$g(n)$,
$$g(n)=\sum _{x\mid n}{x}^2-\sum _{x|n}\sum _{d|x}d\ast \phi (\frac{x}{d})$$
设$T=\frac{x}{d}$,
$$g(n)=\sum _{x\mid n}{x}^2-\sum _{dT|n}\sum _{d|n}d\ast \phi (T)$$
$$=\sum _{x\mid n}{x}^2-\sum _{T|\frac{n}{d}}\sum _{d|n}d\ast \phi (T)$$
$$=\sum _{x\mid n}{x}^2-\sum _{d\mid n}d\sum _{T|\frac{n}{d}}\phi (T)$$
注意到${\sum }_{T|\frac{n}{d}}\phi (T)=\frac{n}{d}$，
$$g(n)=\sum _{x\mid n}{x}^2-\sum _{d|n}d\ast \frac{{}_n}{{}^d}=\sum _{x\mid n}{x}^2-\sum _{d|n}n$$
$$g(n)=\sum _{d\mid n}{x}^2-n$$
这个式子是求n的所有因子的平方和减去所有因子的个数乘N的差，质因子分解后求就行。还有记得使用unsigned long long。

#include <iostream>
#include <cstdio>
#define ll long long
#define llu unsigned ll
using namespace std;

const int N = 1e5+10;

llu p[N], prime[N], m;
ll n;

void primes(int n){
	for(int i = 2; i <= n; i++){
		if(!p[i]) p[i] = i, prime[++m] = i;
		for(int j = 1; j <= m; j++){
			if(prime[j] > p[i] || prime[j] > n/i) break;
			p[i*prime[j]] = prime[j];
		}
	}
	//for(int i = 1; i <= m; i++) cout << prime[i] << endl;
}

int main(){
	primes(1e5);
	int t;
	scanf("%d", &t);
	while(t--){
		scanf("%lld", &n);
		llu ans1 = 1, ans2 = n;
		for(int i = 1; i <= m && prime[i] * prime[i]<= n; i++){
			if(n%prime[i] == 0){
				llu cnt = 1, tmp = 1, base = 1;
				while(n%prime[i] == 0){
					base *= prime[i];
					tmp += base*base;
					cnt++;
					n /= prime[i];
				}
				ans1 *= tmp;
				ans2 *= cnt;
				//123
			}
		}
		if(n > 1){
			ans1 *= n*n+1;
			ans2 *= 2;
		}
		printf("%llu\n", ans1-ans2);
	}
	return 0;
}