f1,i f_{1,i} 的贡献是 f1,i∗an−i∗bn−1∗Cn−i2n−i−2 f_{1,i}*a^{n-i}*b^{n-1}*C_{2n-i-2}^{n-i}
fi,1 f_{i,1} 同理 注意 f1,1 f_{1,1} 没什么卵用
然后 c c 还有贡献

c∑i=2n∑j=2nan−i∗bn−j∗Cn−i2n−i−j=c∑i=0n−2∑j=0n−2ai∗bj∗Cii+j

c\sum_{i=2}^n\sum_{j=2}^n a^{n-i}*b^{n-j}*C_{2n-i-j}^{n-i}=c\sum_{i=0}^{n-2}\sum_{j=0}^{n-2} a^i*b^j*C_{i+j}^{i}
写成卷积的形式然后FFT就好了

正解似乎是递推

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;const int P=1e6+3;
const int M[]={998244353,1004535809,469762049};
const int G[]={3,3,3};
const ll _M=(ll)M[0]*M[1];inline ll Pow(ll a,int b,int p){ll ret=1;for (;b;b>>=1,a=a*a%p)if (b&1)ret=ret*a%p;return ret;
}
inline ll mul(ll a,ll b,ll p){a%=p; b%=p;return ((a*b-(ll)((ll)((long double)a/p*b+1e-3)*p))%p+p)%p;
}const int m1=M[0],m2=M[1],m3=M[2];
const int inv1=Pow(m1%m2,m2-2,m2),inv2=Pow(m2%m1,m1-2,m1),inv12=Pow(_M%m3,m3-2,m3);
inline int CRT(int a1,int a2,int a3){ll A=(mul((ll)a1*m2%_M,inv2,_M)+mul((ll)a2*m1%_M,inv1,_M))%_M;ll k=((ll)a3+m3-A%m3)*inv12%m3;return (k*(_M%P)+A)%P;
}const int N=600005;int n,A,B,C;
ll ans;#define read(x) scanf("%d",&(x))ll fac[N],inv[N];
inline void Pre(int n){fac[0]=1; for (int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;inv[1]=1; for (int i=2;i<=n;i++) inv[i]=(ll)(P-P/i)*inv[P%i]%P;inv[0]=1; for (int i=1;i<=n;i++) inv[i]=inv[i]*inv[i-1]%P;
}
inline ll _C(int n,int m){return fac[n]*inv[m]%P*inv[n-m]%P;
}int m;
int R[N];struct NTT{int P,G;int num,w[2][N];int R[N];void Pre(int _P,int _G,int m){num=m; P=_P; G=_G;int g=Pow(G,(P-1)/num,P);w[1][0]=1; for (int i=1;i<num;i++) w[1][i]=(ll)w[1][i-1]*g%P;w[0][0]=1; for (int i=1;i<num;i++) w[0][i]=w[1][num-i];int L=0; while (m>>=1) L++;for (int i=1;i<=num;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));}void FFT(int *a,int n,int r){for (int i=0;i<n;i++) if (i<R[i]) swap(a[i],a[R[i]]);for (int i=1;i<n;i<<=1)for (int j=0;j<n;j+=(i<<1))for (int k=0;k<i;k++){int x=a[j+k],y=(ll)a[j+i+k]*w[r][num/(i<<1)*k]%P;a[j+k]=(x+y)%P; a[j+i+k]=(x+P-y)%P;}if (!r) for (int i=0,inv=Pow(n,P-2,P);i<n;i++) a[i]=(ll)a[i]*inv%P;}
}ntt[3];int a[N],b[N],c[N],d[N];int _a[3][N];int main(){int x; freopen("t.in","r",stdin);freopen("t.out","w",stdout);read(n); read(A); read(B); read(C); Pre(2*n);for (int i=1;i<=n;i++)read(x),ans+=(i>1)*(ll)x*Pow(A,n-1,P)%P*Pow(B,n-i,P)%P*_C(2*n-i-2,n-i)%P;for (int i=1;i<=n;i++)read(x),ans+=(i>1)*(ll)x*Pow(A,n-i,P)%P*Pow(B,n-1,P)%P*_C(2*n-i-2,n-i)%P;ans%=P;int L=0; m=1;while (m<=(2*n-4)) m<<=1,L++;for (int i=1;i<m;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));for (int i=0;i<=n-2;i++)a[i]=Pow(A,i,P)*inv[i]%P,b[i]=Pow(B,i,P)*inv[i]%P;for (int i=0;i<3;i++){ntt[i].Pre(M[i],G[i],m);memcpy(c,a,sizeof(int)*(m+5)); memcpy(d,b,sizeof(int)*(m+5));ntt[i].FFT(c,m,1); ntt[i].FFT(d,m,1);for (int j=0;j<m;j++) c[j]=(ll)c[j]*d[j]%ntt[i].P;ntt[i].FFT(c,m,0);for (int j=0;j<m;j++) _a[i][j]=c[j];}ll ret=0;for (int i=0;i<=2*n-4;i++)ret+=fac[i]*CRT(_a[0][i],_a[1][i],_a[2][i])%P;ret=(ret%P)*C%P;ans=(ans+ret)%P;printf("%lld\n",ans);return 0;
}

[FFT || 递推] BZOJ 4451 [Cerc2015]Frightful Formula相关推荐

  1. bzoj 4451 : [Cerc2015]Frightful Formula FFT

    4451: [Cerc2015]Frightful Formula Time Limit: 10 Sec  Memory Limit: 64 MB Submit: 177  Solved: 57 [S ...

  2. BZOJ 4451: [Cerc2015]Frightful Formula

    Description 给你一个n*n矩阵的第一行和第一列,其余的数通过如下公式推出:  \(f_{i,j}=a\times f_{i,j-1}+b\times f_{i-1,j}+c\) 求\(f_ ...

  3. bzoj 4451: [Cerc2015]Frightful Formula 数学+排列组合

    题意 给你一个n*n矩阵的第一行和第一列,其余的数通过如下公式推出: F[i,j]=a*f[i,j-1]+b*f[i-1,j]+c 求f[n][n]%(10^6+3) 2<=n<=2000 ...

  4. 【bzoj 4451】[Cerc2015]Frightful Formula - 递推

    才没有在做cerc2015呢 看到好像不少人这题写fft卡得死死的啊,不如 O(n) O(n) 递推(雾) 首先可以观察出 (i,1) (i,1)这个格子为 x x时对(n,n)(n,n) a,b a ...

  5. BZOJ4451 [Cerc2015]Frightful Formula 多项式 FFT 递推 组合数学

    原文链接http://www.cnblogs.com/zhouzhendong/p/8820963.html 题目传送门 - BZOJ4451 题意 给你一个$n\times n$矩阵的第一行和第一列 ...

  6. BZOJ4451 : [Cerc2015]Frightful Formula

    $(i,1)$对答案的贡献为$l_iC(2n-i-2,n-i)a^{n-1}b^{n-i}$. $(1,i)$对答案的贡献为$t_iC(2n-i-2,n-i)*a^{n-i}b^{n-1}$. $(i ...

  7. CERC2015 Frightful Formula 神奇的模意义下分数

    上网查了下这道题的正解是FFT......然而机智的xushu用待定系数法A了.... f[i,j]+k=a(f[i][j-1]+k)+b(f[i-1,j]+k) 解得 k=c/(a+b-1) 于是令 ...

  8. HYSBZ(BZOJ) 4300 绝世好题(位运算,递推)

    HYSBZ(BZOJ) 4300 绝世好题(位运算,递推) Description 给定一个长度为n的数列ai,求ai的子序列bi的最长长度,满足bi&bi-1!=0(2<=i<= ...

  9. Bzoj 1046: [HAOI2007]上升序列 二分,递推

    1046: [HAOI2007]上升序列 Time Limit: 10 Sec  Memory Limit: 162 MB Submit: 3671  Solved: 1255 [Submit][St ...

最新文章

  1. php 中continue break exit return 的区别
  2. 2021年春季学期-信号与系统-第四次作业参考答案-第二小题
  3. maven整合jar包下载地址
  4. UltraISO软碟通U盘安装Centos7 的各种报错及解决方案
  5. 【每日一包0029】merge-descriptors
  6. 二叉树的前中后序遍历之迭代法(非统一风格迭代方式)
  7. 64位jdk连接32位的mysql_在64位客户端使用32位的ODBC配置
  8. CodeForces 1110H. Modest Substrings
  9. 编写程序在窗口中写出自己名字的拼音缩写_各类英文缩写:全称居然这么朴素?网友:最后一个我笑了...
  10. FTP 两种传输模式 Binary 和 ASCII 的区别
  11. 代码实现——MapReduce实现Hadoop序列化
  12. oracle创建索引01652,建立数据表快照导致ora-01652异常
  13. 基础知识:BT1120
  14. 11 月全国程序员平均工资出炉
  15. Xtts v4 xttdriver.pl xtt.properties
  16. 《第一行代码——Android》封面诞生记
  17. 网页里面的空格的代码怎么写
  18. 【Nova】nova-consoleauth学习
  19. 网站服务器采用CDN+专线 ,完美加速
  20. 王爽汇编语言 实验8

热门文章

  1. python抢票代码_Python自动化xpath实现自动抢票抢货
  2. 三星Gear VR体验:看 · 未来
  3. Storm zk目录结构
  4. ElasticSearch(六)【分词器】
  5. 什么是DFU模式和什么是DUF模式,有什么区别?ipad如何进入DFU模式和DUF模式?
  6. linux jq 遍历数组,jquery怎么遍历数组?
  7. Oracle数据库还原命令
  8. “剥皮”法求区域中心 (转)
  9. Linux 清空缓存命令
  10. linux chmod命令参数及用法详解