LOJ 6053. 简单的函数 Min_25筛

传送门：https://loj.ac/p/6053

题意

有一个神奇的函数f(x)，满足以下性质：

  • $f(1)=1$
  • $f(p^c)=p \oplus c$
  • $f(ab)=f(a)f(b)\;gcd(a,b)=1的情况下$

输出$\sum_{i=1}^nf(i)\;mod\;1e9+7$

思路

因为n有1e10,所有用线性筛会T。 只能用亚线性筛，而f明显是一个积性函数，可以洲阁筛，当然，Min_25更完爆它。

那么怎么构造出一个Min_25筛呢？转移方程这么写？

观察$f(p^c)=p \oplus c$,对于一个质数p:

  1. $if(p=2) \;f(p)=3$
  2. $else\;f(p)=p-1$

所以我们要处理质数和，设$g1[j]=\sum_{i=2}^j1,g2[j]=\sum_{i=2}^ji.$ 不管第一个f(2)=3，就统一认为f(p)=p-1，最后如果有2的话，ans+2即可。

设prs[j]为前j个质数和，prz[j]为前j个质数个数。

g函数状态转移： $$g_1(j,n)=g(j-1,n)-[g(j-1,\frac{n}{p_j})-prz[j-1]]$$

$$g_2(j,n)=g(j-1,n)-pr[j]*(g[j-1,\frac{n}{p_j}]-prs[j-1])$$

S(x,y)的状态转移： $$S(x,y)=\sum{g_t[id(x)]-sum_t[j-1]}+\sum_{kge y\;p_k^e\leq x}F(p_k^e)*S(\frac{x}{p_k^e},k+1)+F(p_k^{e+1})$$

$$S(x,y)=\sum{g_t[id(x)]-sum_t[j-1]}+\sum_{kge y\;p_k^e\leq x}(p \oplus e)*S(\frac{x}{p_k^e},k+1)+p \oplus (e+1)$$

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
##include "bits/stdc++.h"
using namespace std;

typedef long long ll;
const ll mod = 1e9 + 7;

const int N = 1e6 + 10;

ll quick_pow(ll a, ll b) {
    ll ans = 1;
    while(b) {
        if(b & 1) ans = ans * a % mod;
         a = a * a % mod;
         b >>= 1;
    }
    return ans % mod;
}

namespace Min_25 {
    int is_pr[N], id1[N], id2[N], cnt;
    ll n, m, T;
    ll pr[N], prs[N], a[N], g1[N], g2[N];

    ll inv2;

    inline ll ID(ll x) {
        return x <= T ? id1[x] : id2[n / x];
    }

    inline ll calc1(ll x) {
        return (x - 1) % mod;
    }

    inline ll calc2(ll x) {
        return x % mod * (x + 1) % mod * inv2 % mod - 1;
    }

    void init() {
        cnt = m = 0;
        inv2 = quick_pow(2, mod - 2);
        T = sqrt(n + 0.5);
        for(int i = 2;i <= T; i++) {
            if(!is_pr[i]) pr[++cnt] = i, prs[cnt] = (prs[cnt - 1] + i) % mod;
            for(int j = 1;j <= cnt && i * pr[j] <= T; j++) {
                is_pr[i * pr[j]] = true;
                if(i % pr[j] == 0) break;
            }
        }

        for(ll l = 1;l <= n; l = n / (n / l) + 1) {
            a[++m] = n / l;
            if(a[m] <= T) id1[a[m]] = m;
            else id2[n / a[m]] = m;
            g1[m] = calc1(a[m]);
            g2[m] = calc2(a[m]);
        }

        for(int i = 1;i <= cnt; i++) {
            for(int j = 1;j <= m && pr[i] * pr[i] <= a[j]; j++) {
                g1[j] = (g1[j] - (g1[ID(a[j] / pr[i])] - (i - 1))) % mod;
                g2[j] = (g2[j] - pr[i] * (g2[ID(a[j] / pr[i])] - prs[i - 1]) % mod) % mod;
            }
        }
    }

    ll S(ll x, ll y) {
        if(x <= 1  x < pr[y]) return 0;
        ll ans = (g2[ID(x)] - prs[y - 1] - (g1[ID(x)] - (y - 1)) + (y == 1 ? 2 : 0) + mod) % mod;
        for(int i = y;i <= cnt && pr[i] * pr[i] <= x; i++) {
            ll pe = pr[i];
            for(int e = 1;pe * pr[i] <= x; e++, pe *= pr[i]) {
                ans = (ans + (pr[i] ^ e) * S(x / pe, i + 1) % mod + (pr[i] ^ e + 1) % mod) % mod;
            }
        }
        return ans % mod;
    }

    ll Solve(ll x) {
        return n = x, init(), S(n, 1) + 1;
    }
}

void solve() {
    ll n; cin >> n;
    ll ans = Min_25::Solve(n);
    cout << (ans % mod + mod) % mod << endl;
}

signed main() {
    solve();
}

本文作者：jujimeizuo
本文地址： https://blog.jujimeizuo.cn/2021/03/05/loj-6053/
本博客所有文章除特别声明外，均采用 CC BY-SA 3.0 协议。转载请注明出处！