题解 P9091 积性函数求和

· · 题解

题目传送门

前言

尝试介绍下一种容易实现且常数很小的积性函数求和法。

Solution

式子的推导请见出题人 Leasier 的题解，进一步转化后得到原式 =mS_{I_{m+1}}(n)，其中 I_n 表示 I*I\,*···*\,I (卷 n 次)，为积性函数。

对于积性函数 f，构造 g 满足 g(p^k)=f(p)^k，则 f=g*h，易得 h(p^k)=f(p^k)-f(p)f(p^{k-1})。因此我们可以先算较容易求出的 g，再借助 PN 得到 f。

回忆一下 min25 筛的第一部分，g 为满足 g(p^k)=g(p)^k 的任一积性函数，令 $S(n, a)=\sum\limits{i=1}^{n}[min_i>p_a \,or\;i\in P]\,g(i)$，有 :

S(n,a)=S(n,a-1)-g(p_a)(S(\frac{n}{p_a},a-1)-S(p_{a-1},a-1))

那么其实它也能写成这样 :

S(n,a-1)=S(n,a)+g(p_a)(S(\frac{n}{p_a},a-1)-S(p_{a-1},a-1))

也就是说对于我们要求的 g，我们只需先求出其在素数处的和（即得到 S(m,\pi(\sqrt{n}))，m=\left\lfloor\dfrac{n}{i}\right\rfloor，i=1,2,···,\sqrt{n}），便可通过后一式得到 g 的块筛。

之后就可以在 O(\sqrt n) 内得到 S_f(n)，或在 O(\sqrt{n}\log^{1+\varepsilon} n) 内得到 f 的块筛。

实现

远古代码直接粘贴过来了，欠佳。

#include <iostream>
#include <vector>
#include <cmath>
#include <functional>
#include <algorithm>
#include <time.h>
#include <fstream>
using namespace std;
using u32 = unsigned int;
using u64 = unsigned long long;
using i64 = long long;
double inv[100005];
u32 c[70][70];
vector<int> primes;
u32 solve(i64 N, u32 k) {
    int v = log2(N) + k;

    const auto div = [&](i64 n, int d) -> i64 {return n * inv[d]; };
    k++;
    v = sqrt(N + 0.5);
    for (int i = 1; i <= v; ++i) inv[i] = 1.0 / i * (1 + 1e-15);
    vector<u32> s(v + 1), l(v + 1);
    for (int i = 1; i <= v; ++i) s[i] = i - 1, l[i] = div(N, i) - 1;
    for (int p = 2; p <= v; ++p)
        if (s[p] != s[p - 1]) {
            primes.push_back(p);
            i64 q = i64(p) * p, M = div(N, p);
            u32 t0 = s[p - 1];
            int t = v / p, u = min<i64>(v, div(M, p));
            for (int i = 1; i <= t; ++i) l[i] -= l[i * p] - t0;
            for (int i = t + 1; i <= u; ++i) l[i] -= s[div(M, i)] - t0;
            for (int i = v, j = t; j >= p; --j)
                for (int l = j * p; i >= l; --i)
                    s[i] -= s[j] - t0;
        }
    for (int i = 1; i <= v; ++i) s[i] *= k, l[i] *= k;
    for (auto it = primes.rbegin(); it != primes.rend(); ++it) {
        int p = *it;
        i64 q = i64(p) * p, M = div(N, p);
        int t = v / p, u = min<i64>(v, div(M, p));
        u32 t0 = s[p - 1];
        if(q <= v) {
            u32 s0 = k * k;
            for (int j = p, i = q; j < t; s0 = (s[++j] - t0) * k)
                for (int l = j * p + p; i < l; ++i) s[i] += s0;
            for (int i = t * p; i <= v; ++i) s[i] += s0;
        }
        for (int i = u; i > t; --i) l[i] += (s[div(M, i)] - t0) * k;
        for (int i = t; i >= 1; --i) l[i] += (l[i * p] - t0) * k;
    }
    for (int i = 1; i <= v; ++i) s[i] += 1, l[i] += 1;

    k--;
    const auto divide = [](i64 n, i64 d) -> i64 {return double(n) / d; };
    function< u32(i64, int, u32)> rec = [&] (i64 n, size_t beg, u32 coeff) -> u32 {
        if (!coeff) return 0;
        u32 ret = coeff * (n > v ? l[divide(N, n)] : s[n]);
        for (size_t i = beg; i < primes.size(); ++i) {
            int p = primes[i];
            i64 q = i64(p) * p;
            if (q > n) break;
            i64 nn = divide(n, q);
            for (u32 e = 1; nn; nn = div(nn, p), ++e)
                ret += rec(nn, i + 1, -coeff * e * c[e + k][e + 1]);
        }
        return ret;
    };
    return k * rec(N, 0, 1);
}
int main() {
    i64 N;
    u32 k;
    cin >> N >> k;
    cout << solve(N, k);
    return 0;
}