【题解】51nod-1965-奇怪的式子

题目

http://www.51nod.com/onlineJudge/problemSolution.html#!problemId=1965

题意

\(\prod_{i=1}^n \sigma_0(i)^{\mu(i)+i}\)

题解

\(S_1=\prod_{i=1}^n \sigma_0(i)^i,S_2=\prod_{i=1}^n \sigma_0(i)^{\mu(i)}\)

\(x=\prod p_i^{\alpha_i}\),那么有 \(\sigma_0(x) = \prod (\alpha_i+1)\)

枚举素数 \(p\) 的指数 \(k\) 的贡献,有 \[ \begin{eqnarray} S_1 &=& \prod_p \prod_k (k+1)^{\sum_{i=1}^n i[p^k|i][p^{k+1}\nmid i]} \\ &=& \prod_p \prod_k (k+1)^{\sum_{i=1}^n i[p^k|i] -\sum_{i=1}^ni[p^k|i]} \\ &=& \prod_p \prod_k (k+1)^{p^ks(\lfloor \frac n {p^k}\rfloor)-p^{k+1}s(\lfloor \frac n {p^{k+1}}\rfloor)} ,s(x)=\frac{x(x+1)}2 \\ S_2 &=& \prod_p 2^{\sum_{i=1}^n \mu(i)[p|i][p^{2}\nmid i]} \\ f(n) &=& \sum_{i=1}^n \mu(i)[p|i][p^{2}\nmid i] \\ &=& \sum_{i=1}^{\lfloor \frac np \rfloor} \mu(ip) [i \perp p] \\ &=& -\sum_{i=1}^{\lfloor \frac np \rfloor} \mu(i) [i \perp p] \\ &=& -\sum_{k\ge 1}\sum_{i=1}^{\lfloor \frac n{p^k} \rfloor} \mu(i) \end{eqnarray} \] 对于 \(S_1\)\(S_2\)\(p \le \sqrt n\) 的部分直接暴力计算,设剩下的部分是 \(S_1'\)\(S_2'\)\[ \begin{eqnarray} S_1'&=& 2^{\sum_{p>\sqrt n} p\cdot s(\lfloor \frac n p \rfloor)} \\ S_2'&=& 2^{-\sum_{p>\sqrt n}\sum_{i=1}^{\lfloor \frac np \rfloor}\mu(i)} \end{eqnarray} \] 还需要完成的部分有:

  • 杜教筛预处理 \(\sum_{i=1}^{\lfloor \frac n x \rfloor} \mu(i)\)
  • Min_25 筛 / 洲阁筛 中前半部分预处理 \(\sum_{p \le \lfloor \frac nx \rfloor} p\) 以及 \(\sum_{p \le \lfloor \frac n x \rfloor}\)

代码(在 51nod 上会 CE,但不给 CE 信息,没办法。保证不 TLE,测过大数据)

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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
#include <bits/stdc++.h>
using namespace std;
//using LL = long long;
#define FOR(i, x, y) for (decay<decltype(y)>::type i = (x), _##i = (y); i < _##i; ++i)
#define FORD(i, x, y) for (decay<decltype(x)>::type i = (x), _##i = (y); i > _##i; --i)
#ifdef zerol
#define dbg(args...) do { cout << "\033[32;1m" << #args << " -> "; err(args); } while (0)
#else
#define dbg(...)
#endif
void err() { cout << "\033[39;0m" << endl; }
template<template<typename...> class T, typename t, typename... Args>
void err(T<t> a, Args... args) { for (auto x: a) cout << x << ' '; err(args...); }
template<typename T, typename... Args>
void err(T a, Args... args) { cout << a << ' '; err(args...); }
// -----------------------------------------------------------------------------
using LL = __int128;
using ll = long long;
const LL MOD = 1000000000039, MOD1 = MOD - 1, INV2 = (MOD + 1) / 2;
namespace dujiao {
const int M = 2E7;
LL f[M] = {0, 1};
ll pr[2000000], p_sz, d;
void init() {
static bool vis[M];
FOR (i, 2, M) {
if (!vis[i]) { pr[p_sz++] = i; f[i] = -1; }
FOR (j, 0, p_sz) {
if ((d = pr[j] * i) >= M) break;
vis[d] = 1;
if (i % pr[j] == 0) {
f[d] = 0;
break;
} else f[d] = -f[i];
}
}
FOR (i, 2, M) f[i] += f[i - 1];
}
inline LL s_fg(LL n) { return 1; }
inline LL s_g(LL n) { return n; }

LL N, rd[M];
bool vis[M];
LL go(LL n) {
if (n < M) return f[n];
LL id = N / n;
if (vis[id]) return rd[id];
vis[id] = true;
LL& ret = rd[id] = s_fg(n);
for (LL l = 2, v, r; l <= n; l = r + 1) {
v = n / l; r = n / v;
ret -= (s_g(r) - s_g(l - 1)) * go(v);
}
return ret;
}
LL solve(LL n) {
N = n;
memset(vis, 0, sizeof vis);
return go(n);
}
}

namespace min25 {
const int M = 1E6 + 100;
LL B, N;

// g(x)
inline LL pg(LL x) { return 1; }
inline LL ph(LL x) { return x % MOD1; }
// Sum[g(i),{x,2,x}]
inline LL psg(LL x) { return x % MOD1 - 1; }
inline LL psh(LL x) {
return x * (x + 1) / 2 % MOD1 - 1;
}

ll pr[M], pc;
LL sg[M], sh[M];
void get_prime(ll n) {
static bool vis[M]; pc = 0;
FOR (i, 2, n + 1) {
if (!vis[i]) {
pr[pc++] = i;
sg[pc] = (sg[pc - 1] + pg(i)) % MOD1;
sh[pc] = (sh[pc - 1] + ph(i)) % MOD1;
}
FOR (j, 0, pc) {
if (pr[j] * i > n) break;
vis[pr[j] * i] = 1;
if (i % pr[j] == 0) break;
}
}
}

LL w[M];
LL id1[M], id2[M], h[M], g[M];
inline LL id(LL x) { return x <= B ? id1[x] : id2[N / x]; }

void init() { get_prime(M); }
void solve(LL _N) {
N = _N;
B = sqrt(N + 0.5);
int sz = 0;
for (LL l = 1, v, r; l <= N; l = r + 1) {
v = N / l; r = N / v;
w[sz] = v; g[sz] = psg(v); h[sz] = psh(v);
if (v <= B) id1[v] = sz; else id2[r] = sz;
sz++;
}
FOR (k, 0, pc) {
LL p = pr[k];
FOR (i, 0, sz) {
LL v = w[i]; if (p * p > v) break;
LL t = id(v / p);
g[i] = (g[i] - (g[t] - sg[k]) * pg(p)) % MOD1;
h[i] = (h[i] - (h[t] - sh[k]) * ph(p)) % MOD1;
}
}
}
}

inline LL s(LL x) {
return x * (x + 1) / 2 % MOD1;
}
LL bin(LL x, LL n) {
LL ret = 1;
for (n = (n % MOD1 + MOD1) % MOD1; n; n >>= 1, x = x * x % MOD)
if (n & 1) ret = ret * x % MOD;
return ret;
}

int main() {
#ifdef zerol
freopen("in", "r", stdin);
#endif
int T; cin >> T;
dujiao::init(); min25::init();
while (T--) {
long long _n; cin >> _n;
LL n = _n;
dujiao::solve(n);
min25::solve(n);

LL b = sqrt(n), bb = (n / (n / b) + 1);

LL ans = 1;
FOR (i, 0, dujiao::p_sz) {
LL p = dujiao::pr[i];
if (p >= bb) break;
for (LL e = 1, pp = p; pp <= n; ++e, pp *= p) {
ans = ans *
bin(e + 1, pp * s(n / pp) - pp * p * s(n / pp / p))
% MOD;
ans = ans * bin(INV2, dujiao::go(n / pp)) % MOD;
}
}
for (LL l = bb, v, r; l <= n; l = r + 1) {
using namespace min25;
v = n / l; r = n / v;
ans = ans * bin(2, s(v) * (h[id(r)] - h[id(l - 1)])) % MOD;
ans = ans * bin(2, -dujiao::go(v) * (g[id(r)] - g[id(l - 1)])) % MOD;
}

long long Ans = (ans % MOD + MOD) % MOD;
cout << Ans << endl;
}
}
w