51NOD1847:奇怪的数学题

传送门

Sol

\(f(d)\) 表示 \(d\) 所有约数中第二大的,\(low_d\) 表示 \(d\) 的最小质因子
\[f(d)=\frac{d}{low_d}\]
那么
\[\sum_{i=1}^{n}\sum_{j=1}^{n}sgcd^k(i,j)\]
\[=\sum_{i=1}^n\sum_{j=1}^{n}f^k(gcd(i,j))\]
\[=\sum_{d=1}^{n}f^k(d)\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=d]\]
\[=\sum_{d=1}^{n}f^k(d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}[gcd(i,j)=1]\]
\[=\sum_{d=1}^{n}f^k(d)(2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\varphi(i)-1)\]
可以数论分块,后面的 \(\varphi\) 直接杜教筛
考虑计算
\[\sum_{d=2}^{n}f^k(d)=\sum_{d=2}^{n}\frac{d}{low_d}^k\]
\(p_j\) 表示第 \(j\) 个质数
\(g(x,i)\) 表示 \(2\)\(x\) 之间最小质因子大于等于 \(p_i\) 的或者质数的 \(f\)\(k\) 次方和
\(g'(x,i)\) 表示 \(2\)\(x\) 之间最小质因子大于等于 \(p_i\) 的或者质数的 \(k\) 次方和
\(s(i)\) 表示小于等于 \(p_i\) 的质数的 \(k\) 次方和
那么就是要求 \(g(n, 1)\)
\(s\) 直接 \(min25\)
\[g(x,i)=g(x,i+1)+\sum_{e=1}^{p_i^{e+1}\le x}p_i^{k(e-1)}(g'(\lfloor\frac{x}{p_i^{e}}\rfloor,i+1)-s(i)+p_i^{k})\]
\[g'(x,i)=g'(x,i+1)+\sum_{e=1}^{p_i^{e+1}\le x}p_i^{ke}(g'(\lfloor\frac{x}{p_i^{e}}\rfloor,i+1)-s(i)+p_i^{k})\]
这里要用到自然幂数和求和,用第二类斯特林数就好了

# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned int uint;

inline uint Pow(uint x, int y) {
    register uint ret = 1;
    for (; y; y >>= 1, x = x * x) if (y & 1) ret = ret * x;
    return ret;
}

const int maxn(1e6 + 5);

int pr[maxn], tot, k, d, id1[maxn], id2[maxn], cnt, phi[maxn];
uint n, s[55][55], f[maxn], g[maxn], val[maxn], sk[maxn];
uint s1[maxn], s2[maxn], sphi[maxn], sp[maxn];
bitset <maxn> ispr;

inline void Sieve(int mx) {
    register int i, j;
    for (phi[1] = 1, ispr[1] = 1, i = 2; i <= mx; ++i) {
        if (!ispr[i]) pr[++tot] = i, sk[tot] = sk[tot - 1] + Pow(i, k), phi[i] = i - 1;
        for (j = 1; j <= tot && pr[j] * i <= mx; ++j) {
            ispr[i * pr[j]] = 1;
            if (i % pr[j]) phi[i * pr[j]] = phi[i] * (pr[j] - 1);
            else {
                phi[i * pr[j]] = phi[i] * pr[j];
                break;
            }
        }
    }
    for (i = 1; i <= mx; ++i) sphi[i] = sphi[i - 1] + phi[i];
}

# define ID(x) ((x) <= d ? id1[x] : id2[n / (x)])

inline uint Sum(uint x) {
    register uint i, j, v = 0, t, r, tmp = k <= x ? k : x;
    for (i = 1; i <= tmp; ++i) {
        t = i + 1, r = s[k][i];
        for (j = x - i + 1; j <= x + 1; ++j)
            if (t > 1 && j % t == 0) r *= j / t, t = 1;
            else r = r * j;
        v += r;
    }
    return v;
}

uint Sumphi(uint x) {
    if (x <= d) return sphi[x];
    if (sp[ID(x)]) return sp[ID(x)];
    register uint ans = (x & 1) ? ((x + 1) >> 1) * x : (x >> 1) * (x + 1), i, j;
    for (i = 2; i <= x; i = j + 1) j = x / (x / i), ans -= Sumphi(x / i) * (j - i + 1);
    return sp[ID(x)] = ans;
}

int main() {
    register uint i, j, e, ans = 0, lst = 0, cur, r, v, tmp, now;
    scanf("%u%d", &n, &k), Sieve(d = sqrt(n));
    for (i = 1; i <= 50; ++i)
        for (s[i][1] = 1, j = 2; j <= i; ++j)
            s[i][j] = s[i - 1][j] * j + s[i - 1][j - 1];
    for (i = 1; i <= n; i = j + 1) {
        val[++cnt] = n / i, j = n / (n / i);
        val[cnt] <= d ? id1[val[cnt]] = cnt : id2[n / val[cnt]] = cnt;
        f[cnt] = val[cnt] - 1, g[cnt] = Sum(val[cnt]) - 1;
    }
    for (i = 1; i <= tot && pr[i] * pr[i] <= n; ++i)
        for (j = 1; j <= cnt && pr[i] * pr[i] <= val[j]; ++j) {
            f[j] -= f[ID(val[j] / pr[i])] - i + 1;
            g[j] -= (sk[i] - sk[i - 1]) * (g[ID(val[j] / pr[i])] - sk[i - 1]);
        }
    for (i = 1; i <= cnt; ++i) s1[i] = f[i], s2[i] = g[i];
    for (r = 1; r <= tot && pr[r] * pr[r] <= n; ++r);
    for (i = r - 1; i; --i)
        for (tmp = sk[i] - sk[i - 1], j = 1; j <= cnt && pr[i] * pr[i] <= val[j]; ++j)
            for (cur = e = 1, v = pr[i]; v <= val[j] / pr[i]; ++e, v *= pr[i], cur *= tmp)
                now = (s2[ID(val[j] / v)] - sk[i] + tmp) * cur, s1[j] += now, s2[j] += tmp * now;
    for (i = 1; i <= n; i = j + 1) {
        j = n / (n / i), cur = s1[ID(j)];
        ans += (cur - lst) * (Sumphi(n / i) * 2 - 1), lst = cur;
    }
    printf("%u\n", ans);
    return 0;
}

转载于:https://www.cnblogs.com/cjoieryl/p/10162187.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值