【题解】[湖北省队互测2014week1 B] 一个人的数论

给定 $n$ 的质因数分解 $n= \prod _{i=1}^kp^{a_{i}}$,求所有小于 $n$ 且与 $n$ 互质的正整数的 $m$ 次方之和模 $10^9 +7$ 的值。其中 $p_i$ 为质数且不超过 $10^9$ 。

$a_i ≤ 10^9,k ≤ 1000,m ≤ 100$。

写出式子来

然后发现这东西好像不是容易提前预处理。根据伯努利数的推论,可以知道 $\sum i^m$ 这东西是一个关于上界 $n$ 的 $m+1$ 次多项式。发现 $m$ 并不大于是可以 $O(m^2)$ 插出来。 若记这个多项式是 $f$ ,那么原式就等价于

稍微化一下就是

发现对于后面的 $\sum $ 只需要求出每个 $p_i^{a_i}$ 处的值,然后每次 $O(k)$ 暴力合并,似乎也没啥问题。考虑对于每个 $p_i^{a_{i}}$ 怎么求。这个地方大概需要涨个经验,就是 $\mu(d)$ 这东西只有当 $d=1$ 和 $d=p_i$ 的时候才有值(无平方因子,剩下的因子都是 $p_i$ 的某个次数 $>2$ 的幂),所以每个只需要算两次,是 $O(1)$ 的。于是最后的复杂度 $O(m(m+k))$ 。

好神啊好神啊。

哈哈哈哈哈草这个题居然允许用高斯消元来代替插值。毕竟 $m$ 只有 $100$ 。

然后关于拉插,感觉很gg,因为一共有三个版本的拉插,求值和求系数,求值的又分为 $x$ 连续/不连续的。于是今天就把三个版本都写了一遍……

不过有一说一,求系数的拉插感觉就是在模拟。然后我还对着多项式除法(整除以一个形如 $(x+t)$ 的多项式)摸了好久,一直感觉很迷乱,后来用了用大除法,发现就是在模拟大除法的过程罢了。

然后一开始我还憨憨的写了个线性筛,然后发现 $\mu(d)$ 的 $d$ 是 $1e9$ 的,就懵圈了,然后又发现 $\mu(d)$ 只会是 $1/-1$ 并且很好判断,于是感觉自己是个sb。

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
#include <bits/stdc++.h>

using namespace std ;

typedef long long ll ;

const int N = 1000010 ;
const int P = 1000000007 ;

int n, k ;

void debug(int *tp, int minn, int maxn, char c){
for (int i = minn ; i <= maxn ; ++ i)
cout << tp[i] << " " ; cout << c ;
}

namespace Interpolation{
int ans ;
int now ;
int x[N] ;
int y[N] ;
int fac[N] ;
int inv[N] ;
int pres[N] ;
int sufs[N] ;
int expow(int a, int b){
int res = 1 ;
a = (a % P + P) % P ;
while (b){
if (b & 1)
res = (ll)res * a % P ;
a = (ll)a * a % P ; b >>= 1 ;
}
return res ;
}
int fz[N] ;
int fm[N] ;
int tmp[N] ;
int res[N] ;
void add(int &a, int b){
a += b ;
if (a > P) a -= P ;
}
void dec(int &a, int b){
(a -= b) %= P ;
if (a < 0) a += P ;
}
void fmul(int *t, int deg, int opt){
for (int i = deg + 1 ; i >= 1 ; -- i)
tmp[i] = t[i], t[i] = t[i - 1] ;
for (int i = 1 ; i <= deg + 1 ; ++ i)
add(t[i], (ll)opt * tmp[i] % P) ;
}
void fdiv(int *t, int *ret, int deg, int opt){
for (int i = 1 ; i <= deg ; ++ i) tmp[i] = t[i] ;
for (int i = deg - 1 ; i >= 1 ; -- i)
ret[i] = tmp[i + 1], dec(tmp[i], (ll)tmp[i + 1] * opt % P) ;
}
void get_xs(int n){
fz[1] = 1 ;
for (int i = 1 ; i <= n ; ++ i)
fmul(fz, i, (-x[i] + P) % P) ;
//debug(fz, 1, n + 1, '\n') ;
for (int i = 1 ; i <= n ; ++ i){
int fenmu = 1 ;
for (int j = 1 ; j <= n ; ++ j)
if (i != j) fenmu = (ll)fenmu * (x[i] - x[j] + P) % P ;
fdiv(fz, fm, n + 1, -x[i]) ;
fenmu = (ll)y[i] * expow(fenmu, P - 2) % P ;
// cout << fenmu << endl ;
for (int j = 1 ; j <= n ; ++ j)
add(res[j], (ll)fenmu * fm[j] % P) ;
// debug(res, 1, n, '\n') ;
}
}
void pre_do(int U){
fac[0] = 1 ;
for (int i = 1 ; i <= U ; ++ i)
fac[i] = (ll)fac[i - 1] * i % P ;
inv[U] = expow(fac[U], P - 2) ;
for (int i = U ; i >= 1 ; -- i)
inv[i - 1] = (ll)inv[i] * i % P ;
//debug(fac, 1, U, '\n') ;
//debug(inv, 1, U, '\n') ;
}
int evenmark(int x){
return (x & 1) ? -1 : 1 ;
}
void get_dnx(int n, int k, bool m){
if (m){
pre_do(n) ;
pres[0] = 1, sufs[n + 1] = 1 ;
for (int i = 1 ; i <= n ; ++ i)
pres[i] = ((ll)pres[i - 1] * (k - x[i]) % P + P) % P ;
for (int i = n ; i >= 1 ; -- i)
sufs[i] = ((ll)sufs[i + 1] * (k - x[i]) % P + P) % P ;
for (int i = 1 ; i <= n ; ++ i){
now = (ll)pres[i - 1] * sufs[i + 1] % P ;
now = (ll)now * expow((ll)fac[i - 1] * fac[n - i] % P, P - 2) % P ;
now = (ll)evenmark(n - i) * y[i] % P * now % P ; ans = (ll)(ans + now) % P ;
}
}
else {
int inow = 1 ;
pres[0] = 1, sufs[n + 1] = 1 ;
for (int i = 1 ; i <= n ; ++ i)
pres[i] = ((ll)pres[i - 1] * (k - x[i]) % P + P) % P ;
for (int i = n ; i >= 1 ; -- i)
sufs[i] = ((ll)sufs[i + 1] * (k - x[i]) % P + P) % P ;
for (int i = 1 ; i <= n ; ++ i){
inow = 1, now = (ll)pres[i - 1] * sufs[i + 1] % P ;
for (int j = 1 ; j <= n ; ++ j)
if (i != j) inow = (ll)inow * ((x[i] - x[j]) % P + P) % P ;
ans = (ans + (ll)now * expow(inow, P - 2) % P * y[i] % P) % P ;
}
}
}
}/*
namespace Linear_sieve{
int cnt ;
int pr[N] ;
int mu[N] ;
int vis[N] ;
void sieve(int U){
mu[1] = 1 ;
for (int i = 2 ; i <= U ; ++ i){
if (!vis[i]) mu[i] = -1, pr[++ cnt] = i ;
for (int j = 1 ; j <= cnt ; ++ j){
if (i * pr[j] > U) break ;
vis[i * pr[j]] = 1 ;
if (i % pr[j] == 0) break ;
mu[i * pr[j]] = -mu[i] ;
}
}
}
}*/
int d, w, num ;
int base[N], cs[N] ;
int main(){
cin >> d >> w ; n = d + 2 ;
// using namespace Linear_sieve ;
using namespace Interpolation ;
for (int i = 1 ; i <= n ; ++ i)
x[i] = i, y[i] = (y[i - 1] + expow(x[i], d)) % P ;
for (int i = 1 ; i <= w ; ++ i)
cin >> base[i] >> cs[i] ; get_xs(n) ;
// debug(y, 1, n + 1, '\n') ;
// debug(res, 1, n + 1, '\n') ;
for (int i = 1 ; i <= n + 1 ; ++ i){
now = 0, num = 1 ;
for (int j = 1 ; j <= w ; ++ j){
now = expow(expow(base[j], cs[j] - 1), i - 1) ;
now = ((-1ll * now * expow(base[j], d) % P) + P) % P ;
now = (now + expow(expow(base[j], cs[j]), i - 1)) % P ;
num = 1ll * num * now % P ;
}
ans = (ans + (ll)num * res[i] % P) % P ;
}
cout << ans << '\n' ; return 0 ;
}