Problem:
 Score : 600 points    Problem Statement  There is an arithmetic progression with L terms: s_0, s_1, s_2, ... , s_{L-1} .  The initial term is A , and the common difference is B . That is, s_i = A + B \times i holds.  Consider the integer obtained by concatenating the terms written in base ten without leading zeros. For example, the sequence 3, 7, 11, 15, 19 would be concatenated into 37111519 . What is the remainder when that integer is divided by M ?      Constraints   All values in input are integers.  1 \leq L, A, B < 10^{18}  2 \leq M \leq 10^9  All terms in the arithmetic progression are less than 10^{18} .         Input Input is given from Standard Input in the following format:  L  A  B  M       Output  Print the remainder when the integer obtained by concatenating the terms is divided by M .        Sample Input 1 5 3 4 10007      Sample Output 1 5563  Our arithmetic progression is 3, 7, 11, 15, 19 , so the answer is 37111519 mod 10007 , that is, 5563 .       Sample Input 2 4 8 1 1000000      Sample Output 2 891011       Sample Input 3 107 10000000000007 1000000000000007 998244353      Sample Output 3 39122908

Solution (Python):
def examA():
    import numpy as np
    def powDetMod(n,arr,k,MOD=-1):
        repeat = 60
        local_doubling = [[] for _ in range(repeat+1)]
        local_doubling[0] = arr
        for i in range(repeat):
            local_doubling[i+1] = np.dot(local_doubling[i],local_doubling[i])
            if MOD>0:
                for h in range(n):
                    for w in range(n):
                        local_doubling[i+1][h][w] %= MOD
        now = [[1 if i==j else 0 for i in range(n)]for j in range(n)]
        for i in range(repeat+1):
            if k&(1<<i)>0:
                now = np.dot(now,local_doubling[i])
                if MOD > 0:
                    for h in range(n):
                        for w in range(n):
                            now[h][w] %= MOD
        return now
    L, A, B, M = LI()
    C = [0]*18
    cnt = L
    for i in range(18):
        if (pow(10,(i+1))-1-A)<0:
            continue
        C[i] = 1 + (pow(10,(i+1))-1-A)//B
        if (pow(10,i)-1-A)>=0:
            C[i] -= 1 + (pow(10,i)-1-A)//B
        if C[i]>=cnt:
            C[i] = cnt
            break
        cnt -= C[i]
    
    X =  np.array([0,A,B])
    for i in range(18):
        now = [[pow(10,i+1,M),0,0],[1,1,0],[0,1,1]]
        now = np.array(now)
        arr = powDetMod(3,now,C[i],M)
        X = np.dot(X,arr)%M

    ans = X[0]%M
    print(ans)
    return

def examB():
    N = I()
    R = [inf,-inf]
    L = [inf,-inf]
    U = [inf,-inf]
    D = [inf,-inf]
    X = [inf,-inf]
    Y = [inf,-inf]
    for _ in range(N):
        x,y,s = LSI()
        x,y = int(x),int(y)
        if s=="R":
            if x<R[0]: R[0] = x
            if x>R[1]: R[1] = x
            if y<Y[0]: Y[0] = y
            if y>Y[1]: Y[1] = y
        elif s=="L":
            if x<L[0]: L[0] = x
            if x>L[1]: L[1] = x
            if y<Y[0]: Y[0] = y
            if y>Y[1]: Y[1] = y
        elif s=="U":
            if y<U[0]: U[0] = y
            if y>U[1]: U[1] = y
            if x<X[0]: X[0] = x
            if x>X[1]: X[1] = x
        elif s=="D":
            if y<D[0]: D[0] = y
            if y>D[1]: D[1] = y
            if x<X[0]: X[0] = x
            if x>X[1]: X[1] = x
    def calc_area(t):
        x = max(R[1]+t,L[1]-t,X[1])-min(R[0]+t,L[0]-t,X[0])
        y = max(U[1]+t,D[1]-t,Y[1])-min(U[0]+t,D[0]-t,Y[0])
        if x<=0 or y<=0:
            return 0
        return x*y
    ans = inf
    for i in [0,1]:
        for j in [0,1]:
            if R[i] <= X[j]:
                T = (X[j] - R[i])
                ans = min(ans, calc_area(T))
            if L[i] >= X[j]:
                T = (L[i] - X[j])
                ans = min(ans, calc_area(T))
            if R[i] <= L[j]:
                T = (L[j] - R[i]) / 2
                ans = min(ans, calc_area(T))

            if U[i] <= Y[j]:
                T = (Y[j] - U[i])
                ans = min(ans, calc_area(T))
            if D[i] >= Y[j]:
                T = (D[i] - Y[j])
                ans = min(ans, calc_area(T))
            if U[i] <= D[j]:
                T = (D[j] - U[i]) / 2
                ans = min(ans, calc_area(T))

    print(ans)
    return

from decimal import Decimal as dec
import sys,bisect,itertools,heapq,math,random
from copy import deepcopy
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
def I(): return int(input())
def LI(): return list(map(int,sys.stdin.readline().split()))
def DI(): return dec(input())
def LDI(): return list(map(dec,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet,_ep
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
_ep = 10**(-12)
alphabet = [chr(ord('a') + i) for i in range(26)]

sys.setrecursionlimit(10**7)

if __name__ == '__main__':
    examA()

"""
142
12 9 1445 0 1
asd dfg hj o o
aidn
"""

Solution (Python):
n, a, b, MOD = map(int, input().split())


def mult(A, B, N):
    C = [[0] * N for _ in range(N)]
    for i in range(N):
        for j in range(N):
            for k in range(N):
                C[i][j] += A[i][k] * B[k][j]
            C[i][j] %= MOD
    return C


def mat_power(A, N, n):
    if n == 0:
        res = [[0] * N for _ in range(N)]
        for i in range(N):
            res[i][i] = 1
        return res
    X = mat_power(A, N, n // 2)
    X = mult(X, X, N)
    return mult(A, X, N) if n & 1 else X


li = [0]
t = 1
while li[-1] < n:
    p = pow(10, t)
    l, r = li[-1] - 1, n
    while r - l > 1:
        m = (l + r) // 2
        if a + b * m < p:
            l = m
        else:
            r = m
    li.append(r)
    t += 1
ans = 0
for i in range(len(li) - 1):
    p = pow(10, i + 1)
    t = li[i + 1] - li[i]
    am = [[p, 1, 0], [0, 1, 1], [0, 0, 1]]
    mat = mat_power(am, 3, t)
    ans *= pow(10, t * (i + 1), MOD)
    ans += mat[0][1] * ((a + b * li[i])) + mat[0][2] * b
    ans %= MOD
print(ans)

Solution (Python):
import numpy as np
import sys

def II(): return int(sys.stdin.readline())
def MI(): return map(int, sys.stdin.readline().split())
def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(rows_number): return [LI() for _ in range(rows_number)]
def SI(): return sys.stdin.readline()[:-1]


def main():
    l,a,b,md=MI()

    
    k=0
    i=0
    ans=[0,a%md,b%md]
    ans=np.array(ans,dtype="i8")
    
    while i<l:
        k+=1
        j=(10**k-a+b-1)//b
        if j<=0:continue
        if j>l:j=l
        
        e=j-i
        
        
        
        
        
        bb=np.zeros((3,3),dtype="i8")
        bb[0,0]=pow(10,k,md)
        bb[1,0]=bb[1,1]=bb[2,1]=bb[2,2]=1
        while e:
            if e&1:ans=np.dot(ans,bb)%md
            bb=np.dot(bb,bb)%md
            e>>=1
        i=j

    print(ans[0])

main()

Solution (Python):

import sys
from collections import deque
import numpy as np
sys.setrecursionlimit(10**9)
INF=10**18
MOD=10**9+7
input=lambda: sys.stdin.readline().rstrip()
YesNo=lambda b: bool([print('Yes')] if b else print('No'))
YESNO=lambda b: bool([print('YES')] if b else print('NO'))
int1=lambda x:int(x)-1

def main():
    L,A,B,M=map(int,input().split())
    S=lambda i:A+B*i
    sec=[]
    
    def nibutan(ok,ng,k):
        while abs(ok-ng) > 1:
            mid = (ok + ng) // 2
            if solve(mid,k):
                ok = mid
            else:
                ng = mid
        return ok
    
    def solve(mid,k):
        return len(str(S(mid)))<k
    
    left=L-1
    for k in range(19,0,-1):
        if len(str(S(left)))==k:
            right=nibutan(-1,left,k)
            sec.append((k,left,right))
            left=right
    
    X=np.array([0,A%M,1],dtype=np.int64)
    for k,left,right in sec[::-1]:
        n=left-right
        doubling=[np.array([[pow(10,k,M),0,0],[1,1,0],[0,B%M,1]],dtype=np.int64)]
        for i in range((10**18).bit_length()):
            doubling.append(np.dot(doubling[i],doubling[i])%M)
        n=left-right
        for i in range(len(doubling)):
            if n&1:
                X=np.dot(X,doubling[i])%M
            n>>=1
            if n==0:
                break
    print(X[0])

if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np
l, a, b, m = map(int, input().split())





s_max = a + b * (l - 1)


ten_pow = [1]
while ten_pow[-1] <= s_max:
    ten_pow.append(ten_pow[-1] * 10)


D_MAX = len(ten_pow)
d = [0] * D_MAX
for i in range(D_MAX):
    if ten_pow[i] >= a:
        d[i] = max((ten_pow[i] - 1 - a) // b, 0)
d[-1] = l - 1


cd = [d[i + 1] - d[i] for i in range(len(d) - 1)]




y = np.array([a, a + b, 1]).T % m
for i in range(len(cd)):
    x = np.array([[ten_pow[i + 1] % m, 1, 0], [0, 1, b], [0, 0, 1]])
    x %= m
    
    
    k = cd[i]
    while k > 0:
        if k & 1:
            y = np.dot(x, y) % m
            
        x = np.dot(x, x) % m
        k = k >> 1

ans = y[0] % m
print(ans)

Solution (Python):
L,A,B,M=map(int,input().split())
MOD=M



def imat(n):
  ret=[[1 if i==j else 0 for j in range(n)] for i in range(n)]
  return ret


def prod_mat(amat,bmat):
  res_mat = [[sum([amat[i][j]*bmat[j][k] for j in range(len(bmat))]) for k in range(len(bmat[0]))] for i in range(len(amat))]
  return res_mat  

def powmod_mat(amat,p):
  if p==0:
    return imat(len(amat))
  else:
    pow2=powmod_mat(amat,p//2)
    if p%2==0:
      res_mat=prod_mat(pow2,pow2)
    else:
      res_mat=prod_mat(amat,prod_mat(pow2,pow2))

    for i in range(len(amat)):
      for j in range(len(amat)):
        res_mat[i][j]%=MOD
    return res_mat

slist=[]
for d in range(1,19):
  s=-(-(10**d-A)//B)
  if s>=0:
    slist.append(min(s,L))
  else:
    slist.append(0)


dlist=[slist[0]]
for d in range(1,18):
  dlist.append(slist[d]-slist[d-1])


mat=imat(3)
for d in range(18):
  e=dlist[d]
  pmat=powmod_mat([[10**(d+1),1,0],[0,1,B],[0,0,1]],e)
  mat=prod_mat(pmat,mat)
  for i in range(3):
    for j in range(3):
      mat[i][j]%=MOD


vec=[[0],[A],[1]]
vec=prod_mat(mat,vec)
print(vec[0][0]%MOD)

Solution (Python):
import sys
import time

readline = sys.stdin.readline
MOD = 10 ** 9 + 7
INF = float('INF')
sys.setrecursionlimit(10 ** 5)


def main():
    import numpy as np
    l, a, b, m = list(map(int, readline().split()))
    d_min = len(str(a))
    d_max = len(str(a + (l - 1) * b))
    count_d = [-1] * 19

    prev = -1
    for i in range(d_min, d_max + 1):
        power = 10 ** i
        ok = 0
        ng = l + 1
        while abs(ng - ok) > 1:
            mid = (ok + ng) // 2
            x = a + mid * b
            if x < power:
                ok = mid
            else:
                ng = mid
        ok = min(ok, l - 1)
        count_d[i] = ok - prev
        prev = ok

    res = np.array([0, a, 1], dtype="int64")
    res = np.mod(res, m)

    for i in range(19):
        count = count_d[i]
        if count == -1:
            continue
        mat = np.array([[10 ** i, 0, 0],
                        [1, 1, 0],
                        [0, b, 1]], dtype="int64")
        mat = np.mod(mat, m)
        cur = np.eye(3, dtype="int64")
        count_bit = format(count, "b")
        for bit in count_bit[::-1]:
            if bit == "1":
                cur = np.dot(cur, mat)
                cur = np.mod(cur, m)
            mat = np.dot(mat, mat)
            mat = np.mod(mat, m)
        res = np.dot(res, cur)
        res = np.mod(res, m)

    print(res[0])


if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np

L, A, B, M = map(int, input().split())

def dot(A, B, mod):
    r = len(A)
    c = len(B[0])
    m = len(A[0])
    res = [[0] * c for _ in range(r)]
    
    for i in range(r):
        for j in range(c):
            for k in range(m):
                res[i][j] += A[i][k] * B[k][j] % mod
                res[i][j] %= mod
    
    return res

def rec_pow(k, n, mod):
    if n == 1:
        return k
    
    k2 = dot(k, k, mod)
    if n % 2 ==0:
        return rec_pow(k2, n//2, mod)
    else:
        return dot(rec_pow(k2, n//2, mod), k, mod)

s = [[0, A%M, 1]]

for d in range(1, 19):
    n1 = min(max(0, (10**d-1-A)//B+1), L)
    n0 = min(max(0, (10**(d-1)-1-A)//B+1), L)
    n = n1 - n0
    if n == 0:
        continue
        
    k = [[pow(10, d, M), 0, 0], [1, 1, 0], [0, B%M, 1]]
    
    s = dot(s, rec_pow(k, n, M), M)
    
    if n1 >= L:
        break
        
print(s[0][0])

Solution (Python):
L,A,B,M=map(int,input().split())

def partsum(a,b,l,n):
    doublingconst=[0 for i in range(0,60)]
    doublingconst[0]=1
    doublingline=[0 for i in range(0,60)]
    doublingline[0]=0
    for i in range(1,60):
        doublingconst[i]=((pow(10,l*2**(i-1),M)+1)*doublingconst[i-1]%M)
        doublingline[i]=((pow(10,l*2**(i-1),M)+1)*doublingline[i-1]+pow(2,i-1,M)*pow(10,l*2**(i-1),M)*doublingconst[i-1])%M

    ansconst=0
    chousei=0
    for i in range(0,60):
        if n>>i &1==1:
            ansconst+=doublingconst[i]*pow(10,chousei*l,M)
            ansconst%=M
            chousei+=2**i

    ansline=0
    chousei=0
    for i in range(0,60):
        if n>>i &1==1:
            ansline+=(doublingline[i]+chousei*doublingconst[i])*pow(10,chousei*l,M)
            ansline%=M
            chousei+=2**i

    return ansline*a+ansconst*b

start=len(str(A))
end=len(str(B*L+A-B))
part=[]
for i in range(start-1,end):
    l=1+(10**i-A-1)//B
    r=(10**(i+1)-A-1)//B
    l=max(0,l)
    r=min(L-1,r)
    part.append([i+1,l,r])

part.sort(reverse=True)
ans=0
chousei=0
for i in range(0,len(part)):
    length,l,r=part[i]
    const=B*r+A
    line=B
    ans+=(partsum(-line,const,length,r-l+1)*pow(10,chousei,M))%M
    ans%=M
    chousei+=length*(r-l+1)
print(ans)

Solution (Python):
l,a,b,m = map(int,input().split())

ten = 1
def f(x):

    if x == 0:
        return 0
    if x % 2 == 1:
        t = x-1
        x = f(t)
        return x*ten+1
    else:
        t = x//2
        x = f(t)
        return x*pow(ten,t,m) + x

def g(x):

    if x == 0:
        return 0
    if x % 2 == 1:
        t = x-1
        x = g(t)
        return x*ten + b*(t)
    else:
        t = x//2
        x = g(t)
        return x*pow(ten,t,m) + x + b*(t)*f(t)
last = a+b*(l-1)

ans = 0
for i in range(1,20):
    le = ten
    r = ten*10-1
    
    ten *= 10
    if last < le or a > r:
        continue
    if le < a:
        na = a
    else:
        na = (le-a+b-1)//b*b+a
    if last <= r:
        nl = last
    else:
        nl = (r-a)//b*b+a
    n = (nl-na)//b+1
    
    ans *= pow(ten,n,m)
    ans += ((f(n)%m)*na)%m
    ans += g(n)%m
    ans %= m

print(ans)

Solution (Python):

"""
Writer: SPD_9X2
https://atcoder.jp/contests/abc129/tasks/abc129_f
=====解説をみた=====

なるほど行列式ね…
繰り返し同じ写像を繰り返す場合、行列式による繰り返し2乗法
を考えるべきなのかもしれんね…

"""


def matrix_mul(A,B):

    ans = [ [0] * len(B[0]) for i in range(len(A)) ]

    for ai in range(len(A)):

        for bj in range(len(B[0])):

            now = 0
            for same in range(len(A[0])):
                now += A[ai][same] * B[same][bj]
            ans[ai][bj] = now % mod

    return ans


def matrix_pow(A,x):

    B = [[A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
    ans = [[0] * len(A[0]) for i in range(len(A))]
    for i in range(len(A)):
        ans[i][i] = 1

    while x > 0:
        if x % 2 == 1:
            ans = matrix_mul(ans,B)
        B = matrix_mul(B,B)

        x//=2
    return ans

import sys
L,A,C,mod = map(int,input().split())

now = [[0,A,1]]
LA = len(str(A))
LB = len(str(A+C*(L-1)))
B = A+C*(L-1)

if LA==LB:
    np = [[pow(10,LA,mod),0,0],[1,1,0],[0,C,1]]
    k = (B-A)//C+1
    now = matrix_mul(now , matrix_pow(np,k))
    print (now[0][0] % mod)
    sys.exit()

for i in range(LA,LB+1):

    np = [[pow(10,i,mod),0,0],[1,1,0],[0,C,1]]
    

    if i == LA:
        k = (10**i-1-A)//C+1
        now = matrix_mul(now , matrix_pow(np,k))
    elif i == LB:
        k = (B-A)//C-(10**(i-1)-1-A)//C
        now = matrix_mul(now , matrix_pow(np,k))
    else:
        k = (10**i-1-A)//C - (10**(i-1)-1-A)//C
        now = matrix_mul(now , matrix_pow(np,k))

    
print (now[0][0] % mod)

Solution (Python):
import sys
import itertools
import queue
input = sys.stdin.readline
sys.setrecursionlimit(100000)


def read_values(): return map(int, input().split())
def read_index(): return map(lambda x: int(x) - 1, input().split())
def read_list(): return list(read_values())
def read_lists(N): return [read_list() for n in range(N)]


class V:
    def __init__(self, f, v=None):
        self.f = f
        self.v = v
 
    def __str__(self):
        return str(self.v)
 
    def ud(self, n):
        if self.v is None:
            self.v = n
        else:
            self.v = self.f(self.v, n) 


def mat_mul(X, Y, mod):
    m = len(X)
    return [[sum(X[i][k] * Y[k][j] % mod for k in range(m)) % mod for j in range(m)] for i in range(m)]


def mat_pow(M, k, mod):
    m = len(M)
    res = [[1 if i == j else 0 for j in range(m)] for i in range(m)]
    while k > 0:
        if k & 1:
            res = mat_mul(M, res, mod)
        M = mat_mul(M, M, mod)
        k >>= 1

    return res


def main():
    L, A, B, M = read_values()
    res = 0
    D = [(0, 0)]
    for d in range(1, 20):
        l = D[-1][1]
        r = max(0, min((10 ** d - 1 - A) // B + 1, L))
        D.append((l, r))
    
    res = 0
    c = 0
    R = [0, A, 1]
    Mat = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    for d in range(1, 20):
        num = D[d][1] - D[d][0]

        MM = [[10 ** d, 0, 0], [1, 1, 0], [0, B, 1]]
        Mat = mat_mul(Mat, mat_pow(MM, num, M), M)

    X = sum(R[i] * Mat[i][0] % M for i in range(3)) % M

    print(X)


if __name__ == "__main__":
    main()

Solution (Python):
def main():
    import sys
    input = sys.stdin.readline

    N, a0, d, mod = map(int, input().split())

    def matmul(A, B):
        C = [[0] * len(B[0]) for _ in range(len(A))]
        for i in range(len(A)):
            for j in range(len(B[0])):
                for k in range(len(B)):
                    C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod
        return C

    def matpow(A, p):
        n = len(A)
        B = [[0] * n for _ in range(n)]
        for i in range(n):
            B[i][i] = 1
        while p > 0:
            if p & 1:
                B = matmul(B, A)
            A = matmul(A, A)
            p >>= 1
        return B

    ans = [[a0%mod], [a0%mod], [d%mod]]
    ok_prev = 0
    for i in range(len(str(a0)), 20):
        ok = N
        ng = -1
        mid = (ok+ng)//2
        while ok - ng > 1:
            if len(str(a0 + d * mid)) >= i+1:
                ok = mid
            else:
                ng = mid
            mid = (ok+ng)//2

        mat = [[(10**i)%mod, 1, 1],
               [0, 1, 1],
               [0, 0, 1]]
        ans = matmul(matpow(mat, ok - ok_prev-1), ans)
        ok_prev = ok-1
        if ok == N:
            break
    print(ans[0][0])


if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np

L, A, B, M = (int(i) for i in input().split())
keta = [0] * 18
for i in range(18):
    keta[i] = min((10**i - A - 1)//B, L-1)
keta.reverse()
for i in range(18):
    if keta[i] < 0:
        keta[i] = -1
        break
keta.reverse()

def power(g, N, M):
    if N == 0:
        g = np.array([1, 0, 0, 0, 1, 0, 0, 0, 1]).reshape(3, 3)
    elif N == 1:
        pass
    elif N%2 == 0:
        g = np.linalg.matrix_power(power(g, N//2, M), 2)
    elif N%2 == 1:
        g = np.dot(g, power(g, N-1, M))
    return g%M

ans = 0
for i in range(18):
    N = i + 1
    mi = keta[i] + 1
    if mi < 0:
        continue
    if i != 17:
        ma = keta[i+1]
    else:
        ma = L-1
    if mi >= L:
        break
    t = np.array([ans, (A + B*mi)%M, 1])
    g = np.array([10**N%M, 0, 0, 1, 1, 0, 0, B%M, 1]).reshape(3, 3)
    g = power(g, ma-mi+1, M)
    g = np.dot(t, g)
    ans = g[0]%M
print(ans)

Solution (Python):
from numpy.linalg import matrix_power
from numpy import dot, array

class mint:
    def __init__(self, n):
        self.n = n % mod

    def __add__(self, other):
        return mint(self.n + other.n)

    def __mul__(self, other):
        return mint(self.n * other.n)

    def __repr__(self):
        return str(self.n)


L, A, B, mod = map(int, input().split())

ans = array([mint(0), mint(A), mint(B)])

size = 0
for d in range(1, 19):
    x = pow(10, d, mod)
    T = array([
        [mint(x), mint(0), mint(0)],
        [mint(1), mint(1), mint(0)],
        [mint(0), mint(1), mint(1)],
    ])

    cnt = min(L - size, max(0, 10 ** d - 1 - A) // B - size + 1)
    if cnt > 0:
        ans = dot(ans, matrix_power(T, cnt))

    size += cnt
    if size >= L:
        break

print(ans[0])

Solution (Python):
from numpy.linalg import matrix_power
from numpy import dot

class mint:
    def __init__(self, n):
        self.n = n % mod

    def __add__(self, other):
        return mint(self.n + other.n)

    def __mul__(self, other):
        return mint(self.n * other.n)

    def __repr__(self):
        return str(self.n)


L, A, B, mod = map(int, input().split())

ans = [mint(0), mint(A), mint(B)]

size = 0
for d in range(1, 19):
    x = pow(10, d, mod)
    T = [
        [mint(x), mint(0), mint(0)],
        [mint(1), mint(1), mint(0)],
        [mint(0), mint(1), mint(1)],
    ]

    cnt = min(L - size, max(0, 10 ** d - 1 - A) // B - size + 1)
    if cnt > 0:
        ans = dot(ans, matrix_power(T, cnt))

    size += cnt
    if size >= L:
        break

print(ans[0])

Solution (Python):
l, a, b, mod = map(int, input().split())

def mul(m1, m2):
    m = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
    for i in range(3):
        for j in range(3):
            m[i][j] = sum(m1[i][k] * m2[k][j] for k in range(3)) % mod
    return m






mat = [[0, 0, 0],
       [a, 0, 0],
       [1, 0, 0]]
d = 1
e = a + b * (l - 1)
while e >= d:
    d *= 10
    if a >= d:
        continue
    t = (min(d - 1, e) - a) // b + 1
    a += b * t
    bit = [[d, 1, 0],
           [0, 1, b],
           [0, 0, 1]]
    while t > 0:
        if t & 1:
            mat = mul(bit, mat)
        bit = mul(bit, bit)
        t //= 2


print(mat[0][0])

Solution (Python):
import sys
L, A, B, M = map(int, input().split())
mod = M
sys.setrecursionlimit(pow(10, 8))

def power(x, y):
    if   y == 0     : return 1
    elif y == 1     : return x % mod
    elif y % 2 == 0 : return power(x, y//2)**2 % mod
    else            : return power(x, (y-1)//2)**2 * x % mod
def mul(a, b):
    return ((a % mod) * (b % mod)) % mod

def G(n, i):
    
    
    if n == 1:
        return 1
    tn = n//2
    return (G(tn, i) * (1 + power(10, tn*i)) %mod + (n%2)*power(10, (n-1)*i)) % mod

def G2(n, i):
    
    
    if n ==0 :
        return 0
    if n == 1:
        return power(10, i) % mod
    tn = n//2
    X = G2(tn, i)
    I = power(10, tn*i)
    Y = (X*(I+1)%mod + mul(mul(tn, I), (G(tn+1, i)-1))%mod) + mul(n%2, mul(n, power(10, n*i)))
    return Y%mod


def f(a, b, n, i):
    if n == 0:
        return 0
    aa = mul(a, G(n, i))
    bb = mul(b, G2(n-1, i))
    return (aa - bb) % mod

SL = A + (L-1) * B
ts = 0
ks = [0]
for i in range(len(str(SL))+1):
    
    ks.append(max(min((10**(i+1)-1),SL)-(A-B), 0) // B - ts)
    ts += ks[-1]
ti = 0
a = SL
R = 0




for i in range(len(str(SL))+1,0,-1):
    n = ks[i]
    R = (R + f(a, B, n, i)*power(10, ti)) % mod
    ti += i*n
    a -= n*B
print(R)

Solution (Python):
import sys
L, A, B, M = map(int, input().split())
mod = M
sys.setrecursionlimit(pow(10, 8))

def power(x, y):
    if   y == 0     : return 1
    elif y == 1     : return x % mod
    elif y % 2 == 0 : return power(x, y//2)**2 % mod
    else            : return power(x, (y-1)//2)**2 * x % mod
def mul(a, b):
    return ((a % mod) * (b % mod)) % mod
def modinv(a):
    b, u, v = mod, 1, 0
    while b:
        t = a//b
        a, u = a-t*b, u-t*v
        a, b, u, v = b, a, v, u
    u %= mod
    return u
def div(a, b):
    return a * modinv(b) % mod

def G(n, i):
    if n == 1:
        return 1
    tn = n//2
    return (G(tn, i) * (1 + power(10, tn*i)) %mod + (n%2)*power(10, (n-1)*i)) % mod

def G2(n, i):
    if n ==0 :
        return 0
    if n == 1:
        return power(10, i) % mod
    tn = n//2
    X = G2(tn, i)
    I = power(10, tn*i)
    Y = (X*(I+1)%mod + mul(mul(tn, I), (G(tn+1, i)-1))%mod) %mod
    Y = Y + mul(n%2, mul(n, power(10, n*i)))
    return Y%mod


def f(a, b, n, i):
    if n == 0:
        return 0
    r = power(10, i)
    P = power(r,n)
    
    aa = mul(a, G(n, i))
    bb2 = mul(b, div(r-n*P + (n-1)*P*r, mul(1-r, 1-r)))
    bb = mul(b, G2(n-1, i))
    
    
    return (aa - bb) % mod

sl = A + (L-1) * B
t = 10
k = 0
ts = 0
ks = []
mk = 18
for i in range(mk):
    k = max(min((t-1),sl)-(A-B), 0) // B - ts
    
    
    ts += k
    ks.append(k)
    t *= 10
ti = 0
ti2 = 0
b = B
R = 0
for i in range(mk):
    i = mk - i
    n = ks[i-1]
    a = sl - ti2*b
    if n:
        R = (R + (f(a, b, n, i))*power(10, ti)) % mod



    ti += i*n
    ti2 += n
print(R)

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor
from operator import mul
from functools import reduce


sys.setrecursionlimit(2147483647)
INF = 10 ** 20
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 1000000007



def mat_mul(a, b) :
    I, J, K = len(a), len(b[0]), len(b)
    c = [[0] * J for _ in range(I)]
    for i in range(I) :
        for j in range(J) :
            for k in range(K) :
                c[i][j] += a[i][k] * b[k][j]
            c[i][j] %= m
    return c


def mat_pow(x, n):
    y = [[0] * len(x) for _ in range(len(x))]

    for i in range(len(x)):
        y[i][i] = 1

    while n > 0:
        if n & 1:
            y = mat_mul(x, y)
        x = mat_mul(x, x)
        n >>= 1

    return y


l, a, b, m = LI()
d0 = 0
ret = [[0], [a], [1]]
for i in range(1, 19):
    if 10 ** i - 1 - a < 0:
        continue
    d1 = min((10 ** i - 1 - a) // b + 1, l)
    mat = [[10 ** i, 1, 0], [0, 1, b], [0, 0, 1]]
    ret = mat_mul(mat_pow(mat, d1 - d0), ret)
    if d1 == l:
        break
    d0 = d1


print(ret[0][0])

Solution (Python):
from operator import mul

def getMatrixProduct(Ass, Bss, MOD):
    BssTr = [list(Bs) for Bs in zip(*Bss)]  
    ansss = [[sum(map(mul, As, Bs)) % MOD for Bs in BssTr] for As in Ass]
    return ansss

def getMatrixPower(Ass, n, MOD):
    sizeA = len(Ass)
    ansss = [[0]*(sizeA) for _ in range(sizeA)]
    for i in range(sizeA):
        ansss[i][i] = 1
    while n:
        if n & 1:
            ansss = getMatrixProduct(ansss, Ass, MOD)
        Ass = getMatrixProduct(Ass, Ass, MOD)
        n //= 2
    return ansss

L, A, B, M = map(int, input().split())

As = [[A%M], [A%M], [1]]
iPrev = 0
d = 1
while True:
    i = (10**d-1-A) // B
    if i > L-1:
        i = L-1
    num = i-iPrev
    if num > 0:
        Css = [[10**d,1,B%M], [0,1,B%M], [0,0,1]]
        PowCss = getMatrixPower(Css, num, M)
        As = getMatrixProduct(PowCss, As, M)
    if i == L-1:
        break
    iPrev = max(0, i)
    d += 1

print(As[0][0])

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

L,A,B,MOD = map(int,read().split())

def geom_seq_sum(r,n,MOD):
    """
    return \sum_{i=0}^{n-1} r^i mod MOD
    """
    if r == 1:
        return n%MOD
    r1 = r-1
    return ((pow(r,n,r1*MOD)-1)//r1)% MOD

def geom_seq_sum_1(r,n,MOD):
    """
    return \sum_{i=0}^{n-1} ir^i mod MOD
    """
    if r == 1:
        return n*(n-1)//2%MOD
    r1 = r-1
    x = pow(r,n,r1*r1*MOD)
    return (-(x-1)//r1 + (n-1)*x + 1)//r1 % MOD

S = 0
for d in range(1,20):
    
    lower = 10**(d-1)
    upper = 10**d-1
    
    nlower = (lower-A+B-1)//B
    nupper = (upper-A)//B
    if nlower < 0:
        nlower = 0
    if nupper >=L:
        nupper = L-1
    if nlower > nupper:
        continue
    items = nupper - nlower + 1
    r = (10**d) % MOD
    last = A+B*nupper
    x = last * geom_seq_sum(r,items,MOD)
    x -= B * geom_seq_sum_1(r,items,MOD)
    length = d*items
    S *= pow(10,length,MOD); S += x
    S %= MOD

print(S)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

from functools import lru_cache

L,A,B,MOD = map(int,read().split())

@lru_cache(None)
def F(r,n,MOD):
    
    if n == 0:
        return 0
    q = n//2
    x = F(r,q,MOD)
    
    x *= 1 + pow(r,q,MOD); x %= MOD
    if n&1:
        x *= r; x += 1; x %= MOD
    return x

@lru_cache(None)
def G(r,n,MOD):
    
    if n == 0:
        return 0
    q = n//2
    x = G(r,q,MOD)
    
    x *= 1 + pow(r,q,MOD); x %= MOD
    
    x += F(r,q,MOD) * pow(r,q,MOD) * q % MOD
    if n&1:
        x *= r; x += F(r,n,MOD); x -= 1; x %= MOD
    return x

S = 0
for d in range(1,20):
    
    lower = 10**(d-1)
    upper = 10**d-1
    
    nlower = (lower-A+B-1)//B
    nupper = (upper-A)//B
    if nlower < 0:
        nlower = 0
    if nupper >=L:
        nupper = L-1
    if nlower > nupper:
        continue
    items = nupper - nlower + 1
    r = (10**d) % MOD
    last = A+B*nupper
    x = last * F(r,items,MOD)
    x -= B * G(r,items,MOD)
    length = d*items
    S *= pow(10,length,MOD); S += x
    S %= MOD

print(S)