【51nod】1822 序列求和 V5
題解
我是zz吧
nonprime[i * prime[j]] = 0
= =
還以為是要卡常，卡了半天就是過不掉
我們來說這道題……
首先，我們考慮一個\(K^2\)做法
\(f_{k}(N) = \sum_{i = 1}^{N} i^{k}R^{i}\)
\((R - 1)f_{k}(N) = \sum_{i = 1}^{N}i^{k}R^{i + 1} - \sum_{i = 1}^{N} i^{k}R^{i}\)
\((R - 1)f_{k}(N) = N^{k}R^{N + 1} + \sum_{i = 1}^{N} [(i - 1)^{k} - i^{k}]R^{i}\)
\((R - 1)f_{k}(N) = N^{k}R^{N + 1} + \sum_{i = 1}^{N} [\sum_{j = 0}^{k}(-1)^{k - j}\binom{k}{j}i^{j} - i^{k}]R^{i}\)
\((R - 1)f_{k}(N) = N^{k}R^{N + 1} + \sum_{i = 1}^{N} \sum_{j = 0}^{k - 1}(-1)^{k - j}\binom{k}{j}i^{j}R^{i}\)
\((R - 1)f_{k}(N) = N^{k}R^{N + 1} + \sum_{j = 0}^{k - 1}(-1)^{k - j}\binom{k}{j}\sum_{i = 1}^{N} i^{j}R^{i}\)
\((R - 1)f_{k}(N) = N^{k}R^{N + 1} + \sum_{j = 0}^{k - 1}(-1)^{k - j}\binom{k}{j} f_{j}(N)\)
似乎。。。沒法搞了
然而我們可以試試暴力展開一下
\(f_{0}(N) = \frac{R^{N + 1} - R}{R - 1} = R^{N + 1}\cdot (\frac{1}{R - 1}) - R\cdot (\frac{1}{R - 1})\)
$f_{1}(N) = \frac{N \cdot R^{N + 1}}{R - 1} - \frac{f_{0}(N)}{R - 1} = R^{N + 1}(\frac{N}{R - 1} - \frac{1}{(R - 1)^{2}}) - R\cdot (-\frac{1}{(R - 1)^2}) $
$f_{2}(N) = \frac{N^{2}R^{N + 1}}{R - 1} - \frac{2f_{1}(N)}{R - 1} + \frac{f_{0}(N)}{R - 1}= R^{N + 1}(\frac{N^{2}}{R - 1} - \frac{2N}{(R - 1)^{2}} + \frac{2}{(R - 1)^{3}} + \frac{1}{(R - 1)^{2}}) - R\cdot (\frac{1}{(R - 1)^2} + \frac{2}{(R - 1)^{3}}) $
我們可以猜測，並且歸納證明這個結論
\(f_{k}(N) = R^{N} F_{k}(N) - F_{k}(0)\)
其中\(F_{k}(N)\)是一個k次多項式
證明了多項式，那麽就考慮插值
如何插值，我們發現按照定義有$F_{k}(N + 1)= \frac{F_{k}(N)}{R} + (N + 1)^{k} $
我們設\(F_{k}(0) = x\)，然後用x表示剩下\(K + 1\)個多項式
同時，我們嘗試用插值法表示出\(F_{k}(K + 1)\)
式子經過整理後也就是
\(F_{k}(x) = \sum_{i = 0}^{k} (-1)^{k - i}\binom{x}{i}\binom{x - 1 - i}{k - i} F_{k}(i)\)
當\(x = K + 1\)時
有
\(F_{k}(k + 1) = \sum_{i = 0}^{k} (-1)^{k - i}\binom{x}{i}F_{k}(i)\)
也就是
\(\sum_{i = 0}^{k + 1} (-1)^{k - i}\binom{x}{i}F_{k}(i) = 0\)
代入即可解出\(F_{k}(0)\)從而求出F的每個點值
代碼
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <vector>
#include <set>
//#define ivorysi
#define eps 1e-8
#define mo 974711
#define pb push_back
#define mp make_pair
#define pii pair<int,int>
#define fi first
#define se second
#define MAXN 200005
#define space putchar(‘ ‘)
#define enter putchar(‘\n‘)
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef unsigned long long u64;
typedef double db;
template<class T>
void read(T &res) {
res = 0;char c = getchar();T f = 1;
while(c < ‘0‘ || c > ‘9‘) {
if(c == ‘-‘) f = -1;
c = getchar();
}
while(c >= ‘0‘ && c <= ‘9‘) {
res = res * 10 + c - ‘0‘;
c = getchar();
}
res *= f;
}
template<class T>
void out(T x) {
if(x < 0) putchar(‘-‘);
if(x >= 10) {
out(x / 10);
}
putchar(‘0‘ + x % 10);
}
const int64 MOD = 985661441;
int64 fpow(int64 x,int64 c) {
int64 res = 1,t = x;
while(c) {
if(c & 1) res = res * t % MOD;
t = t * t % MOD;
c >>= 1;
}
return res;
}
int T,K;
int prime[MAXN],tot;
int64 N,R,F[MAXN];
int64 a[MAXN],b[MAXN];
int64 MK[MAXN];
int64 fac[MAXN],invfac[MAXN],Le[MAXN],Ri[MAXN];
bool nonprime[MAXN];
int64 C(int n,int m) {
if(n < m) return 0;
return fac[n] * invfac[n - m] % MOD * invfac[m] % MOD;
}
void Solve() {
read(K);read(N);read(R);
R %= MOD;
if(R == 0) {puts("0");return;}
MK[1] = 1;
memset(nonprime,0,sizeof(nonprime));
tot = 0;
for(int i = 2 ; i <= K + 2 ; ++i) {
if(!nonprime[i]) {
prime[++tot] = i;
MK[i] = fpow(i,K);
}
for(int j = 1 ; j <= tot ; ++j) {
if(prime[j] > (K + 2) / i) break;
nonprime[i * prime[j]] = 1;
MK[i * prime[j]] = MK[i] * MK[prime[j]] % MOD;
if(i % prime[j] == 0) break;
}
}
if(R == 1) {
F[0] = 0;
for(int i = 1 ; i <= K + 2 ; ++i) F[i] = (F[i - 1] + MK[i]) % MOD;
if(N <= K + 2) {out(F[N]);enter;return;}
int64 t = 1,ans = 0;
Le[0] = 1;
N %= MOD;
for(int i = 1 ; i <= K + 2 ; ++i) {
Le[i] = Le[i - 1] * (N + MOD - i) % MOD;
}
Ri[K + 3] = 1;
for(int i = K + 2 ; i >= 1 ; --i) {
Ri[i] = Ri[i + 1] * (N + MOD - i) % MOD;
}
for(int i = K + 2 ; i >= 1 ; --i) {
ans += t * invfac[i - 1] % MOD * invfac[K + 2 - i] % MOD * Le[i - 1] % MOD * Ri[i + 1] % MOD * F[i] % MOD;
t = t * (MOD - 1) % MOD;
}
ans %= MOD;
out(ans);enter;return;
}
a[0] = 1,b[0] = 0;
int64 InvR = fpow(R,MOD - 2);
for(int i = 1 ; i <= K + 1; ++i) {
a[i] = 1LL * a[i - 1] * InvR % MOD;
b[i] = (1LL * b[i - 1] * InvR + MK[i]) % MOD;
}
int64 suma = 0,sumb = 0,t = 1;
if(K & 1) t = MOD - 1;
for(int i = 0 ; i <= K + 1; ++i) {
int64 h = t * C(K + 1,i) % MOD;
suma = (suma + a[i] * h) % MOD;
sumb = (sumb + b[i] * h) % MOD;
t = t * (MOD - 1) % MOD;
}
F[0] = (MOD - sumb) * fpow(suma,MOD - 2) % MOD;
for(int i = 1 ; i <= K + 1; ++i) F[i] = (a[i] * F[0] + b[i]) % MOD;
int64 ans = 0;
t = 1;
if(N <= K) {
ans = (fpow(R,N) * F[N] - F[0] + MOD) % MOD;
out(ans);enter;
return ;
}
int64 T = N % (MOD - 1);
N %= MOD;
Le[0] = N;
for(int i = 1 ; i <= K ; ++i) {
Le[i] = Le[i - 1] * (N + MOD - i) % MOD;
}
Ri[K + 1] = 1;
for(int i = K ; i >= 0 ; --i) {
Ri[i] = Ri[i + 1] * (N + MOD - i) % MOD;
}
if(K & 1) t = MOD - 1;
for(int i = 0 ; i <= K; ++i) {
int64 h;
if(i == 0) h = t * invfac[i] % MOD * invfac[K - i] % MOD * Ri[i + 1] % MOD * F[i] % MOD;
else h = t * invfac[i] % MOD * invfac[K - i] % MOD * Le[i - 1] % MOD * Ri[i + 1] % MOD * F[i] % MOD;
t = t * (MOD - 1) % MOD;
ans += h;
}
ans %= MOD;
ans = (ans * fpow(R,T) - F[0] + MOD) % MOD;
out(ans);enter;
}
int main() {
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
fac[0] = 1;
for(int i = 1 ; i <= 200001 ; ++i) fac[i] = fac[i - 1] * i % MOD;
invfac[200001] = fpow(fac[200001],MOD - 2);
for(int i = 200000 ; i >= 0 ; --i) invfac[i] = invfac[i + 1] * (i + 1) % MOD;
read(T);
while(T--) {
Solve();
}
return 0;
}【51nod】1822 序列求和 V5