2020牛客暑期多校训练营（第六场）Interesting Stiriling【第一类斯特林数】【lucas】【反演】-CSDN博客

本文链接：https://blog.csdn.net/Zayin___/article/details/108242401

题目大意

求第一类斯特林数的区间和，即：
$\Bigg(\sum_{i=l}^{r}\begin{bmatrix}n\\i\end{bmatrix} \Bigg)\% p$

$0\leq l \leq r \leq n \leq 10^{18}$ ， $p\leq 10^6$

题解

~~首先把问题转化成前缀和这就不用说了吧~~
总所周知，第一类斯特林数的生成函数
$\sum_{i=0}^n\begin{bmatrix}n\\i\end{bmatrix}x^i=x^{\overline{n}}=\prod_{i=0}^{n-1}(x+i)$

那么在 $\%p$ 意义下可以写成 $\prod_{i=0}^{n-1}(x+i)=\prod_{i=0}^{n-1}(x+i\%p)$

令 $n=d\cdot p+r$ ，则有

$\prod_{i=0}^{n-1}(x+i\%p)=\prod_{k=0}^{d-1}\prod_{i=0}^{p-1}(x+i)\cdot\prod_{i=0}^{r-1}(x+i)=\bigg[\prod_{i=0}^{p-1}(x+i)\bigg]^d\cdot\prod_{i=0}^{r-1}(x+i)$

而 $\prod_{i=0}^{p-1}(x+i)$ 在 $\%p$ 意义下是等于 $x^p-x$ 的（考虑两个多项式的根都是 $0,1,\cdots,p-1$ ，并且最高项系数也相同），即：

$\sum_{i=0}^n\begin{bmatrix}n\\i\end{bmatrix}x^i=x^d\bigg[x^{p-1}-1\bigg]^d\cdot\prod_{i=0}^{r-1}(x+i)$

则左边的第 $k$ 项，则对应右边的 $d + k^{'} (p - 1) + i$ 的指数，其中 $i=(k-d)\%(p-1),k'=\lfloor\frac{k-d}{p-1}\rfloor$ ，即：
$\begin{bmatrix}n\\k\end{bmatrix}=(-1)^{d-k'}\binom{d}{k'} \cdot \begin{bmatrix}r\\i\end{bmatrix}$
那么则有：
$\begin{aligned} \sum_{k=0}^{m}\begin{bmatrix}n\\k\end{bmatrix}&=\sum_{k=0}^{m}(-1)^{d-k'}\binom{d}{k'} \cdot \begin{bmatrix}r\\i\end{bmatrix}\\ &=\sum_{i=0}^{p-2}\begin{bmatrix}r\\i\end{bmatrix}\sum_{k'=0}^{\lfloor\frac{m-d-i}{p-1}\rfloor}(-1)^{d-k'}\binom{d}{k'} \\ &=(-1)^d\sum_{i=0}^{p-2}\begin{bmatrix}r\\i\end{bmatrix}\sum_{k'=0}^{\lfloor\frac{m-d-i}{p-1}\rfloor}(-1)^{k'}\binom{d}{k'} \end{aligned}$

不难看出， $\lfloor\frac{m-d-i}{p-1}\rfloor$ 其实只有两种取值，分别对应 $i\in\Big[0,\cdots,(m-d)\%(p-1)\Big]$ 以及 $i\in\Big((m-d)\%(p-1),p-2\Big]$ ，分开这两段来求的话，现在的问题就变成了两个子问题：
$\begin{aligned} 1、&\sum_{i=0}^{m}\begin{bmatrix}n\\i\end{bmatrix}&,m\leq n<p-1 \\ 2、&\sum_{i=0}^{m}(-1)^i\binom{n}{i} &,m\leq n<10^{18} \end{aligned}$

先说第二个， $l u c a s$ 定理告诉我们：
$\binom{n}{m}\%p=\binom{n/p}{m/p}\cdot\binom{n\%p}{m\%p}\%p$
所以在模意义下有：
$\begin{aligned} \sum_{i=0}^{m}(-1)^i\binom{n}{i}&=\sum_{i=0}^{m}(-1)^i\binom{n/p}{i/p}\cdot\binom{n\%p}{i\%p}\\ &=\sum_{k=0}^{\lfloor\frac{m}{p}\rfloor-1}\binom{n/p}{k}\sum_{i=0}^{p-1}(-1)^i\cdot\binom{n\%p}{i}+\binom{n/p}{m/p}\sum_{i=0}^{m\%p}(-1)^i\cdot\binom{n\%p}{i} \end{aligned}$

注意 $n\%p\leq p-1$ ，所以：
$\begin{aligned} &\sum_{k=0}^{\lfloor\frac{m}{p}\rfloor-1}\binom{n/p}{k}\sum_{i=0}^{p-1}(-1)^i\cdot\binom{n\%p}{i}+\binom{n/p}{m/p}\sum_{i=0}^{m\%p}(-1)^i\cdot\binom{n\%p}{i}\\ =&\sum_{k=0}^{\lfloor\frac{m}{p}\rfloor-1}\binom{n/p}{k}\sum_{i=0}^{n\%p}(-1)^i\cdot\binom{n\%p}{i}+\binom{n/p}{m/p}\sum_{i=0}^{m\%p}(-1)^i\cdot\binom{n\%p}{i}\\ =&\sum_{k=0}^{\lfloor\frac{m}{p}\rfloor-1}\binom{n/p}{k}(1+(-1))^{n\%p}+\binom{n/p}{m/p}\sum_{i=0}^{m\%p}(-1)^i\cdot\binom{n\%p}{i}\\ =&\binom{n/p}{m/p}\sum_{i=0}^{m\%p}(-1)^i\cdot\binom{n\%p}{i}\\ \end{aligned}$

这样子就可以 $O (p)$ 算了（题解使用了一种更为普适的方法，有兴趣的读者可以去看题解）

接下来说怎么算 $\sum_{i=0}^{m}\begin{bmatrix}n\\i\end{bmatrix},m\leq n<p-1$ ，虽然n的范围已经降到1e6，但是第一类斯特林数并没有特别显式的公式，最多只有 $\sum_{i=0}^n\begin{bmatrix}n\\i\end{bmatrix}x^i=x^{\overline{n}}=\prod_{i=0}^{n-1}(x+i)$ ，常数小的选手可以尝试暴力分治FFT，标程就是这样 $O (n l o g n)$ 冲过去的，~~但我是大常数选手被卡也很正常~~，不过还有另外一种更为数学的方法——反演。

根据第一类斯特林的生成函数，可以写出下面的矩阵（为了排版好看，就用 $S_1(n,i)$ 来代替 $\begin{bmatrix}n\\i\end{bmatrix}$ ）

$\begin{bmatrix} {x_0}^0 & {x_0}^1 & {x_0}^2 & \cdots & {x_0}^{p-1}\\ {x_1}^0 & {x_1}^1 & {x_1}^2 & \cdots & {x_1}^{p-1}\\ {x_2}^0 & {x_2}^1 & {x_2}^2 & \cdots & {x_2}^{p-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ {x_{p-1}}^0 & {x_{p-1}}^1 & {x_{p-1}}^2 & \cdots & {x_{p-1}}^{p-1}\\ \end{bmatrix}\cdot \begin{bmatrix} S_1(n,0)\\ S_1(n,1)\\ S_1(n,2)\\ \vdots\\ S_1(n,p-1) \end{bmatrix}= \begin{bmatrix} {x_0}^{\overline{n}} \\ {x_1}^{\overline{n}}\\ {x_2}^{\overline{n}}\\ \vdots\\ {x_{p-1}}^{\overline{n}} \end{bmatrix}$

这就是用斯特林数算上升幂的直接应用嘛。
又因为 $S_1(n,0)=0$ ，可以直接去掉对应矩阵的列：
$\begin{bmatrix} {x_1}^1 & {x_1}^2 & {x_1}^3 & \cdots & {x_1}^{p-1}\\ {x_2}^1 & {x_2}^2 & {x_2}^3 & \cdots & {x_2}^{p-1}\\ {x_3}^1 & {x_3}^2 & {x_3}^3 & \cdots & {x_3}^{p-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ {x_{p-1}}^1 & {x_{p-1}}^2 & {x_{p-1}}^3 & \cdots & {x_{p-1}}^{p-1}\\ \end{bmatrix}\cdot \begin{bmatrix} S_1(n,1)\\ S_1(n,2)\\ S_1(n,3)\\ \vdots\\ S_1(n,p-1) \end{bmatrix}= \begin{bmatrix} {x_1}^{\overline{n}}\\ {x_2}^{\overline{n}}\\ {x_3}^{\overline{n}}\\ \vdots\\ {x_{p-1}}^{\overline{n}} \end{bmatrix}$

~~那么我们随便找几个数，O(n^3)高斯消元就能解出S1(n,i)了~~
看这个矩阵的形式，是不是感觉很眼熟，有点像 $F F T$ 的转移矩阵，那我们干脆一不做二不休，直接选择 $x_i=g^i$ 好了：

$\begin{bmatrix} g^1 & g^2 & g^3 & \cdots & g^{(p-1)}\\ g^2 & g^4 & g^6 & \cdots & g^{2(p-1)}\\ g^3 & g^6 & g^9 & \cdots & g^{3(p-1)}\\ \vdots & \vdots & \ddots & \vdots\\ g^{(p-1)} & g^{2(p-1)} & g^{3(p-1)} &\cdots & g^{(p-1)(p-1)} \end{bmatrix}\cdot \begin{bmatrix} S_1(n,1)\\ S_1(n,2)\\ S_1(n,3)\\ \vdots\\ S_1(n,(p-1)) \end{bmatrix}= \begin{bmatrix} g^{\overline{n}} \\ ({g^2})^{\overline{n}}\\ (g^3)^{\overline{n}}\\ \vdots\\ (g^{(p-1)})^{\overline{n}} \end{bmatrix}$

简写成：

$G S = V$

根据 $F F T$ 我们知道，如果 $g$ 是 $p$ 的原根的话， $G$ 是有逆矩阵的，即：
$S=G^{-1}V$

具体就是：

$\frac{1}{p-1} \begin{bmatrix} g^{-1} & g^{-2} & g^{-3} & \cdots & g^{-(p-1)}\\ g^{-2} & g^{-4} & g^{-6} & \cdots & g^{-2(p-1)}\\ g^{-3} & g^{-6} & g^{-9} & \cdots & g^{-3(p-1)}\\ \vdots & \vdots & \ddots & \vdots\\ g^{-(p-1)} & g^{-2(p-1)} & g^{-3(p-1)} &\cdots & g^{-(p-1)(p-1)} \end{bmatrix}\cdot \begin{bmatrix} g^{\overline{n}} \\ ({g^2})^{\overline{n}}\\ (g^3)^{\overline{n}}\\ \vdots\\ (g^{(p-1)})^{\overline{n}} \end{bmatrix}= \begin{bmatrix} S_1(n,1)\\ S_1(n,2)\\ S_1(n,3)\\ \vdots\\ S_1(n,(p-1)) \end{bmatrix}$

即：
$\begin{aligned} S_1(n,i)&=\sum_{j=1}^{p-1}g^{-ij}\cdot(g^j)^{\overline{n}}\\ \end{aligned}$

接着求前缀和：
$\begin{aligned} \sum_{i=1}^{m}S_1(n,i)&=\sum_{i=1}^{m}\sum_{j=1}^{p-1}g^{-ij}\cdot(g^j)^{\overline{n}}\\ &=\sum_{j=1}^{p-1}(g^j)^{\overline{n}}\sum_{i=1}^{m}g^{-ij}\\ &=\sum_{j=1}^{p-1}(g^j)^{\overline{n}}\cdot\Big[g^j=1?m:\frac{g^{-j(m+1)}-g^{-j}}{g^{-j}-1}\Big] \end{aligned}$

上升幂可以用阶乘来算，所以这样的复杂度是 $O (p)$ ，总复杂度也为 $O (p)$ ，~~比题解不知道要高多少倍~~。

#include<bits/stdc++.h>
#define maxn 2000050
using namespace std;
typedef long long LL;

LL n,l,r,modu;

void add(LL &a,LL b)    {   a=(a+b)%modu;   }
void sub(LL &a,LL b)    {   a=(a-b%modu+modu)%modu;  }
void mul(LL &a,LL b)    {   a=a*b%modu; }

LL pw(LL a,LL k=modu-2)   {
    LL ans=1;
    for (a%=modu;k;k>>=1)  {
        if (k&1)
            ans=ans*a%modu;
        a=a*a%modu;
    }
    return ans;
}

LL F[maxn],I[maxn];

LL sgn(LL x)    {   
    return x&1?modu-1:1;
}

LL C(LL n,LL m) {
    if (!n&&!m) return 1;
    LL _n=n%modu,_m=m%modu;
    LL ans=_n<_m?0:F[_n]*I[_m]%modu*I[_n-_m]%modu;
    return ans*C(n/modu,m/modu);
}

int g;

vector<int> factor(LL x)    {
    vector<int> F;
    for (int i=2;i<=x;++i)  {
        if (x%i) continue;
        while (x%i==0) x/=i;
        F.push_back(i);
    }
    return F;
}

void predo()    {
    for (int i=F[0]=I[0]=I[1]=1;i<modu;++i) F[i]=F[i-1]*i%modu;
    for (int i=2;i<modu;++i)    I[i]=modu-modu/i*I[modu%i]%modu;
    for (int i=1;i<modu;++i)    I[i]=I[i-1]*I[i]%modu;

    auto F=factor(modu-1);
    for (g=2;;++g)  {
        bool flag=1;
        for (LL x:F)
            flag&=(pw(g,(modu-1)/x)!=1);
        if (flag) break;
    }
}

LL sum_Stiriling(LL n,LL m) {
    m=min(n,m);
    LL ans=0;
    for (int j=1;j<modu;++j)    {
        auto sum=[&](LL x,LL m)    {
            if (x==1)
                return m%modu;
            else    
                return (pw(x,m+1)-x+modu)*pw(x-1+modu)%modu;
        };
        auto upow=[&](LL x,LL n)    {
            return F[x+n-1]*I[x-1]%modu;
        };
        LL x=pw(g,j);
        add(ans,upow(x,n)*sum(pw(x),m));
    }
    return (modu-ans)%modu;
}

LL sum_Binomial(LL n,LL m) {//with (-1)^i
    if (m<0) return 0;
    LL ans=0;
    for (int i=0;i<=m%modu;++i)
        add(ans,C(n%modu,i)*sgn(i));
    return ans*C(n/modu,m/modu)%modu*sgn(m/modu);
}

LL solve(LL n,LL m)   {
    LL _n=n%modu;
    if (m<=n/modu) return 0;
    LL _m=m-n/modu,x=_m/(modu-1);

    LL ans=sum_Binomial(n/modu,x-1)*F[_n]+sum_Stiriling(_n,_m%(modu-1))*C(n/modu,x)%modu*sgn(x);
    ans%=modu;
    return ans*sgn(n/modu)%modu;
}

int main()  {
    scanf("%lld%lld%lld%lld",&n,&l,&r,&modu);
    predo();
    printf("%lld\n",(solve(n,r)-solve(n,l-1)+modu)%modu);
    return 0;
}