Language: C++

#include<bits/stdc++.h>
using namespace std;
using namespace ProconLib;

namespace ProconLib{
    
    template<typename T>
    struct Detail{
        static constexpr T EPS = T(1e-9);
        static bool isZero(T v){
            return abs(v)<=EPS;
        }
    };

    template<typename T>
    class Vector{
        int N;
        std::vector<T> dat;
        public:
        Vector(int n):N(n),dat(n){}
        Vector(int n,T x):N(n),dat(n,x){}
        Vector(std::vector<T> vec):N(vec.size()),dat(vec){}
        Vector(const std::vector<T>& vec):N(vec.size()),dat(vec){}
        Vector(const Vector& vec){dat=vec.dat;}
        int size(){return N;}
        T& operator[](int pos){return dat[pos];};
        const T& operator[](int pos) const{return dat[pos];};
        T& at(int pos){return dat.at(pos);}
        const T& at(int pos) const {return dat.at(pos);}
    };
    
    template<typename T>
    class Matrix{
        int r,c;
        std::vector<Vector<T>> dat;
        public:
        Matrix(int r,int c):r(r),c(c),dat(r,Vector<T>(c)){}
        Matrix(int r,int c,T x):r(r),c(c),dat(r,Vector<T>(c,x)){}

        Vector<T>& operator[](int pos){return dat[pos];};
        const Vector<T>& operator[](int pos) const{return dat[pos];};
        Vector<T>& at(int pos){return dat.at(pos);}
        const Vector<T>& at(int pos) const {return dat.at(pos);}
        int rowSize() const {return r;}
        int colSize() const {return c;}
    };
    template<typename T>
    Matrix<T> operator+(const Matrix<T>& lhs,const Matrix<T>& rhs);
    template<typename T>
    Matrix<T> operator-(const Matrix<T>& lhs,const Matrix<T>& rhs);
    template<typename T>
    Matrix<T> operator+(const Matrix<T>& mat);
    template<typename T>
    Matrix<T> operator-(const Matrix<T>& mat);
    template<typename T>
    Matrix<T> operator*(const Matrix<T>& lhs,const Matrix<T> &rhs);
    template<typename T>
    Matrix<T> Identity(int n);
    template<typename T>
    Matrix<T> pow(const Matrix<T>& mat,int k);
    template<typename T,typename Detail=Detail<T>>
    Matrix<T> gaussianElimination(Matrix<T> mat);
    template<typename T>
    int rank(const Matrix<T>& mat);
    template<typename T>
    Matrix<T> inv(const Matrix<T>& mat);
    
    template<typename T>
    Matrix<T> operator+(const Matrix<T>& lhs,const Matrix<T>& rhs){
        assert(lhs.rowSize()==rhs.rowSize() && lhs.colSize()==rhs.colSize());
        int r=lhs.rowSize(),c=lhs.colSize();
        Matrix<T> res(r,c);
        for(int i=0;i<r;i++){
            for(int j=0;j<c;j++){
                res[i][j]=lhs[i]+rhs[i];
            }
        }
        return res;
    }
    template<typename T>
    Matrix<T> operator-(const Matrix<T>& lhs,const Matrix<T>& rhs){
        assert(lhs.rowSize()==rhs.rowSize() && lhs.colSize()==rhs.colSize());
        int r=lhs.rowSize(),c=lhs.colSize();
        Matrix<T> res(r,c);
        for(int i=0;i<r;i++){
            for(int j=0;j<c;j++){
                res[i][j]=lhs[i]-rhs[i];
            }
        }
        return res;
    }
    template<typename T>
    Matrix<T> operator+(const Matrix<T>& mat){
        int r=mat.rowSize(),c=mat.colSize();
        Matrix<T> res(r,c);
        for(int i=0;i<r;i++)for(int j=0;j<c;j++) res[i][j]=-mat[i][j];
        return res;
    }
    template<typename T>
    Matrix<T> operator-(const Matrix<T>& mat){
        int r=mat.rowSize(),c=mat.colSize();
        Matrix<T> res(r,c);
        for(int i=0;i<r;i++)for(int j=0;j<c;j++) res[i][j]=-mat[i][j];
        return res;
    }
    
    template<typename T>
    Matrix<T> operator*(const Matrix<T>& lhs,const Matrix<T> &rhs){
        assert(lhs.colSize()==rhs.rowSize());
        int r=lhs.rowSize(),c=rhs.colSize(),l=lhs.colSize();
        Matrix<T> res(r,c);
        for(int i=0;i<r;i++){
            for(int k=0;k<l;k++){
                for(int j=0;j<c;j++){
                    res[i][j]+=lhs[i][k]*rhs[k][j];
                }
            }
        }
        return res;
    }

    template<typename T>
    Matrix<T> Identity(int n){
        assert(n>=0);
        Matrix<T> res(n,n);
        for(int i=0;i<n;i++){
            res[i][i]=1;
        }
        return res;
    }

    template<typename T>
    Matrix<T> pow(const Matrix<T>& mat,int k){
        assert(mat.rowSize()==mat.colSize());
        Matrix<T> x=mat;
        Matrix<T> res=Identity<T>(mat.rowSize());
        while(k){
            if(k&1) res=res*x;
            x=x*x;
            k>>=1;
        }
        return res;
    }
    
    template<typename T,typename Detail=Detail<T>>
    Matrix<T> gaussianElimination(Matrix<T> mat){
        int m=mat.rowSize(),n=mat.colSize();
        int row=0;
        for(int j=0;j<n;j++){
            if(row==m) break;
            int tar=-1;
            T v=Detail::EPS;
            for(int i=row;i<m;i++){
                if(!Detail::isZero(mat[i][j]) && abs(v)<abs(mat[i][j])){
                    tar=i;
                    v=mat[i][j];
                }
            }
            if(tar==-1) continue;
            if(row!=tar){ 
                auto tmp=mat[row];
                mat[row]=mat[tar];
                mat[tar]=tmp;
            }
            for(int i=row+1;i<m;i++){
                if(!Detail::isZero(mat[i][j])){
                    T r=mat[i][j]/mat[row][j];
                    for(int k=j;k<n;k++){
                        mat[i][k]-=r*mat[row][k];
                    }
                }
            }
            row++;
        }
        return mat;
    }

    template<typename T,typename Detail=Detail<T>>
    int rank(const Matrix<T>& mat){
        auto tmp=gaussialElimination(mat);
        int m=tmp.rowSize(),n=tmp.colSize();
        int i=0,j=0;
        while(i<m && j<n){
            if(Detail::isZero(tmp[i][j])) j++;
            else i++,j++;
        }
        return i;
    }

    template<typename T>
    Matrix<T> inv(const Matrix<T>& mat){
        assert(mat.rowSize()==mat.colSize());
        int n=mat.rowSize();
        Matrix<T> tmp(n,2*n);
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                tmp[i][j]=mat[i][j];
                tmp[i][j+n]=0;
            }
            tmp[i][i+n]=1;
        }
        mat=gaussianElimination(mat);
        Matrix<T> res(n);
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                res[i][j]=mat[i][j+n]/=mat[i][i];
            }
        }
        return res;
    }
 
    template<typename T>
    void debug(Matrix<T> mat){
        int m=mat.rowSize(),n=mat.colSize();
        std::cerr<<"###Matrix_Debug###"<<std::endl;
        for(int i=0;i<m;i++){
            for(int j=0;j<n;j++){
                std::cerr<<mat[i][j]<<" ";
            }
            std::cerr<<std::endl;
        }
    }
}

int main(){
    int s,n,k;
    cin>>s>>n>>k;
    if(n==1 && s%k!=0){
        cout<<-1<<endl;
        return 0;
    }
    int sz=n*k+1;
    
    vector<vector<double>> dp(k+1,vector<double>(sz));
    dp[0][0]=1;
    for(int i=0;i<k;i++){
        for(int j=0;j<sz;j++){
            for(int k=1;k<=n;k++){
                if(j+k<sz) dp[i+1][j+k]+=dp[i][j]/n;
            }
        }
    }
    
    Matrix<double> B(sz,sz+1);
    {
        B[0][0]=1;
        B[0][sz]=0;
        for(int i=1;i<sz;i++){
            B[i][i]=-1;
            for(int j=0;j<sz;j++){
                // i to j
                int d=i-j;
                if(0<=d && d<sz){
                    B[i][j]+=dp[k][d];
                    B[i][sz]-=dp[k][d];
                }
                d=i+j;
                if(j!=0 && 0<=d && d<sz){
                    B[i][j]+=dp[k][d];
                    B[i][sz]-=dp[k][d];
                }
            }
        }
    }
    B=gaussianElimination(B);
    Vector<double> vec(sz+1);
    for(int i=sz-1;i>=0;i--){
        for(int j=0;j<i;j++){
            B[j][sz]-=B[j][i]/B[i][i]*B[i][sz];
            B[j][i]=0;
        }
    }
    for(int i=0;i<sz;i++){
        vec[sz-i-1]=B[i][sz]/B[i][i];
    }
    vec[sz]=1;

    Matrix<double> A(sz+1,sz+1);
    for(int i=1;i<sz;i++){
        A[i][i-1]=1;
    }
    A[sz][sz]=1;
    for(int i=0;i<sz;i++){
        int d=i+1;
        if(d<sz) A[0][i]=dp[k][d];
        if(d<sz) A[0][sz]+=dp[k][d];
    }
    A=pow(A,abs(s));
    double res=0;
    for(int i=0;i<=sz;i++){
        res+=A[sz-1][i]*vec[i];
    }
    cout<<fixed<<setprecision(10);
    cout<<res<<endl;

    return 0;
}


Compiler output:
a.cc:3:17: error: 'ProconLib' is not a namespace-name
    3 | using namespace ProconLib;
      |                 ^~~~~~~~~
a.cc:145:25: error: redefinition of default argument for 'class Detail'
  145 |     template<typename T,typename Detail=Detail<T>>
      |                         ^~~~~~~~
a.cc:61:25: note: original definition appeared here
   61 |     template<typename T,typename Detail=Detail<T>>
      |                         ^~~~~~~~
a.cc: In function 'int main()':
a.cc:244:5: error: 'Matrix' was not declared in this scope; did you mean 'ProconLib::Matrix'?
  244 |     Matrix<double> B(sz,sz+1);
      |     ^~~~~~
      |     ProconLib::Matrix
a.cc:33:11: note: 'ProconLib::Matrix' declared here
   33 |     class Matrix{
      |           ^~~~~~
a.cc:244:12: error: expected primary-expression before 'double'
  244 |     Matrix<double> B(sz,sz+1);
      |            ^~~~~~
a.cc:246:9: error: 'B' was not declared in this scope
  246 |         B[0][0]=1;
      |         ^
a.cc:265:5: error: 'B' was not declared in this scope
  265 |     B=gaussianElimination(B);
      |     ^
a.cc:265:7: error: 'gaussianElimination' was not declared in this scope; did you mean 'ProconLib::gaussianElimination'?
  265 |     B=gaussianElimination(B);
      |       ^~~~~~~~~~~~~~~~~~~
      |       ProconLib::gaussianElimination
a.cc:146:15: note: 'ProconLib::gaussianElimination' declared here
  146 |     Matrix<T> gaussianElimination(Matrix<T> mat){
      |               ^~~~~~~~~~~~~~~~~~~
a.cc:266:5: error: 'Vector' was not declared in this scope; did you mean 'ProconLib::Vector'?
  266 |     Vector<double> vec(sz+1);
      |     ^~~~~~
      |     ProconLib::Vector
a.cc:16:11: note: 'ProconLib::Vector' declared here
   16 |     class Vector{
      |           ^~~~~~
a.cc:266:12: error: expected primary-expression before 'double'
  266 |     Vector<double> vec(sz+1);
      |            ^~~~~~
a.cc:274:9: error: 'vec' was not declared in this scope
  274 |         vec[sz-i-1]=B[i][sz]/B[i][i];
      |         ^~~
a.cc:276:5: error: 'vec' was not declared in this scope
  276 |     vec[sz]=1;
      |     ^~~
a.cc:278:12: error: expected primary-expression before 'double'
  278 |     Matrix<double> A(sz+1,sz+1);
      |            ^~~~~~
a.cc:280:9: error: 'A' was not declared in this scope
  280 |         A[i][i-1]=1;
      |         ^
a.cc:282:5: error: 'A' was not declared in this scope
  282 |     A[sz][sz]=1;
      |     ^