预处理伯努利数(多项式求逆)
2015-09-22 23:02
337 查看
伯努利数定义:
![](http://img.blog.csdn.net/20150922223144963?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
递推式:
![](http://img.blog.csdn.net/20150922223447728?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
将e^t展开得到:
![](http://img.blog.csdn.net/20150922224202009?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
然后就可以利用多项式逆元求出伯努利数了。。
多项式逆元就是给定一个n次多项式A(x),求另一个多项式B(x),满足
![](http://img.blog.csdn.net/20150922224613488?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
其中mod x^n表示多项式除以x^n剩余的多项式,这里可以看成做卷积之后去掉大于等于n次的项。
假设已经求得B(x)满足:
![](http://img.blog.csdn.net/20150922224850551?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
要求另一个B'(x)满足:
![](http://img.blog.csdn.net/20150922230336362?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
因为:
![](http://img.blog.csdn.net/20150923110311132?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
然后就可以递归的求出B(x)了。
详情见:
Inverse Element of Polynomial
和2015国家集训队论文:生成函数的运算与组合计数
PS:自己实现代码的时候还是坑点好多,暴力对拍了一下10000内的数据还是没问题的。。
1.做NTT的时候要做当前长度的两倍,因为B(x)^2*A(x)最高次会大于n。
2.做完NTT的多项式并不是伯努利数,还要乘上i!。
代码实现用求B
%1004535809
递推式:
将e^t展开得到:
然后就可以利用多项式逆元求出伯努利数了。。
多项式逆元就是给定一个n次多项式A(x),求另一个多项式B(x),满足
其中mod x^n表示多项式除以x^n剩余的多项式,这里可以看成做卷积之后去掉大于等于n次的项。
假设已经求得B(x)满足:
要求另一个B'(x)满足:
因为:
然后就可以递归的求出B(x)了。
详情见:
Inverse Element of Polynomial
和2015国家集训队论文:生成函数的运算与组合计数
PS:自己实现代码的时候还是坑点好多,暴力对拍了一下10000内的数据还是没问题的。。
1.做NTT的时候要做当前长度的两倍,因为B(x)^2*A(x)最高次会大于n。
2.做完NTT的多项式并不是伯努利数,还要乘上i!。
代码实现用求B
%1004535809
#include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <cmath> #include <queue> using namespace std; #define LL long long #define eps 1e-8 #define inf 0x3f3f3f3f #define MP make_pair #define N 300020 #define M 300020 #pragma comment(linker, "/STACK:1024000000,1024000000") #define Pi acos(-1.0) #define mod 1004535809 #define ls (i << 1) #define rs (ls | 1) #define md (ll + rr) >> 1 #define lson (ll, md, ls) #define rson (md + 1, rr, rs); //const LL P = 50000000001507329LL; //190734863287 * 2 ^ 18 + 1 常数巨大 const int P = 1004535809LL; //479 * 2 ^ 21 + 1 //const int P = 998244353; // 119 * 2 ^ 23 + 1 const int G = 3; int wn[25]; int n; LL mul(LL x, LL y) { return (x * y - (LL)(x / (long double)P * y + 1e-3) * P + P) % P; } int qpow(int x, int k, int p) { int ret = 1; while(k) { if(k & 1) ret = 1LL * ret * x % p; k >>= 1; x = 1LL * x * x % p; } return ret; } void getwn() { for(int i = 1; i <= 21; ++i) { int t = 1 << i; wn[i] = qpow(G, (P - 1) / t, P); } } void change(int *y, int len) { for(int i = 1, j = len / 2; i < len - 1; ++i) { if(i < j) swap(y[i], y[j]); int k = len / 2; while(j >= k) { j -= k; k /= 2; } j += k; } } void NTT(int *y, int len, int on) { change(y, len); int id = 0; for(int h = 2; h <= len; h <<= 1) { ++id; for(int j = 0; j < len; j += h) { int w = 1; for(int k = j; k < j + h / 2; ++k) { int u = y[k]; int t = 1LL * y[k+h/2] * w % P; y[k] = u + t; if(y[k] >= P) y[k] -= P; y[k+h/2] = u - t + P; if(y[k+h/2] >= P) y[k+h/2] -= P; w = 1LL * w * wn[id] % P; } } } if(on == -1) { for(int i = 1; i < len / 2; ++i) swap(y[i], y[len-i]); int inv = qpow(len, P - 2, P); for(int i = 0; i < len; ++i) y[i] = 1LL * y[i] * inv % P; } } int B ; int f , nf ; int a ; int tmp ; void get_inv(int A[], int A0[], int t) { if(t == 1) { A0[0] = qpow(A[0], P - 2, P); return; } get_inv(A, A0, t / 2); for(int i = 0; i < t; ++i) tmp[i] = A[i]; for(int i = t; i < 2 * t; ++i) tmp[i] = 0; for(int i = t / 2; i < 2 * t; ++i) A0[i] = 0; NTT(tmp, 2 * t, 1); NTT(A0, 2 * t, 1); for(int i = 0; i < 2 * t; ++i) { tmp[i] = (2 - 1LL * tmp[i] * A0[i] % P) % P; if(tmp[i] < 0) tmp[i] += P; A0[i] = 1LL * A0[i] * tmp[i] % P; } NTT(A0, 2 * t, -1); } void init() { f[0] = 1; for(int i = 1; i < N; ++i) f[i] = 1LL * f[i-1] * i % P; nf[N-1] = qpow(f[N-1], P - 2, P); for(int i = N - 2; i >= 0; --i) { nf[i] = 1LL * nf[i+1] * (i + 1) % P; } for(int i = 0; i < N - 1; ++i) a[i] = nf[i+1]; int len = 1 << 17; get_inv(a, B, len); for(int i = 0; i < len; ++i) B[i] = 1LL * B[i] * f[i] % P; } int main() { getwn(); init(); int n; int cas; scanf("%d", &cas); while(cas--) { scanf("%d", &n); printf("%d\n", B ); } return 0; }
相关文章推荐
- c语言初级小程序
- Android视图SurfaceView的实现原理分析
- mysql--预排序遍历树(数据库级层存储的思考)
- iOS讲解迷惑深入浅出之GCD利用屏障模拟玩家进入游戏读取数据
- 常用正则表达式
- 第九章 IP选路
- JDBC链接SQLServer详解
- 《C++primer(第五版)》学习之路-第七章:类
- JavaScript学习总结(六)数据类型和JSON格式
- React Native通信机制详解
- 数据存储和xml
- 设计模式之单例模式
- swing展开所有的结点
- Storm系列(十)聚流示例
- [转]原生js的String类扩展
- spinner适配器
- hdoj 2444 The Accomodation of Students 【黑白染色判二分图 + 最大匹配】
- 黑马程序员——Java基础---多态、内部类、异常、包
- linux网口驱动实现(待续)
- PHP 性能分析第二篇: Xhgui In-Depth