51nod1847 奇怪的数学题

发表于 | 阅读数

杜教筛一下

  • 其实杜教筛也只有$\sqrt{n}$个状态，并且发现在进行数论分块的时候$n/i,与n/(n/i)$，分布是相同的。所以提前杜娇筛，那么$F(x),x\in{n/i}$都可以得出来。前面$\sum_{i=1}^x f(x)$，再Min_25筛的时候正好存储了$\sum_{i=1}^x f(x),x\in{n/i}$
  • 无法处理$2^{32}$的逆元，又$n+1\rightarrow(n-i+1)$一定存在$\mod i+1 =0$。特胖即可，注意$n<i=0$
  • 由于质数较为特殊，拿出来处理$\sum_i^{Prime} f(i)=\sum[i\in Prime]$ ,再搞一个Min_25处理下即可。

时间复杂度$O(k^2\sqrt{n}+\frac{n^{0.75}}{\log n}+n^{\frac{2}{3}})$

代码 
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned int ui;
#define pii pair<int, int>
#define mk make_pair
const int N = 1e6 + 10;
const int mod = 1e9 + 7;

ui qpow(ui a, ui x)
{
    ui res = 1;
    while (x)
    {
        if (x & 1)
            res = 1ll * res * a;
        a = 1ll * a * a;
        x >>= 1;
    }
    return res;
}
ui pcnt, pr[N], npr[N];

void Prime_init(int X)
{
    npr[1] = 1;
    for (int i = 1; i <= X; i++)
    {
        if (!npr[i])
            pr[++pcnt] = i;
        for (int j = 1; j <= pcnt && pr[j] * i <= X; j++)
        {
            npr[pr[j] * i] = 1;
            if (i % pr[j] == 0)
                break;
        }
    }
}
ui S[60][60];
int cnt;
ui n;
ui p[N], pos1[N], pos2[N], pk[N], spk[N], g0[N], gk[N], sp[N];
ui Lim;
ui m;
int get(int x)
{
    return (x <= Lim) ? pos1[x] : pos2[n / x];
}
ui calc(ui x)
{
    ui ans = 0;
    for (int i = 1; i <= m && i <= x; i++)
    {
        ui res = 1;
        for (ui j = x + 1; j >= (x + 1 - (i + 1) + 1); j--)
        {
            if (j % (i + 1) == 0)
                res *= (j) / (i + 1);
            else
                res *= j;
        }
        ans += res * S[m][i];
    }
    return ans;
}
void init()
{
    Lim = sqrt(n);
    Prime_init(Lim);

    S[0][0] = 1;
    for (int i = 1; i < 55; i++)
        for (int j = 1; j <= i; j++)
            S[i][j] = S[i - 1][j - 1] + S[i - 1][j] * j;
    for (int i = 1; i <= pcnt; i++)
        pk[i] = qpow(pr[i], m), spk[i] = spk[i - 1] + pk[i];

    for (ui i = 1, j; i <= n; i = j + 1)
    {

        j = n / (n / i);
        ui k = n / i;
        p[++cnt] = k;
        if (k <= Lim)
            pos1[k] = cnt;
        else
            pos2[n / k] = cnt;
        g0[cnt] = k - 1;
        gk[cnt] = calc(k) - 1;
    }

    for (int i = 1; i <= pcnt; i++)
    {
        for (int j = 1; j <= cnt && pr[i] * pr[i] <= p[j]; j++)
        {
            int k = get(p[j] / pr[i]);
            g0[j] -= g0[k] - (i - 1);
            gk[j] -= pk[i] * (gk[k] - spk[i - 1]);
            sp[j] += (gk[k] - spk[i - 1]);
        }
    }
}
int vis[N];
ui sum[N];
ui dyh(ui x)
{
    if (x <= 0)
        return 0;
    int k = get(x);
    if (vis[k])
        return sum[k];
    ui res = sp[k] + g0[k];
    for (ui i = 2, j; i <= x; i = j + 1)
    {
        j = x / (x / i);
        res -= (j - i + 1) * dyh(x / i);
    }
    vis[k] = 1;
    return sum[k] = res;
}

int main()
{
    cin >> n >> m;
    init();
    ui ans = 0;

    dyh(n);
    for (ui i = 1, j; i <= n; i = j + 1)
    {

        j = n / (n / i);
        ans += (n / i) * (n / i) * (sum[get(j)] - sum[get(i - 1)]);
    }
    cout << ans << endl;
}