又被这神仙题给坑爆了.
神仙题解.
一开始我把 lcm 变成 ij/gcd 然后按照常规套路去推, 推到最后发现不是 miu * Id 而是 miu . Id...... 这还搞鬼啊.
正解居然跟这个差不多, 先转成求其中一部分的函数, 然后再加和...... 这谁顶得住呀.
大概就是先求这个
一顿操作之后得到了 phi 有关的式子......
然后原式就是这个
然后带进去推一推就出来杜教筛了... 这第一步真是神奇.
最后是这个.
按照套路, 前面分块, 后面配一个 g(x) = x2 即可.
- #include <cstdio>
- #include <map>
- typedef long long LL;
- const int N = 5000010, T = 5000008;
- const LL MO = 1000000007;
- std::map<LL, LL> mp;
- LL inv2, inv6, F[N];
- int p[N], top, phi[N];
- bool vis[N];
- inline LL qpow(LL a, LL b) {
- LL ans = 1;
- while(b) {
- if(b & 1) ans = ans * a % MO;
- a = a * a % MO;
- b = b>> 1;
- }
- return ans;
- }
- inline LL S(LL x) {
- x %= MO;
- return x * (x + 1) / 2 % MO;
- }
- inline LL G(LL x) {
- x %= MO;
- return (x << 1 | 1) % MO * (x + 1) % MO * x % MO * inv6 % MO;
- }
- inline LL H(LL x) {
- LL temp = S(x);
- return temp * temp % MO;
- }
- inline void getp(int n) {
- phi[1] = 1;
- for(int i = 2; i <= n; i++) {
- if(!vis[i]) {
- p[++top] = i;
- phi[i] = i - 1;
- }
- for(int j = 1; j <= top && i * p[j] <= n; j++) {
- vis[i * p[j]] = 1;
- if(i % p[j] == 0) {
- phi[i * p[j]] = phi[i] * p[j];
- break;
- }
- phi[i * p[j]] = phi[i] * (p[j] - 1);
- }
- }
- for(int i = 1; i <= n; i++) {
- F[i] = (F[i - 1] + 1ll * i * i % MO * phi[i] % MO) % MO;
- }
- return;
- }
- LL getF(LL x) {
- if(x <= 0) return 0;
- if(x <= T) return F[x];
- if(mp.count(x)) return mp[x];
- LL ans = H(x);
- for(LL i = 2, j; i <= x; i = j + 1) {
- j = x / (x / i);
- ans -= (G(j) - G(i - 1) + MO) * getF(x / i) % MO;
- ans %= MO;
- }
- return mp[x] = (ans + MO) % MO;
- }
- int main() {
- inv2 = (MO + 1) / 2;
- inv6 = qpow(6, MO - 2);
- getp(T);
- LL ans = 0, n;
- scanf("%lld", &n);
- for(LL i = 1, j; i <= n; i = j + 1) {
- j = n / (n / i);
- ans += S(n / i) * (getF(j) - getF(i - 1) + MO) % MO;
- ans %= MO;
- }
- printf("%lld\n", (ans + MO) % MO);
- return 0;
- }
AC 代码
来源: http://www.bubuko.com/infodetail-2970569.html