upd: typo了一个地方，感谢NaCly_Fish的指正。
____
首先写出式子来
然后发现这东西好像不是容易提前预处理。根据伯努利数的推论，可以知道 $\sum i^m$ 这东西是一个关于上界 $n$ 的 $m+1$ 次多项式。发现 $m$ 并不大于是可以 $O(m^2)$ 插出来。 若记这个多项式是 $f$ ，那么原式就等价于
$$\sum_{d|n}d^m\mu(d)\sum_{i=1}^{m+1} f_i\left(\lfloor\frac{n}{d}\rfloor\right) ^ i$$
稍微化一下就是
$$=\sum_{i=0}^{m+1} f_{i} \sum_{d | n} d^{m} \mu(d)\left(\frac{n}{d}\right)^{i}$$
发现对于后面的 $\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$ 。
```cpp
#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 ;
}
```