Language: C++

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define ft first
#define sd second
#define fr(i,n) for(int i=0;i<(n);++i)
#define Fr(i,n) for(int i=1;i<=(n);++i)
#define ifr(i,n) for(int i=(n)-1;i>=0;--i)
#define iFr(i,n) for(int i=(n);i>0;--i)

template<class T,class U>
T pw(T n,U m){
    if(m==0) return T(1);
    else if(m%2==0) return pw(n*n,m/2);
    else return n*pw(n,m-1);
}

ll MOD;
struct modint{
    ll a;
    modint(ll a_=0){
        a=((a_%MOD)+MOD)%MOD;
    }
    modint inv()const{
        ll n=1,m=MOD-2,A=a;
        while(m){
            if(m&1)(n*=A)%=MOD;
            (A*=A)%=MOD;
            m>>=1;
        }
        modint y(n);
        return y;
    }
    bool operator==(const modint& x){
        return a==x.a;
    }
    bool operator!=(const modint& x){
        return a!=x.a;
    }
    modint& operator=(const modint& x){
        a=x.a;
        return *this;
    }
    modint operator+(const modint& x){
        modint y(a+x.a);
        return y;
    }
    modint operator-(const modint& x){
        modint y(a-x.a);
        return y;
    }
    modint operator*(const modint& x){
        modint y(a*x.a);
        return y;
    }
    modint operator/(const modint& x){
        return *this * x.inv();
    }
    modint& operator+=(const modint& x){
        *this=*this+x;
        return *this;
    }
    modint& operator-=(const modint& x){
        *this=*this-x;
        return *this;
    }
    modint& operator*=(const modint& x){
        *this=*this*x;
        return *this;
    }
    modint& operator/=(const modint& x){
        *this=*this/x;
        return *this;
    }
};
istream& operator>>(istream &in,modint& x){
    ll a_;
    in>>a_;
    modint y(a_);
    x=y;
    return in;
}
ostream& operator<<(ostream &out,const modint& x){
    out<<x.a;
    return out;
}
modint pwr(ll a,ll b){
    modint n(1),A(a);
    while(b){
        if(b&1) n*=A;
        A*=A;
        b>>=1;
    }
    return n;
}

template<class T> struct mat{
    int n,m;
    vector<vector<T>> v;
    const ll MOD=1e9+7;
    mat(){}
    mat(int n_,int m_,T e){
        n=n_;m=m_;
        v.clear();
        v.resize(n);
        fr(i,n){
            v[i].resize(m);
            v[i][i]=e;
        }
    }
    mat(vector<vector<T>>& v_){
        n=v_.size();
        m=v_[0].size();
        swap(v,v_);
    }
    bool operator==(const mat& x){
        if(n!=x.n||m!=x.m) return 0;
        fr(i,n) fr(j,m) if(v[i][j]!=x.v[i][j]) return 0;
        return 1;
    }
    mat& operator=(const mat& x){
        n=x.n;
        m=x.m;
        v=vector<vector<T>>(n);
        copy(x.v.begin(),x.v.end(),v.begin());
        return *this;
    }
    mat operator+(const mat& x){
        if(n!=x.n||m!=x.m) return mat();
        vector<vector<T>> y(n,vector<T>(m));
        fr(i,n) fr(j,m) y[i][j]=v[i][j]+x.v[i][j];
        return mat(y);
    }
    mat operator-(const mat& x){
        if(n!=x.n||m!=x.m) return mat();
        vector<vector<T>> y(n,vector<T>(m));
        fr(i,n) fr(j,m) y[i][j]=v[i][j]-x.v[i][j];
        return mat(y);
    }
    mat operator*(const mat& x){
        if(m!=x.n) return mat();
        vector<vector<T>> y(n,vector<T>(x.m));
        fr(i,n) fr(k,m) fr(j,x.m) y[i][j]+=v[i][k]*x.v[k][j];
        return mat(y);
    }
    mat& operator+=(const mat& x){
        if(n!=x.n||m!=x.m) return mat();
        fr(i,n) fr(j,m) v[i][j]+=x.v[i][j];
        return *this;
    }
    mat& operator-=(const mat& x){
        if(n!=x.n||m!=x.m) return mat();
        fr(i,n) fr(j,m) v[i][j]-=x.v[i][j];
        return *this;
    }
    mat& operator*=(const mat& x){
        *this=*this*x;
        return *this;
    }
    T det(){
        if(n!=m) return T();
        mat y=*this;
        int c=1;
        for(int i=0;i<n;++i){
            for(int j=i;j<n;++j){
                if(v[j][i]!=T()){
                    if(i!=j){
                        swap(v[i],v[j]);
                        c*=-1;
                    }
                    break;
                }
                if(j==n-1) return T();
            }
            for(int j=i+1;j<n;++j){
                T ka=v[j][i]/v[i][i];
                for(int k=i;k<n;k++){
                    v[j][k]-=v[i][k]*ka;
                }
            }
        }
        T d=T(1);
        for(int i=0;i<n;i++) d*=v[i][i];
        d*=c;
        *this=y;
        return d;
    }
};

int main(){
    cin.tie(nullptr);
    ios::sync_with_stdio(false);
    ll l,a,b;
    cin>>l>>a>>b>>MOD;
    modint ans;
    int i=0,t;
    vector<ll> ind(40);
    ll c=1;
    fr(k,40){
        if(k<=18){
            if(a>c){
                c*=10;
                continue;
            }
        }
        ll L=0,R=1e18,M;
        while(R-L>1){
            M=(L+R)/2;
            long double ka=log10((long double)a/M+b)+log10((long double)M);
            if(k<=ka) R=M;
            else L=M;
        }
        ind[k]=min(l,R);
        c*=10;
    }
    modint tmp(1);
    ifr(k,40){
        if(k>0&&ind[k]==ind[k-1]) continue;
        if(ind[k]==0) break;
        int lg=ind[k];
        if(k>0) lg-=ind[k-1];
        cout<<k<<" "<<lg<<endl;
        modint cn=pwr(10,k),D;
        if(k) D=modint(a)+modint(b)*modint(ind[k-1]);
        else D=modint(a);
        vector<vector<modint>> vi(3,vector<modint>(3)),vi2(2,vector<modint>(2));
        vi[0][0]=vi2[0][0]=cn;
        vi[0][1]=vi[1][1]=vi[1][2]=vi[2][2]=vi2[0][1]=vi2[1][1]=1;
        mat<modint> mt(vi),mt2(vi2),an=pw(mt,lg),an2=pw(mt2,lg);
        ans+=tmp*(D*mt2.v[0][1]+mt.v[0][1]*b);
        tmp*=pw(pwr(10,lg),k);
    }
    cout<<ans<<endl;
}

Compiler output:
a.cc: In instantiation of 'T pw(T, U) [with T = mat<modint>; U = int]':
a.cc:229:42:   required from here
  229 |         mat<modint> mt(vi),mt2(vi2),an=pw(mt,lg),an2=pw(mt2,lg);
      |                                        ~~^~~~~~~
a.cc:13:21: error: no matching function for call to 'mat<modint>::mat(int)'
   13 |     if(m==0) return T(1);
      |                     ^~~~
a.cc:111:5: note: candidate: 'mat<T>::mat(std::vector<std::vector<_Tp> >&) [with T = modint]'
  111 |     mat(vector<vector<T>>& v_){
      |     ^~~
a.cc:111:28: note:   no known conversion for argument 1 from 'int' to 'std::vector<std::vector<modint> >&'
  111 |     mat(vector<vector<T>>& v_){
      |         ~~~~~~~~~~~~~~~~~~~^~
a.cc:102:5: note: candidate: 'mat<T>::mat(int, int, T) [with T = modint]'
  102 |     mat(int n_,int m_,T e){
      |     ^~~
a.cc:102:5: note:   candidate expects 3 arguments, 1 provided
a.cc:101:5: note: candidate: 'mat<T>::mat() [with T = modint]'
  101 |     mat(){}
      |     ^~~
a.cc:101:5: note:   candidate expects 0 arguments, 1 provided
a.cc:97:26: note: candidate: 'mat<modint>::mat(const mat<modint>&)'
   97 | template<class T> struct mat{
      |                          ^~~
a.cc:97:26: note:   no known conversion for argument 1 from 'int' to 'const mat<modint>&'