我非常智障的化了一部式子然后直接FFT.

但是因为模数是1000003, 这个模数有点大, 用double的精度显然不是很靠谱。

所以有需要使用什么任意模数FFT.

然后我折腾半天好不容易A了, 发现这个题可以递推做。

传送门

题目大意

\(f[i,j]=a*f[i,j-1]+b*f[i-1,j]+c\)

给出\(f[1,i]\)\(f[i,1]\) , 求\(f[n,n] (\mod 1000003)\)

解题报告

贡献分3部分:

  1. \((i, 1)\) 的贡献为\(f[i,1]C(2n-i-2,n-i)a^{n-1}b^{n-i}\) ;

  2. \((1,i)\)的贡献为\(f[1,i]C(2n-i-2,n-i)*a^{n-i}b^{n-1}\);

  3. \((i, j)\)\(c\)的贡献为

\[ \begin{eqnarray*} &&c\sum_{i=2}^n\sum_{j=2}^nC(2n-i-j,n-i)a^{n-j}b^{n-i}\\ &=&c\sum_{i=2}^n\sum_{j=2}^n(2n-i-j)!\times\frac{a^{n-j}}{(n-j)!}\times\frac{b^{n-i}}{(n-i)!} \end{eqnarray*} \]

然后直接FFT? 把系数写成\(aM+b\)的形式,然后做好几次\(FFT\)合并起来就不会炸精度了。

递推做法在这里

代码

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
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <cmath>
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;
typedef long double ff;
const int N=200100;
const int p=1000003;
const int M=1000;
const ff pi=acos(-1);
int n,_n,a,b,c,fc[N+N],rv[N+N],h[N],l[N],as,ma[N],mb[N],rr[N*4];
int A[N*3],B[N*3],C[N*3],a0[N*3],b0[N*3],a1[N*3],b1[N*3];
inline void in(int &x) {
char c=getchar(); int f=1;
for (;c<'0'||c>'9';c=getchar()) f=(c=='-'?-1:f);
for (x=0;c>='0'&&c<='9';c=getchar()) x=x*10+c-48;
x*=f;
}
inline int fast(int x,int k,int md=p) {
int as=1; for (; k; k>>=1, x=(ll)x*x%md)
if (k&1) as=(ll) as*x%md; return as;
}
inline int mo(int x,int y) { x+=y; if(x>=p) x-=p; if (x<0) x+=p; return x;}
struct cmx {
ff r, i; cmx(ff r=0, ff i=0)
: r(r), i(i) {}
cmx operator + (const cmx &b) const { return cmx(r+b.r, i+b.i); }
cmx operator - (const cmx &b) const { return cmx(r-b.r, i-b.i); }
cmx operator * (const cmx &b) const { return cmx(r*b.r-i*b.i,r*b.i+i*b.r); }
} aa[N*3],bb[N*3];
void fft(cmx *a, int f) {
xep(i,_n) if (rr[i]>i) swap(a[i],a[rr[i]]);
for (int m=1; m<_n; m<<=1) {
cmx wn=cmx(cos(pi/m), f*sin(pi/m));
for (int i=0; i<_n; i+=m<<1) {
cmx w=cmx(1,0);
for (int j=0; j<m; ++j) {
cmx x=a[i+j], y=w*a[i+j+m];
a[i+j]=x+y, a[i+j+m]=x-y;
w=w*wn;
}
}
}
if (f==-1) xep(i,_n) a[i].r/=_n;
}
inline void mul(int *a,int *b,int *c) {
memset(aa,0,sizeof(aa)); memset(bb, 0, sizeof(bb));
xep(i,_n) aa[i]=cmx(a[i],0), bb[i]=cmx(b[i],0);
fft(aa,1), fft(bb,1); xep(i,_n) aa[i]=aa[i]*bb[i];
fft(aa,-1); xep(i,_n) c[i]=(ll)(aa[i].r+0.5)%p;
}
inline void FFT_casual(int *a,int *b,int *c) {
xep(i,_n) a0[i]=a[i]/M, b0[i]=b[i]/M; int i;
for (mul(a0,b0,a0), i=0; i<_n; ++i) {
c[i]=1ll*a0[i]*M%p*M%p;
a1[i]=a[i]%M, b1[i]=b[i]%M;
}
for (mul(a1,b1,a1), i=0; i<_n; ++i) {
c[i]=(a1[i]+c[i])%p, a0[i]=(a0[i]+a1[i])%p;
a1[i]=a[i]/M+a[i]%M, b1[i]=b[i]/M+b[i]%M;
}
for (mul(a1,b1,a1), i=0; i<_n;++i)
c[i]=(1ll*M*(a1[i]-a0[i]+p)%p+c[i])%p;
}
int main() {
in(n),in(a),in(b),in(c);
rep(i,1,n) in(l[i]); rep(i,1,n) in(h[i]);
fc[0]=1; rep(i,1,n+n) fc[i]=(ll)fc[i-1]*i%p;
rv[n+n]=fast(fc[n+n],p-2); vep(i,n+n-1,0) rv[i]=(ll)rv[i+1]*(i+1)%p;
ma[0]=1; rep(i,1,n) ma[i]=(ll)ma[i-1]*a%p;
mb[0]=1; rep(i,1,n) mb[i]=(ll)mb[i-1]*b%p;
rep(i,2,n) as=mo(as,(ll)h[i]*ma[n-i]%p*mb[n-1]%p*fc[n-i+n-2]%p*rv[n-2]%p*rv[n-i]%p);
rep(i,2,n) as=mo(as,(ll)l[i]*ma[n-1]%p*mb[n-i]%p*fc[n-i+n-2]%p*rv[n-2]%p*rv[n-i]%p);
for (_n=1; _n<=n; _n<<=1) ; _n<<=1;
for (int i=0, j=0; i<_n; ++i) {
rr[i]=j; for (int k=_n>>1; (j^=k)<k;k>>=1);
}
xep(i,n-1) {
A[i]=(ll)ma[i]*rv[i]%p;
B[i]=(ll)mb[i]*rv[i]%p;
}
FFT_casual(A,B,C);
xep(i, n+n-3) as=mo(as, (ll)C[i]%p*c%p*fc[i]%p);
return printf("%d\n", as), 0;
}

留言