二分和三分
二分
- 单调性
(1) 单调性：
- 左边成立，右边不成立。
- 在成立者之中，找到最靠右的位置。其实就是找到分界点。
int l=1,r=n,res=-1;
while(l<=r){
int m=(l+r)>>1;
if(check(m)) res=m,l=m+1;
else r=m-1;
}
(2) 单调性：左边不成立，右边成立。
- 在成立者之中，找到最靠左的位置。
- 往外一格就是upper_bound，但这个不是内置lower_bound的位置（内置lower_bound指向第一个成立者）！
int l=1,r=n,res=-1;
while(l<=r){
int m=(l+r)>>1;
if(check(m)) res=m,r=m-1;
l=m+1;
}
(3) 单调性：中间某位置成立。左边过小，右边过大。
- 找中间某位置。
int l=1,r=n,res=-1;
while(l<=r){
int m=(l+r)>>1;
if(a[m]<5) l=m+1;
else if(a[m]>5) r=m-1;
else {res=m;break;}
}
(4) 数组从0开始，单调性：左边成立，右边不成立。
- 找最靠右的成立者。
int l=-1,r=n;
while(r-l>1){
int m=(l+r)>>1;
if(check(m)) l=m;
else r=m;
}
(5) 数组从0开始，单调性：左边过小，右边过大。中间成立。
- 找中间位置
int l=-1,r=n,res=-1;
while(r-l>1){
int m=(l+r)>>1;
if(a[m]<5) l=m;
else if(a[m]>5) r=m;
else {res=m;break;}
}
多项式
前置知识
数论、群论、复数
拉格朗日插值
有限差分法
- 约束：步长为1
- 时间复杂度 O ( n 2 ) O(n^2) O(n2)
每一次对一个n阶多项式做差分，相当于将它求导，经过n次差分之后，我们会得到一组定值。此时，等于 n ! × a n n!\times a_n n!×an。因此我们可以得出 a n a_n an的值。现在，我们把每一对都减去这个 a n a_n an的贡献；再算 n − 1 ! × a n − 1 {n-1}!\times a_{n-1} n−1!×an−1。。。就能把所有的系数求出来。
拉格朗日插值
- 构造
- 时间复杂度 O ( n 2 ) O(n^2) O(n2)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll mod = 998244353;
const int maxn = 2e3+4;
ll x[maxn],y[maxn];
ll quick_p(ll a,ll n){
a%=mod;
ll ans = 1;
while(n){
if(n&1) ans = ans*a%mod;
a=a*a%mod;
n>>=1;
}
return ans;
}
int main(){
int n,k;
scanf("%d%d",&n,&k);
for(int i=1;i<=n;++i){
scanf("%d%d",&x[i],&y[i]);
}
ll ans = 0;
for(int i=1;i<=n;++i){
ll s1=1,s2=1;
for(int j=1;j<=n;++j){
if(i==j) continue;
s1 = s1*(k-x[j])%mod;
s2 = s2*(x[i]-x[j])%mod;
}
ans += y[i]*s1*quick_p(s2,mod-2)%mod;
ans%=mod;
}
cout<<(ans+mod)%mod<<endl;
}
references
快速傅里叶变换(FFT)
- 时间复杂度 O ( n l g n ) O(nlgn) O(nlgn)
- 没什么局限的算法
- 蝴蝶变换（位逆序变换）
- DFT和IDFT过程
- 递归
- 单位复根的特性
- Vandermonde矩阵的逆矩阵
#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e6+5;
const double PI = acos(-1.0);
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x;
y = _y;
}
Complex operator-(const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator+(const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator*(const Complex &b) const {
return Complex(x * b.x - y * b.y,x * b.y + y * b.x);
}
};
int rev[maxn<<2];
// O(n) 位逆序置换
void change(Complex y[], int len) {
for (int i = 0; i < len; ++i) {
rev[i] = rev[i>>1] >> 1;
if (i & 1) rev[i] |= len >> 1;
}
for (int i = 0; i < len; ++i)
if (i < rev[i]) swap(y[i], y[rev[i]]);
}
// 非递归fft
void fft(Complex y[], int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j += h) {
Complex w(1, 0);
for (int k = j; k < j + h / 2; ++k) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1)
for (int i = 0; i < len; ++i) y[i].x /= len;
}
Complex f[maxn<<2], g[maxn<<2];
int main() {
int n,m;scanf("%d%d",&n,&m);
int len = 1;
while (len < (n+1) * 2 || len < (m+1) * 2) len <<= 1;
for (int i = 0; i < n + 1; ++i) {
int tmp;
scanf("%d",&tmp);
f[i] = Complex(tmp,0);
}
for (int i=n+1;i<len;++i) f[i] = Complex(0,0);
for (int i=0;i<m+1;++i) {
int tmp;
scanf("%d",&tmp);
g[i] = Complex(tmp,0);
}
for (int i=m+1;i<len;++i) g[i] = Complex(0,0);
fft(f,len,1); fft(g,len,1);
for (int i=0;i<len;++i) f[i] = f[i]*g[i];
fft(f,len,-1);
for (int i=0; i<n+m+1; ++i)
cout<<int(f[i].x+0.5)<<" ";
}
洛谷3338 力
纪念写的第一道fft。
首先
F
j
=
∑
i
=
1
j
−
1
q
i
(
i
−
j
)
2
−
∑
i
=
j
+
1
n
q
i
(
i
−
j
)
2
F_j = \sum_{i=1}^{j-1}\frac{q_i}{(i-j)^2}-\sum_{i=j+1}^{n}\frac{q_i}{(i-j)^2}
Fj=∑i=1j−1(i−j)2qi−∑i=j+1n(i−j)2qi
加多一项，前减去一个后加上一个，不变化
F
j
=
∑
i
=
1
j
q
i
(
i
−
j
)
2
−
∑
i
=
j
n
q
i
(
i
−
j
)
2
F_j =\sum_{i=1}^j\frac{q_i}{(i-j)^2}-\sum_{i=j}^n\frac{q_i}{(i-j)^2}
Fj=∑i=1j(i−j)2qi−∑i=jn(i−j)2qi
设两个函数。设他们下标为0的情况，扩展为0到j。然后可得前一个部分是一个卷积。
后一个也写成这两个函数的样子。然后把枚举顺序从j到n改成从0到n-j。把i也改成i+j。
得到新的公式。把f倒转过来。得卷积。
然后答案就得到了，带进去算即可。
#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e5+5;
const double PI = acos(-1.0);
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x;
y = _y;
}
Complex operator-(const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator+(const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator*(const Complex &b) const {
return Complex(x * b.x - y * b.y,x * b.y + y * b.x);
}
};
int rev[maxn<<2];
// O(n) 位逆序置换
void change(Complex y[], int len) {
for (int i = 0; i < len; ++i) {
rev[i] = rev[i>>1] >> 1;
if (i & 1) rev[i] |= len >> 1;
}
for (int i = 0; i < len; ++i)
if (i < rev[i]) swap(y[i], y[rev[i]]);
}
// 非递归fft
void fft(Complex y[], int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j += h) {
Complex w(1, 0);
for (int k = j; k < j + h / 2; ++k) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1)
for (int i = 0; i < len; ++i) y[i].x /= len;
}
Complex f[maxn<<2], g[maxn<<2], f2[maxn<<2];
double q[maxn];
int main() {
int n;scanf("%d",&n);
f[0] = Complex(0,0);
g[0] = Complex(0,0);
f2[0] = Complex(0,0);
// 长度为n+1。0~n。
for(int i=1;i<=n;++i){
scanf("%lf",&q[i]);
f[i] = Complex(q[i],0);
g[i] = Complex((double)1/i/i,0);
f2[n-i]=Complex(q[i],0);
}
int len = 1;
while (len < (n+1) * 2 ) len <<= 1;
for (int i=n+1;i<len;++i) {
f[i] = Complex(0, 0);
g[i] = Complex(0,0);
f2[i] = Complex(0,0);
}
fft(f,len,1); fft(g,len,1);
fft(f2,len,1);
for (int i=0;i<len;++i) {
f[i] = f[i]*g[i];
f2[i] = f2[i]*g[i];
}
fft(f,len,-1);
fft(f2,len,-1);
for (int i=1; i<n+1; ++i)
printf("%.3lf\n",f[i].x-f2[n-i].x);
}
NTT
为了学NTT来看的这道题。但是任意模数乘法现在好像都不流行用NTT了，都流行用FFT优化(MTT?)。但不管怎么说，NTT还是要学。
NTT的题目就是限定在了模域之下。此时我们再用FFT来写会发现精度不够（而且会爆）。天才数学家们发现了单位复根的那些特殊性质，模域下的原根也有。NTT就是原根版本的FFT。
原根与单位复根的性质对照
但是我们发现NTT只能用于特殊的模数如998244353.这有很大的限制。
如果模数为任意，怎么办？
我们这时可以考虑中国剩余定理，用三个很大(乘积大于最大的值，故不存在取模)满足NTT要求的质数来算，然后合并结果。
但是CRT直接算会爆ll。毕竟这三个p很大。我们就分步来算。
任意模数NTT介绍
#include<bits/stdc++.h>
using namespace std;
int mod;
namespace Math {
// 一种神奇的快速幂
inline int pw(int base,int p,const int mod){
static int res;
for(res=1;p;p>>=1,base=static_cast<long long> (base)*base%mod) if(p&1) res = static_cast<long long> (res)*base%mod;
return res;
}
inline int inv(int x,const int mod) {return pw(x,mod-2,mod);}
}
const int mod1 = 998244353, mod2 = 1004535809, mod3 = 469762049, G = 3;// 原根均为3
const long long mod_1_2 = static_cast<long long>(mod1)*mod2;
const int inv_1 = Math::inv(mod1,mod2),inv_2 = Math::inv(mod_1_2%mod3,mod3); // 两个逆元求出
struct Int {
int A,B,C;
// 三个模意义下的值。新建结构体，重载运算。
explicit inline Int() {}
explicit inline Int(int __num) :A(__num),B(__num),C(__num){}
explicit inline Int(int __A,int __B,int __C):A(__A),B(__B),C(__C){}
static inline Int reduce(const Int &x){
return Int(x.A+(x.A>>31&mod1),x.B+(x.B>>31&mod2),x.C+(x.C>>31&mod3));
}
inline friend Int operator + (const Int &lhs,const Int &rhs){
return reduce(Int(lhs.A+rhs.A-mod1,lhs.B+rhs.B-mod2,lhs.C+rhs.C-mod3));
}
inline friend Int operator - (const Int &lhs,const Int &rhs){
return reduce(Int(lhs.A-rhs.A,lhs.B-rhs.B,lhs.C-rhs.C));
}
inline friend Int operator * (const Int &lhs,const Int &rhs){
return Int(static_cast<long long>(lhs.A)*rhs.A%mod1,static_cast<long long>(lhs.B)*rhs.B%mod2,static_cast<long long>(lhs.C)*rhs.C%mod3);
}
inline int get(){
long long x = static_cast<long long>(B-A+mod2)%mod2*inv_1%mod2*mod1+A;
return (static_cast<long long> (C-x%mod3+mod3)%mod3*inv_2%mod3*(mod_1_2%mod)%mod+x)%mod;
}
};
#define maxn 140004 // 不知道为什么要开这么大
namespace Poly{
#define N (maxn<<1)
int lim,s,rev[N];
Int Wn[N|1];
inline void init(int n){
s=-1,lim=1;while(lim<n)lim<<=1,++s;
for(register int i=1;i<lim;++i) rev[i]=rev[i>>1]>>1|(i&1)<<s; // 位逆序
const Int t(Math::pw(G,(mod1-1)/lim,mod1),Math::pw(G,(mod2-1)/lim,mod2),Math::pw(G,(mod3-1)/lim,mod3));
*Wn = Int(1); for(register Int *i=Wn;i!=Wn+lim;++i) *(i+1)=*i*t; // 计算原版和逆元版(单位复根)
}
inline void NTT(Int *A,const int op=1){
for(register int i=1;i<lim;++i) if(i<rev[i]) std::swap(A[i],A[rev[i]]);
for(register int mid=1;mid<lim;mid<<=1){
const int t=lim/mid>>1;
for(register int i=0;i<lim;i+=mid<<1){
for(register int j=0;j<mid;++j){
const Int W = op?Wn[t*j]:Wn[lim-t*j];
const Int X = A[i+j],Y = A[i+j+mid]*W;
A[i+j]=X+Y,A[i+j+mid]=X-Y;
}
}
}
if(!op){
const Int ilim(Math::inv(lim,mod1),Math::inv(lim,mod2),Math::inv(lim,mod3));
for(register Int *i=A;i!=A+lim;++i) *i=(*i)*ilim;
}
}
#undef N
}
int n,m;
Int A[maxn<<1],B[maxn<<1];
int main(){
scanf("%d%d%d",&n,&m,&mod);++n,++m;
for(int i=0,x;i<n;++i) scanf("%d",&x),A[i]=Int(x%mod);
for(int i=0,x;i<m;++i) scanf("%d",&x),B[i]=Int(x%mod);
Poly::init(n+m);
Poly::NTT(A),Poly::NTT(B);
for(int i=0;i<Poly::lim;++i) A[i] = A[i]*B[i];
Poly::NTT(A,0);
for(int i=0;i<n+m-1;++i){
printf("%d",A[i].get());
putchar(i==n+m-2?'\n':' ');
}
}
高精度
POJ1737 Connected Graph有标号简单无向连通图计数
每个点是不一样的。
这道题可以用指数生成函数EGF来写，貌似需要多项式求逆之类的。。因为是高精。
但也可以考虑公式递推。
此处的排列组合数如何计算呢？
为了防止重复计数或遗漏，我们考虑只关注n个点的其中某个点，let’s say it is A的情况。仅对于该点的情况进行分类讨论。
我们把问题：寻找无向连通图的种类数。转化为任意图的种类减去不连通的图的种类数。
那么对于该点，必然属于一个大小从1到n-1（设为k）的连通块。这个连通块内的其他点是在n-1个中选择k-1个。这个连通块之外的点是任意图。
那么设f为无向连通图个数，h为任意图个数，g为不连通的图个数，得到公式
f
(
n
)
=
h
(
n
)
−
g
(
n
)
f(n) = h(n)-g(n)
f(n)=h(n)−g(n)
h
(
n
)
=
2
(
n
−
1
)
n
2
h(n)=2^{\frac{(n-1)n}{2}}
h(n)=22(n−1)n
g
(
n
)
=
∑
i
=
1
n
−
1
(
n
−
1
i
−
1
)
f
(
i
)
h
(
n
−
i
)
g(n)=\sum_{i=1}^{n-1}\tbinom{n-1}{i-1}f(i)h(n-i)
g(n)=i=1∑n−1(i−1n−1)f(i)h(n−i)
显然2^250次幂会爆炸大，我们要把该递推式化简。
下面板子是一种vector写法，在使用上述公式处理的话会出现RE。问题在赋值函数处。
所以使用了一个生成函数的公式，即
f
(
n
)
=
∑
i
=
1
n
−
1
(
n
−
2
i
−
1
)
f
(
i
)
f
(
n
−
i
)
(
2
i
−
1
)
f(n)=\sum_{i=1}^{n-1}\tbinom{n-2}{i-1}f(i)f(n-i)(2^{i}-1)
f(n)=i=1∑n−1(i−1n−2)f(i)f(n−i)(2i−1)
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<vector>
typedef long long ll;
const int maxn = 400;
int n;
struct BigInt{
std::vector<ll> v;
BigInt(){*this = 0;}
BigInt(int x){*this = x;}
BigInt &operator=(int x){
v.clear();
do v.push_back(x%10);while(x/=10);
return *this;
}
BigInt &operator=(const BigInt &x){
v.resize(x.v.size());
memcpy(const_cast<ll *>(v.data()),x.v.data(),x.v.size()*sizeof(ll));
return *this;
}
friend BigInt operator+(const BigInt &a,const BigInt &b){
BigInt result;
result.v.clear();
bool flag = false;
for(int i=0;i<(int)std::max(a.v.size(),b.v.size());++i){
int tmp = 0;
if(i<(int)a.v.size()) tmp+=a.v[i];
if(i<(int)b.v.size()) tmp+=b.v[i];
if(flag) tmp++,flag=false;
if(tmp>=10) tmp-=10,flag=true;
result.v.push_back(tmp);
}
if(flag) result.v.push_back(1);
return result;
}
friend BigInt &operator+=(BigInt &a,const BigInt &b){
return a = a+b;
}
friend BigInt operator-(const BigInt &a,const BigInt &b){
BigInt result;
result.v.clear();
bool flag = false;
for(int i=0;i<(int)a.v.size();++i){
int tmp = a.v[i];
if(i<(int)b.v.size()) tmp -= b.v[i];
if(flag) tmp--,flag = false;
if(tmp<0) tmp+=10,flag=true;
result.v.push_back(tmp);
}
int size = result.v.size();
while(size>1&&result.v[size-1]==0) size--;
result.v.resize(size);
return result;
}
friend BigInt &operator-=(BigInt &a,const BigInt &b){
return a = a-b;
}
friend BigInt operator*(const BigInt &a,const BigInt &b){
BigInt temp,result;
temp.v.clear();
temp.v.resize((int)a.v.size()+(int)b.v.size()-1);
result.v.clear();
int flag = 0;
for(int i=0;i<(int)a.v.size();++i)
for(int j=0;j<(int)b.v.size();++j)
temp.v[i+j] += a.v[i]*b.v[j];
for(int i=0;i<(int)temp.v.size();++i){
int tmp = temp.v[i] + flag;
flag = tmp/10;
result.v.push_back(tmp%10);
}
while(flag){
result.v.push_back(flag%10);
flag /= 10;
}
return result;
}
friend BigInt operator*(const BigInt &a,int &b){
BigInt temp,result;
result.v.clear();
temp.v.clear();
int flag = 0;
for(int i=0;i<(int)a.v.size();++i) temp.v.push_back(a.v[i]*b);
for(int i=0;i<(int)temp.v.size();++i) {
int tmp = temp.v[i];
tmp += flag;
flag = tmp/10;
result.v.push_back(tmp%10);
}
while(flag){
result.v.push_back(flag%10);
flag /= 10;
}
return result;
}
};
// 重载输出big integer
std::ostream &operator<<(std::ostream &out, const BigInt &x){
for(int i=x.v.size()-1;i>=0;--i) out << (char)(x.v[i]+'0');
return out;
}
const int MAXN = 50;
BigInt combo[MAXN+1][MAXN+1];
BigInt bin[MAXN+1];
inline void makeComboTable(){
for(int i=1;i<=MAXN;++i){
combo[i][0] = combo[i][i] = 1;
for(int j=1;j<i;++j){
combo[i][j] = combo[i-1][j] + combo[i-1][j-1];
}
}
}
inline BigInt &C(int n,int k){
return combo[n][k];
}
inline void makePowerTwoTable(){
bin[0] = 1;
for(int i=1;i<=MAXN;++i) bin[i] = bin[i-1]*2;
}
BigInt f[MAXN+1];
int main(){
makeComboTable();
makePowerTwoTable();
f[1] = 1;f[2] = 1;
for(int i=3;i<=MAXN;++i){
for(int j=1;j<i;++j){
f[i] += C(i-2,j-1)*f[i-j]*f[j]*(bin[j]-1);
}
}
int n;
while(~scanf("%d",&n),n){
std::cout<<f[n]<<std::endl;
}
}
还是整点阳间的板子吧（数组模拟）
#include <string>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef int ll;
inline ll read(){
ll s=0; bool f=0; char ch=' ';
while(!isdigit(ch)) {f|=(ch=='-'); ch=getchar();}
while(isdigit(ch)) {s=(s<<3)+(s<<1)+(ch^48); ch=getchar();}
return (f)?(-s):(s);
}
#define R(x) x=read()
inline void write(ll x){
if(x<0) {putchar('-'); x=-x;}
if(x<10) {putchar(x+'0'); return;}
write(x/10); putchar((x%10)+'0');
}
#define W(x) write(x),putchar(' ')
#define Wl(x) write(x),putchar('\n')
const int power=4,Base=10000;
struct Int{
int a[205];
Int() {memset(a,0,sizeof a);}
Int(int x){
memset(a,0,sizeof a);
while(x>0){
a[++a[0]]=x%Base; x/=Base;
}
}
inline void print(){
int i;
write(a[a[0]]);
for(i=a[0]-1;i>=1;i--){
if(a[i]<1000) putchar('0');
if(a[i]<100) putchar('0');
if(a[i]<10) putchar('0');
write(a[i]);
}
}
}Bin[2005],C[55][55],Ans[55];
#define Pl(x) x.print(),putchar('\n')
inline Int operator+(Int p,Int q)
{
int i; Int res=p; res.a[0]=max(p.a[0],q.a[0]);
for(i=1;i<=q.a[0];i++){
res.a[i]+=q.a[i];
res.a[i+1]+=res.a[i]/Base;
res.a[i]%=Base;
}
while(res.a[res.a[0]+1]) res.a[0]++;
return res;
}
inline Int operator-(Int p,Int q){
int i; Int res=p;
for(i=1;i<=q.a[0];i++){
res.a[i]-=q.a[i];
if(res.a[i]<0){
res.a[i+1]--; res.a[i]+=Base;
}
}
while(!res.a[res.a[0]]) res.a[0]--;
return res;
}
inline Int operator*(Int p,int q){
int i; Int res=p;
for(i=1;i<=res.a[0];i++) res.a[i]*=q;
for(i=1;i<=res.a[0];i++){
res.a[i+1]+=res.a[i]/Base; res.a[i]%=Base;
}
while(res.a[res.a[0]+1]){
res.a[0]++; res.a[res.a[0]+1]+=res.a[res.a[0]]/Base; res.a[res.a[0]]%=Base;
}
return res;
}
inline Int operator*(Int p,Int q){
int i,j; Int res; res.a[0]=p.a[0]+q.a[0];
for(i=1;i<=p.a[0];i++) for(j=1;j<=q.a[0];j++){
res.a[i+j-1]+=p.a[i]*q.a[j];
res.a[i+j]+=res.a[i+j-1]/Base;
res.a[i+j-1]%=Base;
}
while(!res.a[res.a[0]]) res.a[0]--;
return res;
}
int main(){
int i,j,n;
Bin[0]=Int(1); for(i=1;i<=1500;i++) Bin[i]=Bin[i-1]*2;
C[0][0]=Int(1);
for(i=1;i<=50;i++){
C[i][0]=Int(1);
for(j=1;j<=i;j++) C[i][j]=C[i-1][j]+C[i-1][j-1];
}
Ans[1]=Int(1);
for(i=2;i<=50;i++){
Ans[i]=Bin[i*(i-1)/2];
for(j=1;j<i;j++){
Ans[i]=Ans[i]-Ans[j]*C[i-1][j-1]*Bin[(i-j)*((i-j)-1)/2];
}
}
for(;;){
R(n); if(n==0) break; Pl(Ans[n]);
}
return 0;
}
多项式求逆
公式推导见OI wiki。
注意我们从(n+1)/2的子问题处理n的子问题时，是对一个取了(n+1)/2的模的子问题逆元和长度为2*n的原函数f(why？因为当前是模n意义下。f的长度肯定长于n。那么一个这样的多项式与模(n+1)/2的多项式乘积模n，就将两者补全为长度为2n的多项式进行计算，然后再取模)进行递推。
板子：
// 多项式求逆
#include<bits/stdc++.h>
using namespace std;
const int MAXN = 1e5+5;
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch = getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}
const int P = 998244353;
inline int qpow(int x,int y){int res(1);while(y){if(y&1) res = 1ll*res*x%P;x = 1ll*x*x%P;y>>=1;}return res;}
int rev[MAXN<<2];
void change(int y[],int len){
for(int i=0;i<len;++i){
rev[i] = rev[i>>1]>>1;
if(i&1) rev[i] |= len>>1;
}
for(int i=0;i<len;++i)
if(i<rev[i]) swap(y[i],y[rev[i]]);
}
void ntt(int y[],int len,int on){
change(y,len);
for(int h=2;h<=len;h<<=1){
int gn = qpow(3,(P-1)/h);
for(int j=0;j<len;j+=h){
int g = 1;
for(int k=j;k<j+h/2;++k){
int u = y[k];
int t = 1ll*g*y[k+h/2]%P;
y[k] = (u+t)%P;
y[k+h/2] = (u-t+P)%P;
g = 1ll*g*gn%P;
}
}
}
if(on==-1) {
reverse(y+1,y+len);
register int inv = qpow(len,P-2);
for(int i=0;i<len;++i) y[i] = 1ll*y[i]*inv%P;
}
}
static int inv_t[MAXN<<2];
void solve(int f[],int g[],int n){
// 递归:把两个n/2的合并成一个n的
if(n==1) {f[0] = qpow(g[0],P-2);return;}
solve(f,g,(n+1)>>1);
int len(1);
while(len<(n<<1)) len<<=1;
copy(g,g+n,inv_t);
fill(inv_t+n,inv_t+len,0);
ntt(f,len,1);
ntt(inv_t,len,1);
for(int i=0;i<len;++i){
f[i] = 1ll*f[i]*(2-1ll*f[i]*inv_t[i]%P+P)%P;
}
ntt(f,len,-1);
fill(f+n,f+len,0); // 取模
}
void polyinv(int f[], int h[], const int n) {
fill(f, f + n + n, 0);
f[0] = qpow(h[0], P - 2);
for (int t = 2; t <= n; t <<= 1) {
const int t2 = t << 1;
copy(h, h + t, inv_t);
fill(inv_t + t, inv_t + t2, 0);
ntt(f, t2, 1);
ntt(inv_t, t2, 1);
for (int i = 0; i != t2; ++i)
f[i] = 1ll*f[i]*(2-1ll*f[i]*inv_t[i]%P+P)%P;
ntt(f, t2, -1);
fill(f + t, f + t2, 0);
}
}
int f[MAXN<<2],g[MAXN<<2];
int main(){
//freopen("P4238_16.in","r",stdin);
int n=read();
for(int i=0;i<n;++i) g[i]=read();
//cout<<g[i]<<" ";
solve(f,g,n);
//int len(1);
//while(n>len) len<<=1;
//polyinv(f,g,len);
for(int i=0;i<n;++i)
printf("%d%c",f[i],i==n-1?'\n':' ');
}
多项式开根
洛谷板子
跟多项式求逆异曲同工。
// 多项式开根
#include<bits/stdc++.h>
using namespace std;
const int MAXN = 1e5+5;
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch = getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}
const int P = 998244353;
inline int qpow(int x,int y){int res(1);while(y){if(y&1) res = 1ll*res*x%P;x = 1ll*x*x%P;y>>=1;}return res;}
int rev[MAXN<<2];
void change(int y[],int len){
for(int i=0;i<len;++i){
rev[i] = rev[i>>1]>>1;
if(i&1) rev[i] |= len>>1;
}
for(int i=0;i<len;++i)
if(i<rev[i]) swap(y[i],y[rev[i]]);
}
void ntt(int y[],int len,int on){
change(y,len);
for(int h=2;h<=len;h<<=1){
int gn = qpow(3,(P-1)/h);
for(int j=0;j<len;j+=h){
int g = 1;
for(int k=j;k<j+h/2;++k){
int u = y[k];
int t = 1ll*g*y[k+h/2]%P;
y[k] = (u+t)%P;
y[k+h/2] = (u-t+P)%P;
g = 1ll*g*gn%P;
}
}
}
if(on==-1) {
reverse(y+1,y+len);
register int inv = qpow(len,P-2);
for(int i=0;i<len;++i) y[i] = 1ll*y[i]*inv%P;
}
}
static int sqrt_t[MAXN<<2],inv_t[MAXN<<2];
int invf[MAXN<<2];
int inv_two;
void getinv(int f[],int g[],int n){
// 递归:把两个n/2的合并成一个n的
if(n==1) {f[0] = qpow(g[0],P-2);return;}
getinv(f,g,(n+1)>>1);
int len(1);
while(len<n<<1) len<<=1;
copy(g,g+n,inv_t);
fill(inv_t+n,inv_t+len,0);
ntt(f,len,1);
ntt(inv_t,len,1);
for(int i=0;i<len;++i){
f[i] = 1ll*f[i]*(2-1ll*f[i]*inv_t[i]%P+P)%P;
}
ntt(f,len,-1);
fill(f+n,f+len,0); // 取模
}
void solve(int f[],int g[],int n){
// 递归:把两个n/2的合并成一个n的
if(n==1) {f[0] = 1;return;} // 不需要二次剩余因为题目保证了a0=1
solve(f,g,(n+1)>>1);
int len(1);
while(len<n<<1) len<<=1;
copy(g,g+n,sqrt_t);
fill(sqrt_t+n,sqrt_t+len,0);
fill(invf,invf+len,0);
getinv(invf,f,n);
ntt(f,len,1);
ntt(invf,len,1);
ntt(sqrt_t,len,1);
for(int i=0;i<len;++i){
f[i] = (1ll*inv_two*f[i]%P+1ll*inv_two*invf[i]%P*sqrt_t[i]%P)%P;
}
ntt(f,len,-1);
fill(f+n,f+len,0); // 取模
}
int f[MAXN<<2],g[MAXN<<2];
int main(){
int n=read();
for(int i=0;i<n;++i) g[i]=read();
inv_two = qpow(2,P-2);
solve(f,g,n);
for(int i=0;i<n;++i)
printf("%d%c",f[i],i==n-1?'\n':' ');
}
多项式除法
多项式的带余数除法。
推公式就行。。用了一种神奇的思路。利用了取模。
模板
// 多项式除法
#include<bits/stdc++.h>
using namespace std;
const int MAXN = 1e5+5;
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch = getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}
const int P = 998244353;
inline int qpow(int x,int y){int res(1);while(y){if(y&1) res = 1ll*res*x%P;x = 1ll*x*x%P;y>>=1;}return res;}
int rev[MAXN<<2];
void change(int y[],int len){
for(int i=0;i<len;++i){
rev[i] = rev[i>>1]>>1;
if(i&1) rev[i] |= len>>1;
}
for(int i=0;i<len;++i)
if(i<rev[i]) swap(y[i],y[rev[i]]);
}
void ntt(int y[],int len,int on){
change(y,len);
for(int h=2;h<=len;h<<=1){
int gn = qpow(3,(P-1)/h);
for(int j=0;j<len;j+=h){
int g = 1;
for(int k=j;k<j+h/2;++k){
int u = y[k];
int t = 1ll*g*y[k+h/2]%P;
y[k] = (u+t)%P;
y[k+h/2] = (u-t+P)%P;
g = 1ll*g*gn%P;
}
}
}
if(on==-1) {
reverse(y+1,y+len);
register int inv = qpow(len,P-2);
for(int i=0;i<len;++i) y[i] = 1ll*y[i]*inv%P;
}
}
static int inv_t[MAXN<<2];
void solve(int f[],int g[],int n){
// 递归:把两个n/2的合并成一个n的
if(n==1) {f[0] = qpow(g[0],P-2);return;}
solve(f,g,(n+1)>>1);
int len(1);
while(len<(n<<1)) len<<=1;
copy(g,g+n,inv_t);
fill(inv_t+n,inv_t+len,0);
ntt(f,len,1);
ntt(inv_t,len,1);
for(int i=0;i<len;++i){
f[i] = 1ll*f[i]*(2-1ll*f[i]*inv_t[i]%P+P)%P;
}
ntt(f,len,-1);
fill(f+n,f+len,0); // 取模
}
int f[MAXN<<2], g[MAXN<<2];
int f2[MAXN<<2], g2[MAXN<<2]; // 一个未取模的f
int invgr[MAXN<<2], q[MAXN<<2], r[MAXN<<2];
int main(){
int n=read(),m=read();
n++,m++;
for(int i=0;i<n;++i) f[i]=read(),f2[i] = f[i];
for(int i=0;i<m;++i) g[i]=read(),g2[i] = g[i];
reverse(f,f+n);
reverse(g,g+m);
for(int i=n-m+1;i<m;++i) g[i] = 0; // 取模
fill(f+n-m+1,f+n,0); // 取模
solve(invgr,g,n-m+1);
int len(1);
while(len<((n-m+1)<<1)) len<<=1;
ntt(invgr,len,1), ntt(f,len,1);
for(int i=0;i<len;++i) q[i] = 1ll*invgr[i] * f[i] % P;
ntt(q,len,-1);
for(int i=0;i<n-m+1;++i){
int j=n-m-i;
if(i<j) swap(q[i],q[j]);
printf("%d%c",q[i],i==n-m?'\n':' ');
}
for(int i=n-m+1;i<n;++i) q[i] = 0;
// 不晓得为啥要清空，不是说q的阶小于n-m+1嘛
while(len<(n<<1)) len<<=1;
while(len<(m<<1)) len<<=1;
ntt(f2,len,1), ntt(g2,len,1), ntt(q,len,1);
// q要不要重新ntt呢？len改变了，所以我取了。。
for(int i=0;i<len;++i)
r[i] = (1ll*f2[i] - 1ll*g2[i]*q[i]%P + P)%P;
ntt(r,len,-1);
for(int i=0;i<m-1;++i)
printf("%d%c",r[i],i==m-2?'\n':' ');
}
分治FFT
板子
法一：生成函数+多项式求逆
We have
f
[
i
]
=
∑
j
=
1
i
f
[
i
−
j
]
×
g
[
j
]
,
i
>
0
f[i]=\sum_{j=1}^i f[i-j]\times g[j],i>0
f[i]=∑j=1if[i−j]×g[j],i>0.
define
f
(
x
)
=
∑
i
=
0
∞
f
[
i
]
x
i
f(x)=\sum_{i=0}^\infty f[i]x^i
f(x)=∑i=0∞f[i]xi
and
g
(
x
)
=
∑
i
=
0
∞
g
[
i
]
x
i
g(x)=\sum_{i=0}^\infty g[i]x^i
g(x)=∑i=0∞g[i]xi.
(let
g
[
0
]
=
0
g[0]=0
g[0]=0)
convolute:
f
(
x
)
∗
g
(
x
)
=
∑
i
=
0
∞
∑
j
=
0
∞
f
[
i
]
g
[
j
]
x
i
+
j
f(x)*g(x)=\sum_{i=0}^\infty \sum_{j=0}^\infty f[i] g[j]x^{i+j}
f(x)∗g(x)=∑i=0∞∑j=0∞f[i]g[j]xi+j
let
k
=
i
+
j
k=i+j
k=i+j
f
(
x
)
∗
g
(
x
)
=
∑
k
=
0
∞
(
∑
j
=
0
k
f
[
k
−
j
]
g
[
j
]
)
x
k
f(x)*g(x)=\sum_{k=0}^\infty( \sum_{j=0}^k f[k-j]g[j])x^k
f(x)∗g(x)=∑k=0∞(∑j=0kf[k−j]g[j])xk
when
k
>
0
k>0
k>0,
∑
j
=
0
k
f
[
k
−
j
]
g
[
j
]
=
f
[
k
]
\sum_{j=0}^k f[k-j]g[j]=f[k]
∑j=0kf[k−j]g[j]=f[k]
when
k
=
0
k=0
k=0,
f
[
0
]
g
[
0
]
=
0
f[0]g[0]=0
f[0]g[0]=0,
(
g
[
0
]
=
0
)
(g[0]=0)
(g[0]=0)
therefore,
f
(
x
)
∗
g
(
x
)
=
∑
k
=
1
∞
f
[
k
]
x
k
=
f
(
x
)
−
f
[
0
]
f(x)*g(x)=\sum_{k=1}^\infty f[k]x^k=f(x)-f[0]
f(x)∗g(x)=∑k=1∞f[k]xk=f(x)−f[0]
we have
f
(
x
)
=
f
[
0
]
1
−
g
(
x
)
f(x)=\frac{f[0]}{1-g(x)}
f(x)=1−g(x)f[0].(多项式求逆)
#include<bits/stdc++.h>
using namespace std;
const int MAXN = 1e5+5;
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch = getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}
const int P = 998244353;
inline int qpow(int x,int y){int res(1);while(y){if(y&1) res = 1ll*res*x%P;x = 1ll*x*x%P;y>>=1;}return res;}
int rev[MAXN<<2];
void change(int y[],int len){
for(int i=0;i<len;++i){
rev[i] = rev[i>>1]>>1;
if(i&1) rev[i] |= len>>1;
}
for(int i=0;i<len;++i)
if(i<rev[i]) swap(y[i],y[rev[i]]);
}
void ntt(int y[],int len,int on){
change(y,len);
for(int h=2;h<=len;h<<=1){
int gn = qpow(3,(P-1)/h);
for(int j=0;j<len;j+=h){
int g = 1;
for(int k=j;k<j+h/2;++k){
int u = y[k];
int t = 1ll*g*y[k+h/2]%P;
y[k] = (u+t)%P;
y[k+h/2] = (u-t+P)%P;
g = 1ll*g*gn%P;
}
}
}
if(on==-1) {
reverse(y+1,y+len);
register int inv = qpow(len,P-2);
for(int i=0;i<len;++i) y[i] = 1ll*y[i]*inv%P;
}
}
static int inv_t[MAXN<<2];
void solve(int f[],int g[],int n){
// 递归:把两个n/2的合并成一个n的
if(n==1) {f[0] = qpow(g[0],P-2);return;}
solve(f,g,(n+1)>>1);
int len(1);
while(len<(n<<1)) len<<=1;
copy(g,g+n,inv_t);
fill(inv_t+n,inv_t+len,0);
ntt(f,len,1);
ntt(inv_t,len,1);
for(int i=0;i<len;++i){
f[i] = 1ll*f[i]*(2-1ll*f[i]*inv_t[i]%P+P)%P;
}
ntt(f,len,-1);
fill(f+n,f+len,0); // 取模
}
int invg[MAXN<<2],g[MAXN<<2];
int main(){
int n=read();
for(int i=1;i<n;++i) g[i]=read(),g[i]=(P-g[i])%P;
g[0]=1;
solve(invg,g,n);
for(int i=0;i<n;++i){
printf("%d%c",invg[i],(i==n-1)?'\n':' ');
}
}
6389
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
