数论函数练习记录（二）

参考资料：oi-wiki 杜教筛｜blog｜blog

[luoguP3768]简单的数学题

题目链接

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)&=\sum_{d=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}ijd^3[\gcd(i,j)=1]&&枚举\;d=\gcd(i,j)\\ &=\sum_{d=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}ijd^3\sum_{r\mid \gcd(i,j)}\mu(r)&&\epsilon=\mu*1\\ &=\sum_{d=1}^nd^3\sum_{r=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{i=1}^{\left\lfloor\frac{n}{rd}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{rd}\right\rfloor}ijr^2\mu(r)&&求和号换序\\ &=\sum_{d=1}^nd^3\sum_{r=1}^{\left\lfloor\frac{n}{d}\right\rfloor}r^2\mu(r)\left(\frac{\left\lfloor\frac{n}{rd}\right\rfloor\left(\left\lfloor\frac{n}{rd}\right\rfloor+1\right)}{2}\right)^2\\ &=\sum_{T=1}^n\sum_{d\mid T}d^3\frac{T^2}{d^2}\mu\left(\frac{T}{d}\right)\left(\frac{\left\lfloor\frac{n}{T}\right\rfloor\left(\left\lfloor\frac{n}{T}\right\rfloor+1\right)}{2}\right)^2&&令\;T=rd\\ &=\sum_{T=1}^n\left(\frac{\left\lfloor\frac{n}{T}\right\rfloor\left(\left\lfloor\frac{n}{T}\right\rfloor+1\right)}{2}\right)^2T^2\sum_{d\mid T}d\mu\left(\frac{T}{d}\right)\\ &=\sum_{T=1}^n\left(\frac{\left\lfloor\frac{n}{T}\right\rfloor\left(\left\lfloor\frac{n}{T}\right\rfloor+1\right)}{2}\right)^2T^2\varphi(T)&&\mu*\text{id}=\varphi \end{align}$$

于是问题转化为求出 $f(n)=n^2\varphi(n)$ 的前缀和，如此便可以数论分块了。

设 $g(n)=n^2$，那么 $(f*g)(n)=\sum\limits_{d\mid n}d^2\varphi(d)\frac{n^2}{d^2}=n^2\sum\limits_{d\mid n}\varphi(d)=n^3$，根据杜教筛理论，有：

$$S_f(n)=\sum_{i=1}^ni^3-\sum_{i=2}^ni^2S_f\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$$

即可计算。

#include<bits/stdc++.h>

using namespace std;

template<typename T>void read(T&x){x=0;int fl=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')
fl=-1;ch=getchar();}while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}x*=fl;}
template<typename T,typename...Args>inline void read(T&t,Args&...args){read(t);read(args...);}

typedef long long LL;
typedef vector<int> vi;
typedef pair<int, int> pii;
#define mp(x, y) make_pair(x, y)
#define pb(x) emplace_back(x)

const int N = 5000005;
LL p, inv2, inv6;

int phi[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	phi[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			phi[i] = (i - 1) % p;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				phi[i * pList[j]] = 1ll * phi[i] * pList[j] % p;
				break;
			}
			else	phi[i * pList[j]] = 1ll * phi[i] * (pList[j] - 1) % p;
		}
	}
}

inline LL fpow(LL bs, LL idx, LL mod){
	LL res = 1;
	bs %= mod;
	while(idx){
		if(idx & 1)	(res *= bs) %= mod;
		(bs *= bs) %= mod;
		idx >>= 1;
	}
	return res;
}

unordered_map<LL, LL> s;
LL pres[N];
LL apiadu(LL n){
	if(n <= 5000000)	return pres[n];
	if(s.find(n) != s.end())	return s[n];
	LL res = fpow(n%p*(n+1)%p*inv2%p, 2, p);
	for(LL l = 2, r; l <= n; l = r + 1){
		r = n / (n / l);
		LL rs = r%p*(r+1)%p*(2*r%p+1)%p*inv6%p;
		LL ls = (l-1)%p*l%p*(2*l%p-1)%p*inv6%p;
		res -= (rs - ls) * apiadu(n / l) % p;
		((res %= p) += p) %= p;
	}
	return s[n] = res;
}

LL n;
inline LL solve(){
	LL res = 0;
	for(LL l = 1, r; l <= n; l = r + 1){
		r = n / (n / l);
		res += fpow((n/l)%p*(n/l+1)%p, 2, p) * (apiadu(r) - apiadu(l-1)) % p;
		((res %= p) += p) %= p;
	}
	return res * fpow(4, p-2, p) % p;
}

int main(){
	read(p, n);
	inv2 = fpow(2, p-2, p), inv6 = fpow(6, p-2, p);
	Euler(5000000);
	for(int i = 1; i <= 5000000; i++)
		pres[i] = (pres[i-1] + 1ll * i * i % p * phi[i] % p) % p;
	printf("%lld\n", solve());
	return 0;
}

---

[51nod 1220]约数之和

题目链接

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^n\sigma(ij)&=\sum_{i=1}^n\sum_{j=1}^n\sum_{a\mid i}\sum_{b\mid j}[\gcd(a,b)=1]\frac{aj}{b}\\ &=\sum_{a=1}^n\sum_{b=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{a}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{b}\right\rfloor}[\gcd(a,b)=1]aj\\ &=\sum_{a=1}^n\sum_{b=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{a}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{b}\right\rfloor}aj\sum_{d\mid\gcd(a,b)}\mu(d)\\ &=\sum_{a=1}^n\sum_{b=1}^n\sum_{d\mid\gcd(a,b)}a\mu(d)\left\lfloor\frac{n}{a}\right\rfloor\frac{\left\lfloor\frac{n}{b}\right\rfloor\left(\left\lfloor\frac{n}{b}\right\rfloor+1\right)}{2}\\ &=\sum_{d=1}^nd\mu(d)\left(\sum_{a=1}^{\left\lfloor\frac{n}{d}\right\rfloor}a\left\lfloor\frac{n}{ad}\right\rfloor\right)\left(\sum_{b=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\frac{\left\lfloor\frac{n}{bd}\right\rfloor\left(\left\lfloor\frac{n}{bd}\right\rfloor+1\right)}{2}\right) \end{align}$$

注意到：

$$\sum\limits_{i=1}^ni\left\lfloor\frac{n}{i}\right\rfloor=\sum\limits_{i=1}^n\sigma(i)$$ $$\sum\limits_{i=1}^n\frac{\left\lfloor\frac{n}{i}\right\rfloor\left(\left\lfloor\frac{n}{i}\right\rfloor+1\right)}{2}=\sum_{i=1}^n\sum_{j=1}^{\left\lfloor\frac{n}{i}\right\rfloor}j=\sum_{j=1}^{n}\sum_{i=1}^{\left\lfloor\frac{n}{j}\right\rfloor}j=\sum_{j=1}^nj\left\lfloor\frac{n}{j}\right\rfloor=\sum_{j=1}^n\sigma(j)$$

所以，

$$原式=\sum_{d=1}^nd\mu(d)\left(\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sigma(i)\right)^2$$

杜教筛筛出 $d\mu(d)$ 的前缀和以及 $\sigma(i)$ 的前缀和之后，就可以数论分块了。

$d\mu(d)$ 前缀和：设 $f(n)=n\mu(n),\;g(n)=\text{id}(n)$，则 $(fg)(n)=\sum\limits_{d=1}^nd\mu(d)\frac{n}{d}=n\cdot(\mu1)(n)=n\cdot\epsilon(n)=\epsilon(n)$，根据杜教筛理论，有：

$$S_f(n)=\sum_{i=1}^n\epsilon(i)-\sum_{i=2}^n\text{id}(i)S_f\left(\left\lfloor\frac{n}{i}\right\rfloor\right)=1-\sum_{i=2}^niS_f\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$$

$\sigma(i)$ 前缀和：设 $f(n)=\sigma(n)\;,g(n)=\mu(n)$，则 $(f*g)(n)=\text{id}(n)$，根据杜教筛理论，有：

$$S_\sigma(n)=\sum_{i=1}^n\text{id}(i)-\sum_{i=2}^n\mu(i)S_\sigma\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$$

所以筛出莫比乌斯函数的前缀和就可以了。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N = 1000005;
const LL MOD = 1e9+7;

inline LL fpow(LL bs, LL idx){
	LL res = 1;
	bs %= MOD;
	while(idx){
		if(idx & 1)	(res *= bs) %= MOD;
		(bs *= bs) %= MOD;
		idx >>= 1;
	}
	return res;
}

int mu[N], sigma[N], g[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	mu[1] = 1, sigma[1] = 1, g[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			mu[i] = -1, sigma[i] = 1 + i, g[i] = 1 + i;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				mu[i * pList[j]] = 0;
				sigma[i * pList[j]] = sigma[i] / g[i] * (g[i] * pList[j] + 1);
				g[i * pList[j]] = (g[i] * pList[j] + 1);
				break;
			}
			else{
				mu[i * pList[j]] = -mu[i];
				sigma[i * pList[j]] = sigma[i] * (1 + pList[j]);
				g[i * pList[j]] = 1 + pList[j];
			}
		}
	}
}

unordered_map<LL, LL> muS, sigmaS, fS;
LL pre_muS[N], pre_sigmaS[N], pre_fS[N];
LL apiadu_mu(LL n){
	if(n <= 1000000)	return pre_muS[n];
	if(muS.find(n) != muS.end())	return muS[n];
	LL res = 1;
	for(LL l = 2, r; l <= n; l = r + 1){
		r = n / (n / l);
		res -= apiadu_mu(n / l) * (r - l + 1) % MOD;
		((res %= MOD) += MOD) %= MOD;
	}
	return muS[n] = res;
}
LL apiadu_sigma(LL n){
	if(n <= 1000000)	return pre_sigmaS[n];
	if(sigmaS.find(n) != sigmaS.end())	return sigmaS[n];
	LL res = n * (n + 1) / 2 % MOD;
	for(LL l = 2, r; l <= n; l = r + 1){
		r = n / (n / l);
		res -= (apiadu_mu(r) - apiadu_mu(l-1)) * apiadu_sigma(n / l) % MOD;
		((res %= MOD) += MOD) %= MOD;
	}
	return sigmaS[n] = res;
}
LL apiadu_f(LL n){
	if(n <= 1000000)	return pre_fS[n];
	if(fS.find(n) != fS.end())	return fS[n];
	LL res = 1;
	for(LL l = 2, r; l <= n; l = r + 1){
		r = n / (n / l);
		res -= (l + r) * (r - l + 1) / 2 % MOD * apiadu_f(n / l) % MOD;
		((res %= MOD) += MOD) %= MOD;
	}
	return fS[n] = res;
}

int main(){
	Euler(1000000);
	for(int i = 1; i <= 1000000; i++){
		pre_muS[i] = pre_muS[i-1] + mu[i];
		pre_sigmaS[i] = (pre_sigmaS[i-1] + sigma[i]) % MOD;
		pre_fS[i] = ((pre_fS[i-1] + i * mu[i]) % MOD + MOD) % MOD;
	}
	int n; scanf("%d", &n);
	LL res = 0;
	for(LL l = 1, r; l <= n; l = r + 1){
		r = n / (n / l);
		res += (apiadu_f(r) - apiadu_f(l-1) + MOD) % MOD * apiadu_sigma(n / l) % MOD * apiadu_sigma(n / l) % MOD;
		((res %= MOD) += MOD) %= MOD;
	}
	printf("%lld\n", res);
	return 0;
}

---

[hdu6706] huntian oy

题目链接

$$\begin{align} f(n,a,b)&=\sum_{i=1}^n\sum_{j=1}^i\gcd(i^a-j^a,i^b-j^b)[\gcd(i,j)=1]\\ &=\sum_{i=1}^n\sum_{j=1}^i(i-j)[\gcd(i,j)=1]&&\gcd(i,j)=1\implies \gcd(i^a-j^a,i^b-j^b)=i-j\\ &=\sum_{i=1}^ni\sum_{j=1}^i[\gcd(i,j)=1]-\sum_{i=1}^n\sum_{j=1}^ij[\gcd(i,j)=1]\\ &=\sum_{i=1}^ni\varphi(i)-\sum_{i=1}^n\frac{i\varphi(i)+[i=1]}{2}&&使用定义\;\varphi(n)=\sum_{i=1}^n[\gcd(i,n)=1]\\ &&&和定理：小于等于\,n\,和\,n\,互质的数之和为\,\frac{n\varphi(n)+[n=1]}{2}\\ &=\frac{\sum\limits_{i=1}^ni\varphi(i)-1}{2} \end{align}$$

设 $f(n)=n\varphi(n)$，$g(n)=n$，则 $(f*g)(n)=\sum\limits_{d\mid n}d\varphi(d)\frac{n}{d}=n\sum\limits_{d\mid n}\varphi(d)=n^2$，根据杜教筛理论：

$$S_f(n)=\sum_{i=1}^ni^2-\sum_{i=2}^niS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$$

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N = 1000005;
const LL MOD = 1e9+7;
const LL inv2 = 500000004;
const LL inv3 = 333333336;
const LL inv6 = 166666668;

int phi[N], pList[N], pID;
bool notP[N];
void Euler(int n){
	notP[0] = notP[1] = 1;
	phi[1] = 1;
	for(int i = 1; i <= n; i++){
		if(notP[i] == 0){
			pList[++pID] = i;
			phi[i] = i - 1;
		}
		for(int j = 1; j <= pID; j++){
			if(1ll * i * pList[j] > n)	break;
			notP[i * pList[j]] = 1;
			if(i % pList[j] == 0){
				phi[i * pList[j]] = phi[i] * pList[j];
				break;
			}
			else	phi[i * pList[j]] = phi[i] * (pList[j] - 1);
		}
	}
}

unordered_map<LL, LL> fS;
LL preS[N];
LL apiadu_f(LL n){
	if(n <= 1000000)	return preS[n];
	if(fS.find(n) != fS.end())	return fS[n];
	LL res = n % MOD * (n + 1) % MOD * (2 * n % MOD + 1) % MOD * inv6 % MOD;
	for(LL l = 2, r; l <= n; l = r + 1){
		r = n / (n / l);
		res -= (l + r) % MOD * (r - l + 1) % MOD * inv2 % MOD * apiadu_f(n / l) % MOD;
		((res %= MOD) += MOD) %= MOD;
	}
	return fS[n] = res;
}

int main(){
	Euler(1000000);
	for(int i = 1; i <= 1000000; i++)
		preS[i] = (preS[i-1] + 1ll * phi[i] * i % MOD) % MOD;
	int T; for(scanf("%d", &T); T; T--){
		LL n, res;
		scanf("%lld%*d%*d", &n);
		res = (apiadu_f(n) - 1) % MOD * inv2 % MOD;
		((res %= MOD) += MOD) %= MOD;
		printf("%lld\n", res);
	}
	return 0;
}