从FJWC看到的一个题, 比较容易想到的做法是分析多项式的次数, 并使用LAG插值法进行求解。

但是有比较优美但是颇为繁琐的STD做法? 用的是伯努利数那套理论。

传送门

题目大意

求:

\[\sum_{i=0}^{n} \sum_{j=1}^{a+i \times d} \sum_{l=1}^{j}l^k\]

解题报告

  • 算法一
    • \(f[n]=S(n,k)=\sum_{i=1}^{n} i^k\), 则\(f\)为对一个\(k\)次多项式求前缀和,所以\(f\)\(k+1\)次多项式;
    • 同理\(g[n]=\sum_{i=1}^{n} f[i]\),则\(g\)为对\(f\)求前缀和,\(g\)\(k+2\)次多项式。
    • 继续\(h[x]=g[a+x \times d]\), 则\(h[x]\)次数界与\(g\)相同, 为\(k+2\)次多项式。
    • 最后,\(ans[n]=\sum_{i=0}^{n}h[x]\)\(ans\)为对\(h\)求前缀和,所以\(ans\)\(k+4\)次多项式。
    • 因为\(g\)\(k+2\)次多项式,\(O(k \log{k})\) 预处理处\(g\)的前\(k+3\)项后,可以使用拉格朗日插值\(O(k)\)插出一个\(h[x]\), \(O(k^2)\)插出\(k+4\)\(h\)并求前缀和得到\(k+4\)\(ans[x]\),再使用\(k+4\)\(ans[x]\)进行插值得到\(ans[n]\)
  • 算法二
    • 题目中给出伯努利数,主要目的是提示多项式的次数界。

    • 但实际上可以使用伯努利数求解。

    • 承接算法一中的定义,\(g(x)=\sum_{i=0}^{n+2} g_i x^i\)

    • \[ \begin{aligned} ans[n] &=\sum_{i=0}^{n} g(a+i*d) \\ &= \sum_{i=0}^{n} \sum_{t=0}^{n+2} g_t (a+i*d)^t \\ &= \sum_{i=0}^{n} \sum_{t=0}^{n+2} g_t \sum_{k=0}^{t} \binom{t}{k}a^{t-k}(i \times d)^{k} \\ &= \sum_{t=0}^{n+2} g_t \sum_{k=0}^{t} \binom{t}{k}a^{t-k}d^k \sum_{i=0}^n i^k \end{aligned} \]

    • 然后发现, 如果能求出\(g_t\), 并能预处理\(\sum_{i=0}^{n}i^k\), 就可以\(k^2\)的搞了。

    • 翻一下伯努利数和自然数幂和

    • 摘选重要的公式:

    • \[\sum_{k=0}^{n} \binom{n+1}{k} B_k=0\]

    • 利用这个公式可以\(O(k^2)\)预处理伯努利数\(B\)

    • \[\sum_{i=1}^{x}i^k=\frac{1}{k+1}\sum_{i=1}^{k+1}\binom{k+1}{i}B_{k+1-i}(x+1)^i\]

    • 上式中,\((x+1)^i\)不够优美,因为多项式中需要\(x^i\)的形式, 所以把上式重写:

    • \[\sum_{i=0}^{x-1}i^k=\frac{1}{k+1}\sum_{i=1}^{k+1}\binom{k+1}{i}B_{k+1-i}x^i \quad (1)\]

    • 对g化式子: \[ \begin{aligned} g(x)&=\sum_{i=1}^{x} \sum_{j=1}^{i} j^k \\ &=\sum_{i=1}^{x} (x-i+1)*i^k\\ &=x*\sum_{i=0}^{x-1}i^k-\sum_{i=0}^{x-1}i^{k+1}+\sum_{i=0}^{x-1}i^k+x^k \end{aligned} \ (2) \]

    • 利用这个式子\((1)(2)\),\(O(k)\)预处理\(g_t\)

    • 利用式子\((1)\), \(\sum_{i=0}^{n}i^k\)单次\(O(k)\), 所以可以\(O(k^2)\)预处理\(S(n,k)\).
    • 这样就完结撒花喽~~

代码

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
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=int(a),nn=int(b);i<=nn;++i)
#define vep(i,a,b) for(int i=int(a),nn=int(b);i>=nn;--i)
#define xep(i,b) for(int i=0,nn=int(b);i<nn;++i)
typedef long long ll;
const int p=1234567891;
const int N=3015;
ll f[N],g[N],fc[N],iv[N],pe[N],su[N],k,a,n,d;
inline ll fast(ll x,int k) {ll as=1; for(;k;k>>=1,x=x*x%p) if(k&1) as=as*x%p;return as;}
inline ll mo(ll x, ll y) { x+=y; if (x>=p) x-=p; if (x<0) x+=p; return x;}
inline int P(ll x,ll y) { if ((x-y)&1) return p-1; else return 1;}
inline int lag(ll f[],ll n, int k) {
pe[0]=1; rep(i,1,k) pe[i]=pe[i-1]*(n-i+p)%p;
su[0]=1; rep(i,1,k) su[i]=su[i-1]*(n-k+i-1+p)%p;
ll as=0;
rep(i,1,k) {
ll up=pe[i-1]*su[k-i]%p*f[i]%p;
ll dn=iv[i-1]*iv[k-i]%p;
as=mo(as, up*dn%p*P(k,i)%p);
}
return as;
}
int main() {
fc[0]=1; rep(i,1,3010) fc[i]=fc[i-1]*i%p;
iv[3010]=fast(fc[3010],p-2); vep(i,3009,0) iv[i]=iv[i+1]*(i+1)%p;
int test; scanf("%d",&test);
while (test--) {
scanf("%d%d%d%d",&k,&a,&n,&d);
rep(i,0,k+3) g[i]=fast(i,k);
rep(i,1,k+3) g[i]=mo(g[i-1],g[i]);
rep(i,1,k+3) g[i]=mo(g[i-1],g[i]);
f[0]=lag(g,a,k+3);
rep(i,1,k+5) f[i]=mo(f[i-1],lag(g,(a+d*i)%p,k+3));
printf("%d\n",lag(f,n,k+5));
}
return 0;
}

留言