【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}inom{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}inom{k}{j}i^{j}R^{i})
((R - 1)f_{k}(N) = N^{k}R^{N + 1} + sum_{j = 0}^{k - 1}(-1)^{k - j}inom{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}inom{k}{j} f_{j}(N))

似乎。。。没法搞了

然而我们可以试试暴力展开一下

(f_{0}(N) = frac{R^{N + 1} - R}{R - 1} = R^{N + 1}cdot (frac{1}{R - 1}) - Rcdot (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}}) - Rcdot (-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}}) - Rcdot (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}inom{x}{i}inom{x - 1 - i}{k - i} F_{k}(i))
当(x = K + 1)时
有
(F_{k}(k + 1) = sum_{i = 0}^{k} (-1)^{k - i}inom{x}{i}F_{k}(i))
也就是
(sum_{i = 0}^{k + 1} (-1)^{k - i}inom{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('
')
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;
}