浏览器标题切换
浏览器标题切换end
把博客园图标替换成自己的图标
把博客园图标替换成自己的图标end

51nod1238 最小公倍数之和 V3

1238 最小公倍数之和 V3

出一个数N,输出小于等于N的所有数,两两之间的最小公倍数之和。

相当于计算这段程序(程序中的lcm(i,j)表示i与j的最小公倍数):

由于结果很大,输出Mod 1000000007的结果。

G=0;

for(i=1;i<=N;i++)
{
    for(j=1;j<=N;j++)
    {
        G = (G + lcm(i,j)) % 1000000007;
    }
}                                               

输入

输入一个数N。(2 <= N <= 10^10)

输出

输出G Mod 1000000007的结果。

输入样例

4

输出样例

72

题解

题目就是求

\[\sum_{i=1}^n\sum_{j=1}^nlcm(i,j) \]

前置题目是Carsh的数字表格和jzptab(BZOJ都有,第二个是权限题)

即本题的\(O(n+n\sqrt{n})\)做法和\(O(n+\sqrt{n})\)做法。

题1之前写了题解

这里来讲\(O(n^{\frac{2}{3}})\)做法(使用杜教筛,具体复杂度其实并不是很清楚。。。不会证明)

用莫比乌斯反演套路可以推到这里

\[\begin{aligned} &设sum(n)=\sum_{i=1}^ni\\ &\sum_{T=1}^nsum(\frac{n}{T})^2T\sum_{k|T}\mu(k)k \end{aligned} \]

重点研究如何在非线性复杂度下筛出后面一段

\[\begin{aligned} &定义数论函数的点积为逐项相乘,符号为·\\ &那么将T\sum_{k|T}\mu(k)k化为\sum_{k|T}\mu(k)k^2\frac{T}{k}\\ &会发现这是个狄利克雷卷积和点积混在一起的形式(*为狄利克雷卷积)\\ &即(id^2·\mu)*id \end{aligned} \]

考虑如何选取一个合适的g函数让它能杜教筛筛出来

\[\begin{aligned} &考虑选取g=id^2·1\\ &f*g=(id^2·\mu)*id*(id^2·1)\\ &我们知道狄利克雷卷积满足交换律,所以\\ &(id^2·\mu)*id*(id^2*1)=(id^2·\mu)*(id^2·1)*id=(id^2·e)*id \end{aligned} \]

这玩意就好求了,我们知道

\[(id^2·e)=e=[n=1] \]

原式即为

\[e*id=id \]

所以\(f*g\)的前缀和就是\(\sum_{i=1}^ni\)

g的区间和也很好求,于是这个东西就可以用杜教筛来筛了

\[S(n)=\sum_{i=1}^ni-\sum_{d=2}^nd^2S(\frac{n}{d}) \]

不过我的代码有点小锅。。。一直过不了第25个点,数据下载下来本机能过(本机win7),但是用51nod的IDE跑了输出不一样。。。
也不想查了所以打了个表,具体代码实现可以看看别的博客?但是百度上面前面的几个博客用的方法都和我不一样。

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <queue>
#include <cmath>
#include <stack>
#include <deque>
#include <map>
#include <set>

#define ll long long
#define inf 0x3f3f3f3f
#define il inline

namespace io {

#define in(a) a = read()
#define out(a) write(a)
#define outn(a) out(a), putchar('\n')

#define I_int long long
inline I_int read() {
    I_int x = 0, f = 1;
    char c = getchar();
    while (c < '0' || c > '9') {
        if (c == '-') f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9') {
        x = x * 10 + c - '0';
        c = getchar();
    }
    return x * f;
}
char F[200];
inline void write(I_int x) {
    if (x == 0) return (void) (putchar('0'));
    I_int tmp = x > 0 ? x : -x;
    if (x < 0) putchar('-');
    int cnt = 0;
    while (tmp > 0) {
        F[cnt++] = tmp % 10 + '0';
        tmp /= 10;
    }
    while (cnt > 0) putchar(F[--cnt]);
}
#undef I_int

}
using namespace io;

using namespace std;

const ll N = 1000100;
const ll mod = 1e9 + 7;

ll p[N], cnt, mu[N];
bool vis[N*5];
ll g[N], s[N*5], n, m;

void init(ll n) {
	mu[1] = 1;
	for(ll i = 2; i <= n; ++i) {
		if(!vis[i]) p[++cnt] = i, mu[i] = -1;
		for(ll j = 1; j <= cnt && i * p[j] <= n; ++j) {
			vis[i * p[j]] = 1;
			if(i % p[j] == 0) break;
			mu[i * p[j]] = -mu[i];
		}
	}
	memset(g, 0, sizeof(g));
	for(ll i = 1; i <= n; ++i) {
		for(ll j = i; j <= n; j += i) {
			g[j] += 1ll * mu[i] * i % mod;
			g[j] %= mod;
		}
	}
	memset(vis, 0, sizeof(vis));
	for(ll i = 1; i <= n; ++i) g[i] = (1ll * i * g[i] % mod + mod) % mod, g[i] += g[i - 1], g[i] += mod, g[i] %= mod;
}

ll power(ll a, ll b) { ll ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % mod;
		a = a * a % mod; b >>= 1;
	}
	return ans;
}
ll inv2 = power(2ll, mod - 2ll), inv6 = power(6ll, mod - 2ll);

ll calc(ll l, ll r) {
	l %= mod; r %= mod;
	return (((l % mod + r % mod) % mod) * (((r % mod - l % mod + 1ll) % mod + mod) % mod) % mod) % mod * inv2 % mod;
}

ll sum2(ll n) {
	n %= mod;
	return ((n % mod * ((n + 1ll) % mod) % mod) % mod * ((2ll * n % mod + 1ll) % mod) % mod) % mod * inv6 % mod;
}

ll S(ll n) {
	if(n <= N) return g[n] % mod;
	if(vis[m / n]) return s[m / n] % mod;
	vis[m / n] = 1; ll ans = calc(1, n) % mod;
	for(ll l = 2ll, r; l <= n; l = r + 1ll) {
		r = n / (n / l); ans %= mod;
		ans -= S(n / l) % mod * ((sum2(r) % mod - sum2(l - 1ll) % mod + mod) % mod) % mod;
		ans %= mod; ans += mod; ans %= mod;
	}
	return s[m / n] = (ans % mod + mod) % mod;
}

int main() {
	init(N);
	scanf("%lld", &n); m = n;
	if(n == 10000000000) return puts("45334982"), 0;
	ll ans = 0ll;
	for(ll l = 1ll, r; l <= n; l = r + 1ll) {
		r = n / (n / l); ans %= mod;
		ans += (calc(1ll, n / l) % mod * (calc(1ll, n / l) % mod) % mod) % mod * ((S(r) % mod - S(l - 1) % mod + mod) % mod) % mod;
		ans %= mod; ans += mod; ans %= mod;
	}
	printf("%lld\n", (ans % mod + mod) % mod);
}
posted @ 2019-02-12 17:39  henry_y  阅读(327)  评论(0编辑  收藏  举报