ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

【杜教筛】这是一道简单的数学题

2022-01-01 15:37:19  阅读:160  来源: 互联网

标签:frac limits kd int sum 杜教 mu 一道 数学题


题目大意:

\[\sum\limits_{i=1}^n\sum\limits^i_{j=1}\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)} \]

\(1\le n\le 10^9\)。

题解:

为了方便,考虑求 \(\sum\limits_{i=1}^n\sum\limits^n_{j=1}\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\)

\[\sum\limits_{i=1}^n\sum\limits^n_{j=1}\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\\ =\sum\limits_{i=1}^n\sum\limits^n_{j=1}\frac{ij}{\gcd^2(i,j)}\\ =\sum\limits_{d=1}\frac{\sum\limits^n_{i=1}\sum\limits^n_{j=1}ij[\gcd(i,j)=d]}{d^2}\\ =\sum\limits_{d=1}\frac{\sum\limits_{k=1}\mu(k)(\frac{n}{kd})^2(\frac{n}{kd}+1)^2(kd)^2/4}{d^2}\\ =\sum\limits_{d=1}\sum\limits_{k=1}\mu(k)(\frac{n}{kd})^2(\frac{n}{kd}+1)^2k^2/4\\ \]

以下忽略那个 \(/4\),因为这和整个式子毫无关联。

考虑枚举 \(kd\),得到原式\(=\)

\[\sum\limits_{d=1}(\frac{n}{d})^2(\frac{n}{d}+1)^2\sum_{k\vert d}\frac{\mu(k)k^2}{d} \]

\[\sum\limits_{d=1}(\frac{n}{d})^2(\frac{n}{d}+1)^2\frac{\sum_{k\vert d}\mu(k)k^2}{d} \]

最后一步,考虑使用杜教筛快速求出函数 \(f(n)=\sum_{k\vert d}\mu(k)k^2\) 的前缀和。

令 \(h(n)=\mu(n)n^2\),那么 \(f(n)=h*1\)。通过观察,发现若构造函数 \(g(n)=n^2\),则有 \(h*g=\epsilon\)。那么 \(f*g=h*g*1=\epsilon*1=1\)。

综上,我们可以快速求出函数 \(g\) 与函数 \(f*g\) 的前缀和,套杜教筛即可。

#include <cstdio>
#include <cmath>

const int mod = 1e9 + 7, inv2 = 5e8 + 4, inv6 = 166666668;
int f[1500005], Prime[150005], cnt, n;
bool mark[1500005];
class hash_map {
	static const int mod = 19260817;
	int key[10000005], val[10000005], head[10000005], nxt[10000005], tot;
public:
	inline void insert(int x, int y) {
		key[++ tot] = x, val[tot] = y, nxt[tot] = head[x % mod], head[x % mod] = tot;
	}
	inline int find(int x) const {
		for (int i = head[x % mod]; i; i = nxt[i])
			if (key[i] == x) return val[i];
		return -2000000000;
	}
} Sum;

void init(int n) {
	f[1] = 1;
	Sum.insert(0, 0), Sum.insert(1, 1);
	for (int i = 2, sum = 1; i <= n; ++ i) {
		if (!mark[i]) f[i] = (1ll - 1ll * i * i % mod) % mod, Prime[++ cnt] = i;
		Sum.insert(i, sum = (sum + f[i]) % mod);
		for (int j = 1; j <= cnt && i * Prime[j] <= n; ++ j) {
			mark[i * Prime[j]] = true;
			if (i % Prime[j] == 0) {f[i * Prime[j]] = f[i]; break;}
			f[i * Prime[j]] = 1ll * f[i] * f[Prime[j]] % mod;
		}
	}
}

inline long long sq(int n) {return 1ll * n * n % mod * (n + 1) % mod * (n + 1) % mod;}
inline int g(int n) {return 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;}
int getsum(int n) {
	if (Sum.find(n) != -2000000000) return Sum.find(n);
	int ans = n;
	for (int l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans = (ans - 1ll * (g(r) - g(l - 1)) % mod * getsum(n / l) % mod) % mod;
	}
	Sum.insert(n, ans);
	return ans;
}

int main() {
	init(1500000);
	int n, ans = 0;
	scanf("%d", &n);
	for (int l = 1, r; l <= n; l = r + 1) {	
		r = n / (n / l);
		ans = (ans + sq(n / l) % mod * (getsum(r) - getsum(l - 1)) % mod) % mod;
	}
	printf("%d", ((int)((1ll * ans * inv2 % mod * inv2 % mod + n) % mod * inv2 % mod) + mod) % mod);
	return 0;
}

标签:frac,limits,kd,int,sum,杜教,mu,一道,数学题
来源: https://www.cnblogs.com/stinger/p/15755445.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有