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 rec(r, b, c, a0, n, mod):
    
    a = a0
    n_ = 0
    c += b
    diff = 1
    while n:
        n, m = divmod(n, 2)
        if m:
            a = (a * r + b * n_ + c) % mod
            n_ += diff
        r, b, c = r*r%mod, (r*b+b)%mod, (r*c+b*diff+c)%mod
        diff <<= 1
    return a

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

    def bis(d):
        
        n = 10**(d-1)
        ng = -1
        ok = L
        while ng + 1 < ok:
            c = ok+ng >> 1
            s_c = A + B * c
            if s_c >= n:
                ok = c
            else:
                ng = c
        return ok
    Idxs = [0]
    for d in range(1, 19):
        Idxs.append(bis(d))
    Idxs.append(L)

    ans = 0
    for d, (l, r) in enumerate(zip(Idxs, Idxs[1:])):
        if r - l == 0:
            continue
        s_l = A + B * l
        n = r - l
        ratio = 10**d
        a = rec(ratio, B, s_l, s_l, n-1, M)
        
        ans = (ans * pow(10, d*n, M) + a) % M
    print(ans)

main()

Solution (Python):
L,A,B,mod=map(int,input().split())
def f(l):
  return max(0,(10**l-A-1)//B+1)

def m(a,b):
  r=[[0]*len(b[0]) for i in range(len(a))]
  for i in range(len(a)):
    for k in range(len(b)):
      for j in range(len(b[0])):
        r[i][j]=(r[i][j]+a[i][k]*b[k][j])%mod
  return r

def p(a,n):
  r=[[0]*len(a) for i in range(len(a))]
  b=[]
  for i in range(len(a)):
    r[i][i]=1
    b.append(a[i][:])
  l=n
  while l>0:
    if l&1:
      r=m(b,r)
    b=m(b,b)
    l>>=1
  return r

X=[[0,0,1],[0,0,A],[0,0,0]]
Y,R,Z=0,1,0
while L:
  Y=[[1,0,0],[B,1,0],[0,1,pow(10,R,mod)]]
  Z=min(L,f(R)-f(R-1))
  L-=Z
  Y=p(Y,Z)
  X=m(Y,X)
  R+=1
print(X[2][2]%mod)

Solution (Python):
L,A,B,mod=map(int,input().split())
def f(l):
  return max(0,(10**l-A-1)//B+1)

def m(a,b):
  r=[[0]*len(b[0]) for i in range(len(a))]
  for i in range(len(a)):
    for k in range(len(b)):
      for j in range(len(b[0])):
        r[i][j]=(r[i][j]+a[i][k]*b[k][j])%mod
  return r

def p(a,n):
  r=[[0]*len(a) for i in range(len(a))]
  b=[]
  for i in range(len(a)):
    r[i][i]=1
    b.append(a[i][:])
  l=n
  while l>0:
    if l&1:
      r=m(b,r)
    b=m(b,b)
    l>>=1
  return r

'''
100   001
B10 * 00s
01l   00P
'''

X=[[0,0,1],[0,0,A],[0,0,0]]
Y,R,Z=0,1,0
while L:
  Y=[[1,0,0],[B,1,0],[0,1,pow(10,R,mod)]]
  Z=min(L,f(R)-f(R-1))
  L-=Z
  Y=p(Y,Z)
  X=m(Y,X)
  R+=1
print(X[2][2]%mod)

Solution (Python):
import sys
from math import log10

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 19
MOD = 10 ** 19 + 7
EPS = 10 ** -10

def bisearch_min(mn, mx, func):
    ok = mx
    ng = mn
    while ng+1 < ok:
        mid = (ok+ng) // 2
        if func(mid):
            ok = mid
        else:
            ng = mid
    return ok

def mat_pow(mat, init, K, MOD):
    """ 行列累乗 """

    def mat_dot(A, B, MOD):
        """ 行列の積 """

        if not isinstance(A[0], list) and not isinstance(A[0], tuple):
            A = [A]
        if not isinstance(B[0], list) and not isinstance(A[0], tuple):
            B = [[b] for b in B]
        n1 = len(A)
        n2 = len(A[0])
        _ = len(B)
        m2 = len(B[0])
        res = list2d(n1, m2, 0)
        for i in range(n1):
            for j in range(m2):
                for k in range(n2):
                    res[i][j] += A[i][k] * B[k][j]
                    res[i][j] %= MOD
        return res

    def _mat_pow(mat, k, MOD):
        """ 行列matをk乗する """

        n = len(mat)
        res = list2d(n, n, 0)
        for i in range(n):
            res[i][i] = 1
        while k > 0:
            if k & 1:
                res = mat_dot(res, mat, MOD)
            mat = mat_dot(mat, mat, MOD)
            k >>= 1
        return res

    res = _mat_pow(mat, K, MOD)
    res = mat_dot(res, init, MOD)
    return [a[0] for a in res]

L, a, b, M = MAP()

A = [0] * 20
for i in range(1, 20):
    x = 10 ** i
    
    A[i] = max(min(ceil(x-a, b), L), 0)
C = [0] * 20
for i in range(1, 20):
    C[i] = A[i] - A[i-1]

init = [0, a, 1]
for d in range(1, 20):
    K = C[d]
    if K == 0:
        continue
    mat = [
        
        [pow(10, d, M), 1, 0],
        
        [0, 1, b],
        
        [0, 0, 1],
    ]
    res = mat_pow(mat, init, K, M)
    init[0] = res[0]
    init[1] = res[1]
ans = res[0]
print(ans)

Solution (Python):
import sys
from math import log10

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 19
MOD = 10 ** 19 + 7
EPS = 10 ** -10

def bisearch_min(mn, mx, func):
    ok = mx
    ng = mn
    while ng+1 < ok:
        mid = (ok+ng) // 2
        if func(mid):
            ok = mid
        else:
            ng = mid
    return ok

def mat_pow(mat, init, K, MOD):
    """ 行列累乗 """

    def mat_dot(A, B, MOD):
        """ 行列の積 """

        if not isinstance(A[0], list) and not isinstance(A[0], tuple):
            A = [A]
        if not isinstance(B[0], list) and not isinstance(A[0], tuple):
            B = [[b] for b in B]
        n1 = len(A)
        n2 = len(A[0])
        _ = len(B)
        m2 = len(B[0])
        res = list2d(n1, m2, 0)
        for i in range(n1):
            for j in range(m2):
                for k in range(n2):
                    res[i][j] += A[i][k] * B[k][j]
                    res[i][j] %= MOD
        return res

    def _mat_pow(mat, k, MOD):
        """ 行列matをk乗する """

        n = len(mat)
        res = list2d(n, n, 0)
        for i in range(n):
            res[i][i] = 1
        while k > 0:
            if k & 1:
                res = mat_dot(res, mat, MOD)
            mat = mat_dot(mat, mat, MOD)
            k >>= 1
        return res

    res = _mat_pow(mat, K, MOD)
    res = mat_dot(res, init, MOD)
    return [a[0] for a in res]

L, a, b, M = MAP()

A = [0] * 20
for i in range(1, 20):
    x = 10 ** i
    res = bisearch_min(-1, L, lambda m: ceil(x-a, b) <= m)
    A[i] = res
C = [0] * 20
for i in range(1, 20):
    C[i] = A[i] - A[i-1]

init = [0, a, 1]
for d in range(1, 20):
    K = C[d]
    if K == 0:
        continue
    mat = [
        
        [pow(10, d, M), 1, 0],
        
        [0, 1, b],
        
        [0, 0, 1],
    ]
    res = mat_pow(mat, init, K, M)
    init[0] = res[0]
    init[1] = res[1]
ans = res[0]
print(ans)

Solution (Python):
import sys
from math import log10

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 19
MOD = 10 ** 19 + 7
EPS = 10 ** -10

def bisearch_min(mn, mx, func):
    ok = mx
    ng = mn
    while ng+1 < ok:
        mid = (ok+ng) // 2
        if func(mid):
            ok = mid
        else:
            ng = mid
    return ok

def mat_pow(mat, init, K, MOD):
    """ 行列累乗 """

    def mat_dot(A, B, MOD):
        """ 行列の積 """

        if not isinstance(A[0], list) and not isinstance(A[0], tuple):
            A = [A]
        if not isinstance(B[0], list) and not isinstance(A[0], tuple):
            B = [[b] for b in B]
        n1 = len(A)
        n2 = len(A[0])
        _ = len(B)
        m2 = len(B[0])
        res = list2d(n1, m2, 0)
        for i in range(n1):
            for j in range(m2):
                for k in range(n2):
                    res[i][j] += A[i][k] * B[k][j]
                    res[i][j] %= MOD
        return res

    def _mat_pow(mat, k, MOD):
        """ 行列matをk乗する """

        n = len(mat)
        res = list2d(n, n, 0)
        for i in range(n):
            res[i][i] = 1
        while k > 0:
            if k & 1:
                res = mat_dot(res, mat, MOD)
            mat = mat_dot(mat, mat, MOD)
            k >>= 1
        return res

    res = _mat_pow(mat, K, MOD)
    res = mat_dot(res, init, MOD)
    return [a[0] for a in res]

L, a, b, M = MAP()

A = [0] * 20
for i in range(1, 20):
    x = 10 ** i
    res = bisearch_min(-1, L, lambda m: ceil(x-a, b)-1 < m)
    A[i] = res
C = [0] * 20
for i in range(1, 20):
    C[i] = A[i] - A[i-1]

init = [0, a, 1]
for d in range(1, 20):
    K = C[d]
    if K == 0:
        continue
    mat = [
        [pow(10, d, M), 1, 0],
        [0, 1, b],
        [0, 0, 1],
    ]
    res = mat_pow(mat, init, K, M)
    init[0] = res[0]
    init[1] = res[1]
ans = res[0]
print(ans)

Solution (Python):
import sys
from math import log10

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 19
MOD = 10 ** 19 + 7
EPS = 10 ** -10

def bisearch_min(mn, mx, func):
    ok = mx
    ng = mn
    while ng+1 < ok:
        mid = (ok+ng) // 2
        if func(mid):
            ok = mid
        else:
            ng = mid
    return ok

def mat_pow(mat, init, K, MOD):
    """ 行列累乗 """

    def mat_dot(A, B, MOD):
        """ 行列の積 """

        if not isinstance(A[0], list) and not isinstance(A[0], tuple):
            A = [A]
        if not isinstance(B[0], list) and not isinstance(A[0], tuple):
            B = [[b] for b in B]
        n1 = len(A)
        n2 = len(A[0])
        _ = len(B)
        m2 = len(B[0])
        res = list2d(n1, m2, 0)
        for i in range(n1):
            for j in range(m2):
                for k in range(n2):
                    res[i][j] += A[i][k] * B[k][j]
                    res[i][j] %= MOD
        return res

    def _mat_pow(mat, k, MOD):
        """ 行列matをk乗する """

        n = len(mat)
        res = list2d(n, n, 0)
        for i in range(n):
            res[i][i] = 1
        while k > 0:
            if k & 1:
                res = mat_dot(res, mat, MOD)
            mat = mat_dot(mat, mat, MOD)
            k >>= 1
        return res

    res = _mat_pow(mat, K, MOD)
    res = mat_dot(res, init, MOD)
    return [a[0] for a in res]

L, a, b, M = MAP()

A = [0] * 20
for i in range(1, 20):
    x = 10 ** i
    res = bisearch_min(-1, L, lambda m: ceil(x-a, b)-1 < m)
    A[i] = res
C = [0] * 20
for i in range(1, 20):
    C[i] = A[i] - A[i-1]

init = [0, a, 1]
for d in range(1, 20):
    K = C[d]
    if K == 0:
        continue
    mat = [
        [pow(10, d, M), 1, 0],
        [0, 1, b],
        [0, 0, 1],
    ]
    res = mat_pow(mat, init, K, M)
    init[0] = res[0]
    init[1] = res[1]
ans = res[0]
print(ans)

Solution (Python):
class SquareMatrix():
    def __init__(self, n, mod=1000000007):
        self.n = n
        self.mat = [[0 for j in range(n)] for i in range(n)]
        self.mod = mod

    @staticmethod
    def id(n, mod=1000000007):
        res = SquareMatrix(n, mod)
        for i in range(n):
            res.mat[i][i] = 1
        return res

    @staticmethod
    def modinv(n, mod):
        assert n % mod != 0
        c0, c1 = n, mod
        a0, a1 = 1, 0
        b0, b1 = 0, 1
        while c1:
            a0, a1 = a1, a0 - c0 // c1 * a1
            b0, b1 = b1, b0 - c0 // c1 * b1
            c0, c1 = c1, c0 % c1
        return a0 % mod

    def set(self, arr):
        for i in range(self.n):
            for j in range(self.n):
                self.mat[i][j] = arr[i][j] % self.mod

    def operate(self, vec):
        assert len(vec) == self.n
        res = [0 for _ in range(self.n)]
        for i in range(self.n):
            for j in range(self.n):
                res[i] += self.mat[i][j] * vec[j]
                res[i] %= self.mod
        return res

    def add(self, other):
        assert other.n == self.n
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[i][j] + other.mat[i][j]
                res.mat[i][j] %= self.mod
        return res

    def subtract(self, other):
        assert other.n == self.n
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[i][j] - other.mat[i][j]
                res.mat[i][j] %= self.mod
        return res

    def times(self, k):
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[i][j] * k
                res.mat[i][j] %= self.mod
        return res

    def multiply(self, other):
        assert self.n == other.n
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                for k in range(self.n):
                    res.mat[i][j] += self.mat[i][k] * other.mat[k][j]
                    res.mat[i][j] %= self.mod
        return res

    def power(self, k):
        tmp = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                tmp.mat[i][j] = self.mat[i][j]
        res = SquareMatrix.id(self.n, self.mod)
        while k:
            if k & 1:
                res = res.multiply(tmp)
            tmp = tmp.multiply(tmp)
            k >>= 1
        return res

    def trace(self):
        res = 0
        for i in range(self.n):
            res += self.mat[i][i]
            res %= self.mod
        return res

    def determinant(self):
        res = 1
        tmp = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                tmp.mat[i][j] = self.mat[i][j]
        for j in range(self.n):
            if tmp.mat[j][j] == 0:
                for i in range(j + 1, self.n):
                    if tmp.mat[i][j] != 0:
                        idx = i
                        break
                else:
                    return 0
                for k in range(self.n):
                    tmp.mat[j][k], tmp.mat[idx][k] = tmp.mat[idx][k], tmp.mat[j][k]
                res *= -1
            inv = SquareMatrix.modinv(tmp.mat[j][j], self.mod)
            for i in range(j + 1, self.n):
                c = -inv * tmp.mat[i][j] % self.mod
                for k in range(self.n):
                    tmp.mat[i][k] += c * tmp.mat[j][k]
                    tmp.mat[i][k] %= self.mod
        for i in range(self.n):
            res *= tmp.mat[i][i]
            res %= self.mod
        return res

    def transpose(self):
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[j][i]
        return res

    def inverse(self): 
        res = SquareMatrix.id(self.n, self.mod)
        tmp = SquareMatrix(self.n, self.mod)
        sgn = 1
        for i in range(self.n):
            for j in range(self.n):
                tmp.mat[i][j] = self.mat[i][j]
        for j in range(self.n):
            if tmp.mat[j][j] == 0:
                for i in range(j + 1, self.n):
                    if tmp.mat[i][j] != 0:
                        idx = i
                        break
                else:
                    return 0
                for k in range(self.n):
                    tmp.mat[j][k], tmp.mat[idx][k] = tmp.mat[idx][k], tmp.mat[j][k]
                    res.mat[j][k], res.mat[idx][k] = res.mat[idx][k], res.mat[j][k]
            inv = SquareMatrix.modinv(tmp.mat[j][j], self.mod)
            for k in range(self.n):
                tmp.mat[j][k] *= inv
                tmp.mat[j][k] %= self.mod
                res.mat[j][k] *= inv
                res.mat[j][k] %= self.mod
            for i in range(self.n):
                c = tmp.mat[i][j]
                for k in range(self.n):
                    if i == j:
                        continue
                    tmp.mat[i][k] -= tmp.mat[j][k] * c
                    tmp.mat[i][k] %= self.mod
                    res.mat[i][k] -= res.mat[j][k] * c
                    res.mat[i][k] %= self.mod
        return res

    def linear_equations(self, vec):
        return self.inverse().operate(vec)

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

D = [0 for _ in range(18)]

for i in range(18):
    D[i] = (int('9' * (i + 1)) - A) // B + 1
    D[i] = max(D[i], 0)
    D[i] = min(D[i], L)

for i in range(17)[::-1]:
    D[i + 1] -= D[i]

mat = SquareMatrix.id(3, M)

for i in range(18):
    op = SquareMatrix(3, M)
    op.mat[0][0] = 10**(i + 1)
    op.mat[0][1] = 1
    op.mat[1][1] = 1
    op.mat[1][2] = B
    op.mat[2][2] = 1
    mat = op.power(D[i]).multiply(mat)

print(mat.operate([0, A, 1])[0])

Solution (Python):
def matrix_product(X,Y,M):
    XY = [[0]*3 for _ in range(3)]

    for i in range(3):
        for j in range(3):
            for k in range(3):
                XY[i][j] += X[i][k]*Y[k][j] % M
            XY[i][j] %= M
    
    return XY

def power(X,n,M):
    res = [[1,0,0],[0,1,0],[0,0,1]]

    if n == 0: return res

    while n > 0:
        if n & 1: res = matrix_product(res,X,M)
        X = matrix_product(X,X,M)
        n //= 2
    
    return res

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

X = [0,A,1]

for d in range(1,19):
    if L == 0: break

    
    l = (10**(d-1) - A + B-1) // B
    r = (10**d - A) // B
    if (10**d - A) % B == 0: r -= 1
    
    if l < 0: l = 0
    if r < 0: r = -1
    
    C = min(r-l+1,L)
    L -= C

    Y = power([[10**d,0,0],[1,1,0],[0,B,1]],C,M)
    
    next_X = [0,0,0]
    for i in range(3):
        for j in range(3):
            next_X[i] += X[j] * Y[j][i] % M
        next_X[i] %= M
    X = next_X

print(X[0])

Solution (Python):

def mul_square_matrix(a, b, mod):
    n = len(a)
    ret = [[0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                ret[i][j] = (ret[i][j] + a[i][k] * b[k][j]) % mod
    return ret


def power_square_matrix(mat:list, k, mod):
    n = len(mat)
    ret = [[0] * n for _ in range(n)]
    for i in range(n):
        ret[i][i] = 1
    while k > 0:
        if k % 2 == 1:
            ret = mul_square_matrix(ret, mat, mod)
        mat = mul_square_matrix(mat, mat, mod)
        k //= 2
    return ret

l, a, b, m = map(int, input().split())
s_foot = a + b * (l - 1)
digit_n = len(str(s_foot))
count_digit_n = [0] * (digit_n + 1)
d = 1
for i in range(1, digit_n + 1):
    d *= 10
    x = min(d - 1, s_foot) - a
    if x < 0:
        continue
    n = x // b
    count_digit_n[i] = n + 1

for i in range(digit_n, 0, -1):
    count_digit_n[i]  -= count_digit_n[i - 1]

d = 1
x = [[0], [a], [1]]
for n in count_digit_n[1:]:
    d = d * 10 % m
    if n == 0:
        continue
    mat = [[d, 1, 0], [0, 1, b], [0, 0, 1]]
    mat = power_square_matrix(mat, n, m)
    for i in range(3):
        x[i][0] = sum([mat[i][j] * x[j][0] for j in range(3)]) % m
print(x[0][0])

Solution (Python):


import math
L, A, D, mod = map(int, input().split())
que = []

ans = 0

for i in range(1, 22):
    
    if (10**(i-1)-A+D) % D == 0:
        lowest = max(1, (10**(i-1)-A+D)//D)
    else:
        lowest = max(1, (10**(i-1)-A+D)//D+1)
    highest = 0
    if (10**i-A+D) % D == 0:
        highest = min((10**i-A+D)//D-1, L)

    else:
        highest = min((10**i-A+D)//D, L)

    if 1 <= lowest <= highest:
        que.append((i, lowest, highest))


def poly_sum(X, N, mod):
    
    
    if N == 0:
        return 1
    if N == 1:
        return (X+1) % mod
    else:
        if N % 2 == 0:
            tmp_res = poly_sum(X, N//2, mod)
            return ((pow(X, N//2, mod)+1)*tmp_res-pow(X, N//2, mod)) % mod

        elif N % 2 != 0:
            tmp_res = poly_sum(X, N//2, mod)
            return (pow(X, N//2+1, mod)+1)*tmp_res % mod


def sum_of_poly_sum(X, N, mod):
    
    
    if N == 0:
        return poly_sum(X, 0, mod)
    elif N == 1:
        return (poly_sum(X, 1, mod)+poly_sum(X, 0, mod)) % mod
    else:
        return ((sum_of_poly_sum(X, (N-1)//2, mod))*(pow(X, N//2+1, mod)+1)+(N//2+1)*poly_sum(X, N//2, mod)) % mod


def find_mod(length, lowest, highest, mod):
    
    X = pow(10, length, mod)
    if highest > lowest:
        base = (A+D*(lowest-1))*poly_sum(X, highest-lowest, mod) + \
            D*sum_of_poly_sum(X, highest-lowest - 1, mod)
    elif highest == lowest:
        base = A+D*(lowest-1)
    return base % mod


que.reverse()

digits = 0
for length, lowest, highest in que:
    ans += find_mod(length, lowest, highest, mod)*pow(10, digits, mod)
    ans %= mod
    digits += length*(highest+1-lowest)

print(ans % mod)

Solution (Python):


import math
L, A, D, mod = map(int, input().split())
que = []

ans = 0

for i in range(1, 22):
    
    if (10**(i-1)-A+D) % D == 0:
        lowest = max(1, (10**(i-1)-A+D)//D)
    else:
        lowest = max(1, (10**(i-1)-A+D)//D+1)
    highest = 0
    if (10**i-A+D) % D == 0:
        highest = min((10**i-A+D)//D-1, L)

    else:
        highest = min((10**i-A+D)//D, L)

    if 1 <= lowest <= highest:
        que.append((i, lowest, highest))


def poly_sum(X, N, mod):
    
    
    if N == 0:
        return 1
    if N == 1:
        return (X+1) % mod
    else:
        if N % 2 == 0:
            tmp_res = poly_sum(X, N//2, mod)
            return ((pow(X, N//2, mod)+1)*tmp_res-pow(X, N//2, mod)) % mod

        elif N % 2 != 0:
            tmp_res = poly_sum(X, N//2, mod)
            return (pow(X, N//2+1, mod)+1)*tmp_res % mod


def sum_of_poly_sum(X, N, mod):
    
    
    if N == 0:
        return poly_sum(X, 0, mod)
    elif N == 1:
        return (poly_sum(X, 1, mod)+poly_sum(X, 0, mod)) % mod
    else:
        return ((sum_of_poly_sum(X, (N-1)//2, mod))*(pow(X, N//2+1, mod)+1)+(N//2+1)*poly_sum(X, N//2, mod)) % mod


def find_mod(length, lowest, highest, mod):
    
    X = pow(10, length, mod)
    if highest > lowest:
        base = (A+D*(lowest-1))*poly_sum(X, highest-lowest, mod) + \
            D*sum_of_poly_sum(X, highest-lowest - 1, mod)
    elif highest == lowest:
        base = A+D*(lowest-1)
    return base % mod


que.reverse()

digits = 0
for length, lowest, highest in que:
    ans += find_mod(length, lowest, highest, mod)*pow(10, digits, mod)
    ans %= mod
    digits += length*(highest+1-lowest)

print(ans % mod)

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], [b]]
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, 1], [0, 0, 1]]
    ret = mat_mul(mat_pow(mat, d1 - d0), ret)
    if d1 == l:
        break
    d0 = d1


print(ret[0][0])

Solution (Python):
import sys
sys.setrecursionlimit(10**7) 
import math
from copy import copy, deepcopy
from copy import deepcopy as dcp
from operator import itemgetter
from bisect import bisect_left, bisect, bisect_right

from collections import deque, defaultdict


from collections import Counter

from itertools import accumulate,combinations,permutations

from heapq import heapify,heappop,heappush


from functools import reduce,lru_cache

from decimal import Decimal

def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printe(*x):print("
def printl(li): _=print(*li, sep="\n") if li else None
def argsort(s, return_sorted=False): 
    inds=sorted(range(len(s)), key=lambda k: s[k])
    if return_sorted: return inds, [s[i] for i in inds]
    return inds
def alp2num(c,cap=False): return ord(c)-97 if not cap else ord(c)-65
def num2alp(i,cap=False): return chr(i+97) if not cap else chr(i+65)
def matmat(A,B,mod):
    K,N,M=len(B),len(A),len(B[0])
    return [[sum([(A[i][k]*B[k][j]) for k in range(K)])%mod for j in range(M)] for i in range(N)]
def matvec(M,v):
    N,size=len(v),len(M)
    return [sum([M[i][j]*v[j] for j in range(N)]) for i in range(size)]
def T(M):
    n,m=len(M),len(M[0])
    return [[M[j][i] for j in range(n)] for i in range(m)]
def main():

    

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

    mod=M
    ans=0
    cur=0
    a10=1
    cs=A
    for alp in range(18):
        l=a10-A-1
        r=a10*10-A-1
        a10*=10

        if r<0:
            
            continue
        if l<0:
            lb=0
        else:
            lb=l//B+1
        rb=r//B+1
        
        rrb=min(L,rb)
        n=(rrb-lb)
        mat=[[0]*3 for _ in range(3)]
        mat[0][0]=pow(10,alp+1,mod)
        mat[0][1]=1
        mat[1][1]=1
        mat[1][2]=B%mod
        mat[2][2]=1
        v=[ans,cs,1]
        tot=[[0]*3 for _ in range(3)]
        for i in range(3):tot[i][i]=1
        
        for j in range(n.bit_length()):
            if (n>>j)&1:
                tot=matmat(mat,tot,M)
            mat=matmat(mat,mat,M)
            

        x=matvec(tot,v)
        ans=x[0]%mod
        cs=x[1]%mod
        
        if rb>L:
            break
        


    print(ans)



if __name__ == "__main__":
    main()

Solution (Python):
import math, itertools

L, A, B, mod = [int(_) for _ in input().split()]


def sum_ri(r, n, mod):
    
    if r == 1:
        return n % mod
    else:
        
        
        return (pow(r, n, mod * (r - 1)) - 1) // (r - 1) % mod


def sum_iri(r, n, mod):
    
    if r == 1:
        return n * (n - 1) // 2 % mod
    else:
        r1 = r - 1
        p = pow(r, n, mod * r1**2)
        ret = (n - 1) * p * r1 + r1 - (p - 1)
        ret //= r1**2
        ret %= mod
        return ret


ans = 0

for digit in range(1, 20):
    
    
    
    ileft = max(0, (10**(digit - 1) - A - 1) // B + 1)
    iright = min(L - 1, (10**digit - A - 1) // B)
    if L <= ileft:
        break
    if ileft > iright:
        continue
    idiff = iright - ileft + 1
    
    
    
    pow10digit = pow(10, digit, mod)
    now = (A + B * iright) * sum_ri(pow10digit, idiff, mod)
    now -= B * sum_iri(pow10digit, idiff, mod)
    ans *= pow(10, digit * idiff, mod)
    ans += now
    ans %= mod
print(ans)