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}^nsum_{j=1}^nlcm(i,j)
]
前置题目是Carsh的数字表格和jzptab（BZOJ都有，第二个是权限题）
即本题的(O(n+nsqrt{n}))做法和(O(n+sqrt{n}))做法。
题1之前写了题解。
这里来讲(O(n^{frac{2}{3}}))做法（使用杜教筛，具体复杂度其实并不是很清楚。。。不会证明）
用莫比乌斯反演套路可以推到这里
[egin{aligned}
&设sum(n)=sum_{i=1}^ni\
&sum_{T=1}^nsum(frac{n}{T})^2Tsum_{k|T}mu(k)k
end{aligned}
]
重点研究如何在非线性复杂度下筛出后面一段
[egin{aligned}
&定义数论函数的点积为逐项相乘，符号为·\
&那么将Tsum_{k|T}mu(k)k化为sum_{k|T}mu(k)k^2frac{T}{k}\
&会发现这是个狄利克雷卷积和点积混在一起的形式（*为狄利克雷卷积）\
&即(id^2·mu)*id
end{aligned}
]
考虑如何选取一个合适的g函数让它能杜教筛筛出来
[egin{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('
')
#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
", (ans % mod + mod) % mod);
}