Stirling Number for query
(Math/stirlingquery.hpp)

• View this file on GitHub
• Last update: 2023-01-17 01:58:11+09:00
• Include: #include "Math/stirlingquery.hpp"

Depends on

• Fast Division (Math/fastdiv.hpp)

Verified with

• Verify/LC_stirling_number_of_the_first_kind_small_p_large_n.test.cpp
• Verify/LC_stirling_number_of_the_second_kind_small_p_large_n.test.cpp

Code

#pragma once
#include "Math/fastdiv.hpp"


class StirlingNumberQuery{
    const int p;
    FastDiv ip;
    vector<vector<int>> binom,F,S;
    ll nCr(ll n,ll k){
        if(n<0 or k<0 or n<k)return 0;
        ll res=1;
        while(n){
            res=(res*binom[n%ip][k%ip])%ip;
            n/=p; k/=p;
        }
        return res;
    }
public:
    StirlingNumberQuery(int _p):p(_p),ip(p){
        binom.resize(p,vector<int>(p));
        F.resize(p,vector<int>(p));
        S.resize(p,vector<int>(p));
        binom[0][0]=F[0][0]=S[0][0]=1;
        rep(n,1,p)rep(k,0,n+1){
            if(k)binom[n][k]=binom[n-1][k-1];
            binom[n][k]=(binom[n][k]+binom[n-1][k])%ip;

            if(k)F[n][k]=F[n-1][k-1];
            F[n][k]=(F[n][k]+ll(p-n+1)*F[n-1][k])%ip;
            
            if(k)S[n][k]=S[n-1][k-1];
            S[n][k]=(S[n][k]+ll(k)*S[n-1][k])%ip;
        }
    }
    int FirstKind(ll n,ll k){
        if(n<0 or k<0 or k>n)return 0;
        ll i=n/ip,j=n%ip;
        if(k<i)return 0;
        ll a=(k-i)/(p-1),b=(k-i)%(p-1);
        if(b==0 and j)b+=p-1,a--;
        if(a<0 or a>i or b>j)return 0;
        int res=(nCr(i,a)*F[j][b])%ip;
        if((i+a)&1)res=(p-res)%ip;
        return res;
    }
    int SecondKind(ll n,ll k){
        if(n<0 or k<0 or k>n)return 0;
        if(n==0)return 1;
        ll i=k/ip,j=k%ip;
        if(n<i)return 0;
        ll a=(n-i)/(p-1),b=(n-i)%(p-1);
        if(b==0)b+=p-1,a--;
        if(a<0 or b<j)return 0;
        if(b==p-1 and j==0)return nCr(a,i-1);
        else return (nCr(a,i)*S[b][j])%ip;
    }
};

/**
 * @brief Stirling Number for query
 */

#line 2 "Math/fastdiv.hpp"

struct FastDiv{
    using u64=uint64_t;
    using u128=__uint128_t;
    constexpr FastDiv():m(),s(),x(){}
    constexpr FastDiv(int _m)
        :m(_m),s(__lg(m-1)),x(((u128(1)<<(s+64))+m-1)/m){}
    constexpr int get(){return m;}
    constexpr friend u64 operator/(u64 n,const FastDiv& d){
        return (u128(n)*d.x>>d.s)>>64;
    }
    constexpr friend int operator%(u64 n,const FastDiv& d){
        return n-n/d*d.m;
    }
    constexpr pair<u64,int> divmod(u64 n)const{
        u64 q=n/(*this);
        return {q,n-q*m};
    }
    int m,s; u64 x;
};

/**
 * @brief Fast Division
*/
#line 3 "Math/stirlingquery.hpp"

class StirlingNumberQuery{
    const int p;
    FastDiv ip;
    vector<vector<int>> binom,F,S;
    ll nCr(ll n,ll k){
        if(n<0 or k<0 or n<k)return 0;
        ll res=1;
        while(n){
            res=(res*binom[n%ip][k%ip])%ip;
            n/=p; k/=p;
        }
        return res;
    }
public:
    StirlingNumberQuery(int _p):p(_p),ip(p){
        binom.resize(p,vector<int>(p));
        F.resize(p,vector<int>(p));
        S.resize(p,vector<int>(p));
        binom[0][0]=F[0][0]=S[0][0]=1;
        rep(n,1,p)rep(k,0,n+1){
            if(k)binom[n][k]=binom[n-1][k-1];
            binom[n][k]=(binom[n][k]+binom[n-1][k])%ip;

            if(k)F[n][k]=F[n-1][k-1];
            F[n][k]=(F[n][k]+ll(p-n+1)*F[n-1][k])%ip;
            
            if(k)S[n][k]=S[n-1][k-1];
            S[n][k]=(S[n][k]+ll(k)*S[n-1][k])%ip;
        }
    }
    int FirstKind(ll n,ll k){
        if(n<0 or k<0 or k>n)return 0;
        ll i=n/ip,j=n%ip;
        if(k<i)return 0;
        ll a=(k-i)/(p-1),b=(k-i)%(p-1);
        if(b==0 and j)b+=p-1,a--;
        if(a<0 or a>i or b>j)return 0;
        int res=(nCr(i,a)*F[j][b])%ip;
        if((i+a)&1)res=(p-res)%ip;
        return res;
    }
    int SecondKind(ll n,ll k){
        if(n<0 or k<0 or k>n)return 0;
        if(n==0)return 1;
        ll i=k/ip,j=k%ip;
        if(n<i)return 0;
        ll a=(n-i)/(p-1),b=(n-i)%(p-1);
        if(b==0)b+=p-1,a--;
        if(a<0 or b<j)return 0;
        if(b==p-1 and j==0)return nCr(a,i-1);
        else return (nCr(a,i)*S[b][j])%ip;
    }
};

/**
 * @brief Stirling Number for query
 */

Back to top page