In mathematics you don't understand things, you just get used to them.

标签 polynomials 下的文章

Description

Link.

给你一个序列,你每次可以取 $1\sim3$ 个数然后计算和,问你对于每一种和,方案数是多少。

Solution

设一个 OGF $A(x)=\sum_{i=0}^{+\infty}a_{i}x^{i}$,指数为物品的价值,$a_{i}$ 为出现的次数。

也就是说,$\sum a_{i}$ 就是选个一个数的答案。

再来考虑选两个数。$A^{2}(x)$ 显然会有重复。

重复有两种情况,第一种是 $i\times j$ 和 $j\times i$,这个简单,除个 $2!$ 即可。

还有就是 $i\times i$ 的情况,那么再设一个表示一个斧头选两次的 OGF $B(x)=\sum_{i=0}^{+\infty}b_{i}x^{2i}$。

那么选两个数的答案为 $\frac{A^{2}(x)-B(x)}{2!}$。

再来考虑选三个数,基本同理地设 $C(x)=\sum_{i=0}^{+\infty}c_{i}x^{3i}$。

然后选三个数的答案为 $\frac{A^{3}(x)-3A(x)B(x)+2C(x)}{3!}$,这个容斥就是考虑 $\sf aab,baa$ 和 $\sf abc,acb$ 情况的重复。

做完了。

#include<bits/stdc++.h>
using namespace std;
namespace Poly
{
    typedef complex<double> comp;
    typedef vector<complex<double> > poly;
    #define len(x) (int((x).size()))
    const double bh_pi=acos(-1),eps=1e-3;
    int lim,rev[500010];
    void fft(poly &f,int op)
    {
        for(int i=0;i<lim;++i)    if(i<rev[i])    swap(f[i],f[rev[i]]);
        for(int len=2;len<=lim;len<<=1)
        {
            comp bas(cos(2*bh_pi/len),op*sin(2*bh_pi/len));
            for(int fr=0;fr<lim;fr+=len)
            {
                comp now(1,0);
                for(int ba=fr;ba<fr+(len>>1);++ba,now*=bas)
                {
                    comp tmp=now*f[ba+(len>>1)];
                    f[ba+(len>>1)]=f[ba]-tmp;
                    f[ba]+=tmp;
                }
            }
        }
        if(op==-1)    for(int i=0;i<lim;++i)    f[i]/=lim;
    }
    poly mulPoly(poly f,poly g)
    {
        int n=len(f)+len(g)-1;
        for(lim=1;lim<n;lim<<=1);
        for(int i=0;i<lim;++i)    rev[i]=(rev[i>>1]>>1)|((i&1)?(lim>>1):0);
        f.resize(lim),g.resize(lim);
        fft(f,1),fft(g,1);
        for(int i=0;i<lim;++i)    f[i]*=g[i];
        fft(f,-1),f.resize(n);
        return f;
    }
    poly addPoly(poly f,poly g)
    {
        int n=max(len(f),len(g));
        f.resize(n),g.resize(n);
        for(int i=0;i<n;++i)    f[i]+=g[i];
        return f;
    }
    poly decPoly(poly f,poly g)
    {
        int n=max(len(f),len(g));
        f.resize(n),g.resize(n);
        for(int i=0;i<n;++i)    f[i]-=g[i];
        return f;
    }
    poly mulPoly(poly f,double x)
    {
        for(int i=0;i<len(f);++i)    f[i]*=x;
        return f;
    }
    poly divPoly(poly f,double x)
    {
        for(int i=0;i<len(f);++i)    f[i]/=x;
        return f;
    }
}using namespace Poly;
int main()
{
    int waste,x;
    scanf("%d",&waste);
    poly one,two,thr;
    while(waste--)
    {
        scanf("%d",&x);
        if(x>=len(one))    one.resize(x+1);
        if((x<<1)>=len(two))    two.resize(x<<1|1);
        if((x<<1)+x>=len(thr))    thr.resize((x<<1)+x+1);
        one[x]=comp(one[x].real()+1,0);
        two[x<<1]=comp(two[x<<1].real()+1,0);
        thr[(x<<1)+x]=comp(thr[(x<<1)+x].real()+1,0);
    }
    poly ans=addPoly(addPoly(one,divPoly(decPoly(mulPoly(one,one),two),2)),divPoly(addPoly(decPoly(mulPoly(mulPoly(one,one),one),mulPoly(mulPoly(one,two),3)),mulPoly(thr,2)),6));
    for(int i=0;i<len(ans);++i)    if(int(ans[i].real()+eps))    printf("%d %d\n",i,int(ans[i].real()+eps));
    return 0;
}

Descrption

Link.

对于每一个 $i$,求出:

$$ \sum_{j=1}^{i-1}\frac{a_{j}}{(i-j)^{2}}-\sum_{j=i+1}^{n}\frac{a_{j}}{(i-j)^{2}} $$

Solution

令 $f(i)=a_{i},g(i)=\frac{1}{i^{2}}$。

然后

$$ \sum_{j=1}^{i-1}f(j)\times g(i-j)-\sum_{j=i+1}^{n}f(j)\times g(j-i) $$

可以 FFT 了。

#include<bits/stdc++.h>
using namespace std;
const double bh_pi=acos(-1);
int n;
namespace Poly
{
    typedef complex<double> comp;
    typedef vector<complex<double> > poly;
    #define len(x) (int((x).size()))
    int lim,rev[400010];
    void fft(poly &f,int op)
    {
        for(int i=0;i<lim;++i)    if(i<rev[i])    swap(f[i],f[rev[i]]);
        for(int len=2;len<=lim;len<<=1)
        {
            comp bas(cos(2*bh_pi/len),op*sin(2*bh_pi/len));
            for(int fr=0;fr<lim;fr+=len)
            {
                comp now(1,0);
                for(int ba=fr;ba<fr+(len>>1);++ba,now*=bas)
                {
                    comp tmp=now*f[ba+(len>>1)];
                    f[ba+(len>>1)]=f[ba]-tmp;
                    f[ba]+=tmp;
                }
            }
        }
        if(op==-1)    for(int i=0;i<lim;++i)    f[i]/=lim;
    }
    poly mulPoly(poly f,poly g)
    {
        int n=len(f)+len(g)-1;
        for(lim=1;lim<=n;lim<<=1);
        for(int i=0;i<lim;++i)    rev[i]=(rev[i>>1]>>1)|((i&1)?(lim>>1):0);
        f.resize(lim),g.resize(lim);
        fft(f,1),fft(g,1);
        for(int i=0;i<lim;++i)    f[i]*=g[i];
        fft(f,-1),f.resize(n);
        return f;
    }
}using namespace Poly;
int main()
{
    poly f,g;
    scanf("%d",&n);
    f.resize(n+1),g.resize(n+1);
    for(int i=1;i<=n;++i)
    {
        double x;
        scanf("%lf",&x);
        f[i]=comp(x,0);
    }
    for(int i=1;i<=n;++i)    g[i]=comp(1.0/i/i,0);
    poly onetmp=mulPoly(f,g);
    reverse(next(f.begin()),f.end());
    poly anotmp=mulPoly(f,g);
    for(int i=1;i<=n;++i)    printf("%.3lf\n",onetmp[i].real()-anotmp[n-i+1].real());
    return 0;
}

Part. 1 Stirling Number / FK.

Def. 定义 $\begin{bmatrix}n \\ m\end{bmatrix}$ 表示将 $n$ 个元素分成 $m$ 个环的方案数。

递推式为

$$ \begin{bmatrix}n \\ m\end{bmatrix}=\begin{bmatrix}n-1 \\ m-1\end{bmatrix}+(n-1)\begin{bmatrix}n-1 \\ m\end{bmatrix} $$

即考虑已经放好的 $n-1$ 个数,第一种情况是自成环即 $\begin{bmatrix}n-1 \\ m-1\end{bmatrix}$,第二种情况是放在某一个环内,可以放在任意一个已经放好的数前,即 $(n-1)\begin{bmatrix}n-1 \\ m\end{bmatrix}$。

边界为:

$$ \begin{bmatrix}n \\ 0\end{bmatrix}=[n=0] $$

一个性质:

$$ n!=\sum_{i=0}^{n}\begin{bmatrix}n \\ i\end{bmatrix} $$

多项式形式:$\begin{bmatrix}n \\ m\end{bmatrix}$ 为 $f_{n}(x)=\prod_{i=0}^{n-1}(x+i)$ 的 $k$ 次项系数。

Part. 2 Stirling Number / SK.

Def. 定义 $\begin{Bmatrix}n \\ m\end{Bmatrix}$ 表示将 $n$ 个不同的球放进 $m$ 个相同的盒子的方案数。

递推式为

$$ \begin{Bmatrix}n \\ m\end{Bmatrix}=\begin{Bmatrix}n-1 \\ m-1\end{Bmatrix}+m\begin{Bmatrix}n-1 \\ m\end{Bmatrix} $$

意义即前 $n-1$ 个球放进了 $m-1$ 的盒子里,$n$ 就只有一种方法,如果前 $n-1$ 个球放进了 $m$ 个盒子,那么这个球就有 $m$ 种放法。

容斥一下可以得到通项公式(背吧)

$$ \begin{Bmatrix}n \\ m\end{Bmatrix}=\frac{1}{m!}\sum_{i=0}^{m}\left((-1)^{i}{m\choose i}(m-i)^{n}\right) $$

拆完组合数可以卷积 $\texttt{NTT}$ 做,不过咱多项式学得废,就不整了。

Description

Link.

$$ \sum\prod_{i=1}^{m}F_{a_{i}},(m>0,a_{1},\cdots a_{m}>0,\sum a_{i}=n) $$

Solution

这是一篇不用 $\mathbf{OGF}$ 的题解。

设 $f_{i}$ 为 $i$ 的 $\operatorname{lqp}$ 拆分值。

然后有显然的过不了递推式:

$$ f_{n}=\begin{cases} 1,n=0 \\ \displaystyle \sum_{i=0}^{n-1}F_{n-i}\times f_{i},n\neq0 \end{cases} $$

然后传统艺能错位相减操作一下:

$$ \begin{aligned} f_{n}&=\sum_{i=0}^{n-1}F_{n-i}\times f_{i} \\ f_{n-1}&=\sum_{i=0}^{n-2}F_{n-i-1}\times f_{i} \\ f_{n-2}&=\sum_{i=0}^{n-3}F_{n-i-2}\times f_{i} \end{aligned} \Longrightarrow \begin{aligned} f_{n}-f_{n-1}-f_{n-2}&=\sum_{i=0}^{n-1}F_{n-i}\times f_{i}-\sum_{i=0}^{n-2}F_{n-i-1}\times f_{i}-\sum_{i=0}^{n-3}F_{n-i-2}\times f_{i} \\ f_{n}-f_{n-1}-f_{n-2}&=(F_{2}-F_{1})\times f_{n-2}+F_{1}\times f_{n-1} \end{aligned} \\ \downarrow \\ f_{n}=2f_{n-1}+f_{n-2} $$

递推公式有了,然后矩阵快速幂:

$$ \begin{bmatrix} f_{n} \\ f_{n-1} \end{bmatrix} =\begin{bmatrix} 2f_{n-1}+f_{n-2} \\ f_{n-1} \end{bmatrix} =\begin{bmatrix} f_{n-1} \\ f_{n-2} \end{bmatrix} \times\begin{bmatrix} 2 & 1 \\ 1 & 0 \end{bmatrix} $$

这样就可以做了(吗?):

(code?)

#include <cstdio>
#include <iostream>
#include <cstring>
#include <queue>
#define mod ( 1000000007 )

using namespace std;
typedef long long LL;

template<typename _T, typename _P>
_T qkpow( _T bas, _T one, _P fur ){
    _T res = one;
    while( fur != 0 ){
        if( fur % 2 == ( _P )1 )    res = bas * res;
        bas = bas * bas;
        fur /= 2;
    }
    return res;
}

template<typename _T>
_T add( _T x, _T y ){ if( y >= mod )    y %= mod; x += y; if( x >= mod )    x -= mod; return x; }

struct bigInt : vector<int>{
    bigInt &check( ){
        while( ! empty( ) && ! back( ) ) pop_back( );
        if( empty( ) )    return *this;
        for( unsigned i = 1; i < size( ); ++ i ){ ( *this )[i] += ( *this )[i - 1] / 10; ( *this )[i - 1] %= 10; }
        while( back( ) >= 10 ){ push_back( back( ) / 10 ); ( *this )[size( ) - 2] %= 10; }
        return *this;
    }
    bigInt( int tpN = 0 ){ push_back( tpN ); check( ); }
};
istream &operator >> ( istream &is, bigInt &tpN ){
    string s;
    is >> s; tpN.clear( );
    for( int i = s.size( ) - 1; i >= 0; --i ) tpN.push_back( s[i] - '0' );
    return is;
}
ostream &operator << ( ostream &os, const bigInt &tpN ){
    if( tpN.empty( ) )    os << 0;
    for( int i = tpN.size( ) - 1; i >= 0; --i )    os << tpN[i];
    return os;
}
bool operator != ( const bigInt &one, const bigInt &another ){
    if( one.size( ) != another.size( ) )    return 1;
    for( int i = one.size( ) - 1; i >= 0; --i ){
        if( one[i] != another[i] )    return 1;
    }
    return 0;
}
bool operator == ( const bigInt &one, const bigInt &another ){
    return ! ( one != another );
}
bool operator < ( const bigInt &one, const bigInt &another ){
    if( one.size( ) != another.size( ) )    return one.size( ) < another.size( );
    for( int i = one.size( ) - 1; i >= 0; --i ){
        if( one[i] != another[i] )    return one[i] < another[i];
    }
    return 0;
}
bool operator > ( const bigInt &one, const bigInt &another ){ return another < one; }
bool operator <= ( const bigInt &one, const bigInt &another ){ return ! (one > another ); }
bool operator >= ( const bigInt &one, const bigInt &another ){ return ! (one < another ); }
bigInt &operator += ( bigInt &one, const bigInt &another ){
    if( one.size( ) < another.size( ) )    one.resize(another.size( ) );
    for( unsigned i = 0; i != another.size( ); ++ i ) one[i] += another[i];
    return one.check( );
}
bigInt operator + ( bigInt one, const bigInt &another ){ return one += another; }
bigInt &operator -= ( bigInt &one, bigInt another ){
    if( one < another )    swap( one, another );
    for( unsigned i = 0; i != another.size( ); one[i] -= another[i], ++ i ){
        if( one[i] < another[i] ){
            unsigned j = i + 1;
            while( ! one[j] ) ++ j;
            while( j > i ){ -- one[j]; one[--j] += 10; }
        }
    }
    return one.check( );
}
bigInt operator - ( bigInt one, const bigInt &another ){ return one -= another; }
bigInt operator * ( const bigInt &one, const bigInt &another ){
    bigInt tpN;
    tpN.assign( one.size( ) + another.size( ) - 1, 0 );
    for( unsigned i = 0; i != one.size( ); ++ i ){
        for( unsigned j = 0; j != another.size( ); ++ j ) tpN[i + j] += one[i] * another[j];
    }
    return tpN.check( );
}
bigInt &operator *= ( bigInt &one, const bigInt &another ){ return one = one * another; }
bigInt divMod( bigInt &one, const bigInt &another ){
    bigInt ans;
    for( int t = one.size( ) - another.size( ); one >= another; -- t ){
        bigInt tpS;
        tpS.assign( t + 1, 0 );
        tpS.back( ) = 1;
        bigInt tpM = another * tpS;
        while( one >= tpM ){ one -= tpM; ans += tpS; }
    }
    return ans;
}
bigInt operator / ( bigInt one, const bigInt &another ){ return divMod(one, another ); }
bigInt &operator /= ( bigInt &one, const bigInt &another ){ return one = one / another; }
bigInt &operator %= ( bigInt &one, const bigInt &another ){ divMod( one, another ); return one; }
bigInt operator % ( bigInt one, const bigInt &another ){ return one %= another; }

struct matrixS{
    int mat[2][2];
    matrixS( int x = 0 ){ memset( mat, x, sizeof( mat ) ); }
    matrixS operator * ( const matrixS &another ) const{
        matrixS res;
        for( int i = 0; i < 2; ++ i ){
            for( int j = 0; j < 2; ++ j ){
                for( int k = 0; k < 2; ++ k )    res.mat[i][j] = add( ( LL )res.mat[i][j], ( LL )mat[i][k] * another.mat[k][j] );
            }
        }
        return res;
    }
} unit, erng;

bigInt N;

void progressBaseInformation( ){
    int unitS[2][2] = { { 1, 0 }, { 0, 1 } };
    memcpy( unit.mat, unitS, sizeof( unitS ) );
    int erngS[2][2] = { { 2, 1 }, { 1, 0 } };
    memcpy( erng.mat, erngS, sizeof( erngS ) );
}

signed main( ){
    ios::sync_with_stdio( 0 ); cin.tie( 0 ); cout.tie( 0 );
    progressBaseInformation( );
    cin >> N; cout << qkpow( erng, unit, N ).mat[1][0] << '\n';
    return 0;
}

不,凉心出题人友好地卡了高精的常数,于是你打开题解,发现 $f_{n}=f_{n\bmod (10^{9}+6)}$,于是你又行了。

$\mathcal{Code}$

#include <cstdio>
#include <cstring>
#include <queue>
#define mod ( 1000000007 )

using namespace std;
typedef long long LL;

template<typename _T>
void read( _T &x ){
    x = 0; char c = getchar( ); _T f = 1;
    while( c < '0' || c > '9' ){ if( c == '-' )    f = -1; c = getchar( ); }
    while( c >= '0' && c <= '9' ){ x = ( ( x << 3 ) + ( x << 1 ) + ( c & 15 ) ) % ( mod - 1 ); c = getchar( ); }
    x *= f;
}

template<typename _T>
void write( _T x ){
    if( x < 0 ){ putchar( '-' ); x = -x; }
    if( x > 9 )    write( x / 10 );
    putchar( x % 10 + '0' );
}

template<typename _T, typename _P>
_T qkpow( _T bas, _T one, _P fur ){
    _T res = one;
    while( fur != 0 ){
        if( fur % 2 == ( _P )1 )    res = bas * res;
        bas = bas * bas;
        fur /= 2;
    }
    return res;
}

template<typename _T>
_T add( _T x, _T y ){ if( y >= mod )    y %= mod; x += y; if( x >= mod )    x -= mod; return x; }

struct matrixS{
    int mat[2][2];
    matrixS( int x = 0 ){ memset( mat, x, sizeof( mat ) ); }
    matrixS operator * ( const matrixS &another ) const{
        matrixS res;
        for( int i = 0; i < 2; ++ i ){
            for( int j = 0; j < 2; ++ j ){
                for( int k = 0; k < 2; ++ k )    res.mat[i][j] = add( ( LL )res.mat[i][j], ( LL )mat[i][k] * another.mat[k][j] );
            }
        }
        return res;
    }
} unit, erng;

LL N;

void progressBaseInformation( ){
    int unitS[2][2] = { { 1, 0 }, { 0, 1 } };
    memcpy( unit.mat, unitS, sizeof( unitS ) );
    int erngS[2][2] = { { 2, 1 }, { 1, 0 } };
    memcpy( erng.mat, erngS, sizeof( erngS ) );
}

signed main( ){
    progressBaseInformation( );
    read( N ); write( qkpow( erng, unit, N ).mat[1][0] ), putchar( '\n' );
    return 0;
}

$\mathbf{OGF}$ 的定义

对于一个序列 $a_{1},a_{2},\cdots$,我们称:

$$ G(x)=\sum_{i=0}^{\infty}a_{i}x^{i} $$

为序列 $a$ 的 $\mathbf{OGF}$ 即普通生成函数 ($\texttt{Ordinary Generating Function}$)。

同时因为我们不关心 $x$ 的取值,因此 $\sum_{i=1}^{+\infty}a_{i}x^{i}$ 又称作以 $x$ 为自由元的形式幂级数。 -- 摘自 自为风月马前卒

一个例子

举个例子,序列:

$$ \left(^{n}_{0}\right),\left(^{n}_{1}\right),\cdots,\left(^{n}_{n}\right) $$

的 $\mathbf{OGF}$ 为(二项式定理):

$$ G(x)=(1+x)^{n} $$

由等比数列求和公式,有一个常用的等式:

$$ \sum_{i=0}^{\infty}x^{i}=\frac{1-x^{\infty}}{1-x} $$

因为指数为 $\infty$,所以 $x^{\infty}$ 趋近于 $0$,箭头方向随便打,因为我们并不关心 $x$ 的取值。

$$ \sum_{i=0}^{\infty}x^{i}=\frac{1}{1-x} $$

这个等式还有一个重要的运用,我们把 $x$ 替换成 $kx$ 即可得:

$$ \sum_{i=0}^{\infty}(kx)^{i}=\frac{1}{1-kx} $$

后文的用 $\mathbf{OGF}$ 求序列的通项公式里面这个东西很有用的。

$\texttt{Fibonacci}$ 序列的生成函数求法

定义一个序列

$$ F_{i}=\begin{cases} 1,i\in[0,1] \\ \displaystyle F_{i-1}+F_{i-2},i\in[2,\infty) \end{cases} $$

则我们称 $F$ 为 $\texttt{Fibonacci}$ 序列。

接下来我们来推导其生成函数:

$$ \begin{aligned} G(x)&=\sum_{i=0}^{\infty}F_{i}x^{i} \\ G(x)&=1+x+2x^{2}+3x^{3}+\cdots \\ xG(x)&=x+x^{2}+2x^{3}+3x^{4}+\cdots \\ x^{2}G(x)&=x^{2}+x^{3}+2x^{4}+3x^{5}+\cdots \end{aligned} $$

这里运用初中数学中经常用的到错位相减这一小技巧,可得

$$ G(x)-xG(x)-x^{2}G(x)=1 $$

即可得

$$ G(x)=\frac{1}{1-x-x^{2}} $$

至此,我们已经求出了 $\texttt{Fibonacci}$ 序列的 $\mathbf{OGF}$ 了。

利用生成函数求数列通项

以前文提到的 $\texttt{Fibonacci}$ 为例。

首先我们知道其 $\mathbf{OGF}$ 为:

$$ G(x)=\frac{1}{1-x-x^{2}} $$

待定系数一下分母我们就可以得到:

$$ G(x)=\frac{1}{(1-\frac{1+\sqrt{5}}{2}x)(1-\frac{1-\sqrt{5}}{2}x)} $$

~后面的还没推出来,咕了~