您的位置:首页 > 其它

预处理伯努利数(多项式求逆)

2015-09-22 23:02 337 查看
伯努利数定义:

递推式:


将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;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: