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):

from collections import defaultdict
from collections import deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS():return list(map(list, sys.stdin.readline().split()))
def S(): return list(sys.stdin.readline())[:-1]
def IR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = I()
    return l
def LIR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LI()
    return l
def SR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = S()
    return l
def LSR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LS()
    return l
sys.setrecursionlimit(1000000)
mod = 1000000007


def A():
    p,q,r = LI()
    print(p+q+r-max(p,q,r))
    return


def B():
    n = I()
    w = LI()
    ans = float("inf")
    for i in range(n):
        m = sum(w[:i])
        l = sum(w[i:])
        ans = min(ans,abs(m-l))
    print(ans)
    return


def C():
    n,m = LI()
    f = [1 for i in range(n+10000)]
    for i in range(m):
        a = I()
        f[a] = 0
    dp = [0 for i in range(n+10000)]
    dp[0] = 1
    for i in range(n):
        if f[i+1]:
            dp[i+1] += dp[i]
            dp[i+1] %= mod
        if f[i+2]:
            dp[i+2] += dp[i]
            dp[i+2] %= mod
    print(dp[n])
    return


def D():
    h,w = LI()
    s = SR(h)
    vi = [[0 for j in range(w)] for i in range(h)]
    ho = [[0 for j in range(w)] for i in range(h)]
    q = deque()
    for y in range(h):
        for x in range(w):
            if s[y][x] == "
                l = len(q)
                while q:
                    i,j = q.pop()
                    vi[i][j] = l
            else:
                q.append((y,x))
        l = len(q)
        while q:
            i,j = q.pop()
            vi[i][j] = l
    for x in range(w):
        for y in range(h):
            if s[y][x] == "
                l = len(q)
                while q:
                    i,j = q.pop()
                    ho[i][j] = l
            else:
                q.append((y,x))
        l = len(q)
        while q:
            i,j = q.pop()
            ho[i][j] = l
    ans = 0
    for y in range(h):
        for x in range(w):
            ans = max(ans,vi[y][x]+ho[y][x]-1)
    print(ans)
    return


def E():
    b = S()
    l = len(b)
    for i in range(l):
        b[i] = int(b[i])
    dp = [[0 for j in range(2)] for i in range(1000002)]
    dp[0][0] = 1
    for i in range(l):
        for j in range(2):
            x = 1 if j else b[i]
            for d in range(x+1):
                if d:
                    dp[i+1][j or d < b[i]] += 2*dp[i][j]
                    dp[i+1][j or d < b[i]] %= mod
                else:
                    dp[i+1][j or d < b[i]] += dp[i][j]
                    dp[i+1][j or d < b[i]] %= mod
    ans = 0
    for j in range(2):
        ans += dp[l][j]
        ans %= mod
    print(ans)
    return


def F():
    def dot(a,b,m):
        return sum(list(map(lambda x:x[0]*x[1]%m, zip(a,b))))%m

    def mat_mul(a,b,m):
        b_ = [[b[i][j] for i in range(len(b))] for j in range(len(b[0]))]
        return [[dot(a[i],b_[j],m) for j in range(len(b_))] for i in range(len(a))]

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

    l,a,b,m = LI()
    ans = 0
    p = 1
    s = 0
    for d in range(1,19):
        left = max(0,-((a-p)//b)) 
        if left >= l:
            break
        p *= 10
        right = min(l,-((a-p)//b))
        if right <= left:
            continue
        mat = mat_pow([[p,1,0],[0,1,1],[0,0,1]],right-left,m)
        ans *= pow(10,d*(right-left),m)
        ans += (mat[0][1]*(a+b*left)%m+mat[0][2]*b%m)%m
        ans %= m
    print(ans)
    return


if __name__ == "__main__":
    F()

Solution (Python):

from collections import defaultdict
from collections import deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS():return list(map(list, sys.stdin.readline().split()))
def S(): return list(sys.stdin.readline())[:-1]
def IR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = I()
    return l
def LIR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LI()
    return l
def SR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = S()
    return l
def LSR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LS()
    return l
sys.setrecursionlimit(1000000)
mod = 1000000007


def A():
    p,q,r = LI()
    print(p+q+r-max(p,q,r))
    return


def B():
    n = I()
    w = LI()
    ans = float("inf")
    for i in range(n):
        m = sum(w[:i])
        l = sum(w[i:])
        ans = min(ans,abs(m-l))
    print(ans)
    return


def C():
    n,m = LI()
    f = [1 for i in range(n+10000)]
    for i in range(m):
        a = I()
        f[a] = 0
    dp = [0 for i in range(n+10000)]
    dp[0] = 1
    for i in range(n):
        if f[i+1]:
            dp[i+1] += dp[i]
            dp[i+1] %= mod
        if f[i+2]:
            dp[i+2] += dp[i]
            dp[i+2] %= mod
    print(dp[n])
    return


def D():
    h,w = LI()
    s = SR(h)
    vi = [[0 for j in range(w)] for i in range(h)]
    ho = [[0 for j in range(w)] for i in range(h)]
    q = deque()
    for y in range(h):
        for x in range(w):
            if s[y][x] == "
                l = len(q)
                while q:
                    i,j = q.pop()
                    vi[i][j] = l
            else:
                q.append((y,x))
        l = len(q)
        while q:
            i,j = q.pop()
            vi[i][j] = l
    for x in range(w):
        for y in range(h):
            if s[y][x] == "
                l = len(q)
                while q:
                    i,j = q.pop()
                    ho[i][j] = l
            else:
                q.append((y,x))
        l = len(q)
        while q:
            i,j = q.pop()
            ho[i][j] = l
    ans = 0
    for y in range(h):
        for x in range(w):
            ans = max(ans,vi[y][x]+ho[y][x]-1)
    print(ans)
    return


def E():
    b = S()
    l = len(b)
    for i in range(l):
        b[i] = int(b[i])
    dp = [[0 for j in range(2)] for i in range(1000002)]
    dp[0][0] = 1
    for i in range(l):
        for j in range(2):
            x = 1 if j else b[i]
            for d in range(x+1):
                if d:
                    dp[i+1][j or d < b[i]] += 2*dp[i][j]
                    dp[i+1][j or d < b[i]] %= mod
                else:
                    dp[i+1][j or d < b[i]] += dp[i][j]
                    dp[i+1][j or d < b[i]] %= mod
    ans = 0
    for j in range(2):
        ans += dp[l][j]
        ans %= mod
    print(ans)
    return


def F():
    def dot(a,b,m):
        return sum(list(map(lambda x:x[0]*x[1]%m, zip(a,b))))%m

    def mat_mul(a,b,m):
        b_ = [[b[i][j] for i in range(len(b))] for j in range(len(b[0]))]
        return [[dot(a[i],b_[j],m) for j in range(len(b_))] for i in range(len(a))]

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

    l,a,b,m = LI()
    ans = 0
    p = 1
    s = 0
    for d in range(1,19):
        left = max(0,-((a-p)//b))
        if left >= l:
            break
        p *= 10
        right = min(l,-((a-p)//b))
        if right <= left:
            continue
        mat = mat_pow([[p,1,0],[0,1,1],[0,0,1]],right-left,m)
        ans *= pow(10,d*(right-left),m)
        ans += (mat[0][1]*(a+b*left)%m+mat[0][2]*b%m)%m
        ans %= m
    print(ans)
    return


if __name__ == "__main__":
    F()

Solution (Python):
import numpy as np
def matrix_power(Mat, n):
    if n == 0:
        return I
    x = matrix_power(Mat, n>>1)
    x = np.dot(x, x)
    if not n&1:
        return x%M
    x = np.dot(x, Mat)
    return x%M
I = np.identity(3, dtype = 'object') 
L, A, B, M = map(int,  input().split())
ans = np.array([[0],
                          [A],
                          [1]], dtype = 'object')
s_L = A+B*(L-1)
for digit in range(len(str(A)), len(str(s_L))+1):
    cnt = (min(s_L, 10**digit-1)-A)//B-max(-1, (10**(digit-1)-1-A)//B)
    Mat = np.array([[10**digit, 1, 0],
                            [0, 1, B],
                            [0, 0, 1]])
    ans = np.dot(matrix_power(Mat, cnt), ans)%M
print(ans[0, 0])

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

def productMatrix(N, A, B):
    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] += A[i][k]*B[k][j]
    return Ret
def modMatrix(N, A, Q): 
    for i in range(N):
        for j in range(N):
            A[i][j] %= Q
    return

def powOfMatrix(N, X, n, Q): 
    Ret = [[1,0,0],[0,1,0],[0,0,1]]
    power = '{:060b}'.format(n)[::-1] 
    for p in power:
        if p == "1":
            Ret = productMatrix(N,Ret, X)
            modMatrix(N, Ret, Q)
        X = productMatrix(N,X,X)
        modMatrix(N, X, Q)
    return Ret

def main():
    L, A, B, M = map( int, input().split())
    s = A
    ANS = [[1,0,0],[0,1,0],[0,0,1]]
    for i in range(1, 37):
        if s >= pow(10,i):
            continue
        P = [[pow(10, i, M), 0, 0], [1,1,0], [0,B,1]]
        step = (pow(10,i)-s+B-1)//B
        if L <= step:
            ANS = productMatrix(3,ANS,powOfMatrix(3,P,L,M))
            modMatrix(3,ANS,M)
            break
        ANS = productMatrix(3,ANS, powOfMatrix(3,P, step,M))
        modMatrix(3,ANS, M)
        L -= step
        s += step*B
    print( (ANS[1][0]*A + ANS[2][0])%M)

if __name__ == '__main__':
    main()

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

def productMatrix(N, A, B):
    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] += A[i][k]*B[k][j]
    return Ret
def modMatrix(N, A, Q): 
    for i in range(N):
        for j in range(N):
            A[i][j] %= Q
    return

def powOfMatrix(N, X, n, Q): 
    Ret = [[1,0,0],[0,1,0],[0,0,1]]
    power = '{:060b}'.format(n)[::-1] 
    for p in power:
        if p == "1":
            Ret = productMatrix(N,Ret, X)
            modMatrix(N, Ret, Q)
        X = productMatrix(N,X,X)
        modMatrix(N, X, Q)
    return Ret

def main():
    L, A, B, M = map( int, input().split())
    s = A
    ANS = [[1,0,0],[0,1,0],[0,0,1]]
    for i in range(1, 37):
        if s >= pow(10,i):
            continue
        P = [[pow(10, i, M), 0, 0], [1,1,0], [0,B,1]]
        step = (pow(10,i)-s+B-1)//B
        if L <= step:
            ANS = productMatrix(3,ANS,powOfMatrix(3,P,L,M))
            modMatrix(3,ANS,M)
            break
        ANS = productMatrix(3,ANS, powOfMatrix(3,P, step,M))
        modMatrix(3,ANS, M)
        L -= step
        s += step*B
    print( (ANS[1][0]*A + ANS[2][0])%M)

if __name__ == '__main__':
    main()

Solution (Python):



def productMatrix(N, A, B):
    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] += A[i][k]*B[k][j]
    return Ret
def modMatrix(N, A, Q): 
    for i in range(N):
        for j in range(N):
            A[i][j] %= Q
    return

def powOfMatrix(N, X, n, Q): 
    Ret = [[1,0,0],[0,1,0],[0,0,1]]
    power = '{:060b}'.format(n)[::-1] 
    for p in power:
        if p == "1":
            Ret = productMatrix(N,Ret, X)
            modMatrix(N, Ret, Q)
        X = productMatrix(N,X,X)
        modMatrix(N, X, Q)
    return Ret

def main():
    L, A, B, M = map( int, input().split())
    s = A
    ANS = [[1,0,0],[0,1,0],[0,0,1]]
    for i in range(1, 37):
        if s >= pow(10,i):
            continue
        P = [[pow(10, i, M), 0, 0], [1,1,0], [0,B,1]]
        step = (pow(10,i)-s+B-1)//B
        if L <= step:
            ANS = productMatrix(3,ANS,powOfMatrix(3,P,L,M))
            modMatrix(3,ANS,M)
            break
        ANS = productMatrix(3,ANS, powOfMatrix(3,P, step,M))
        modMatrix(3,ANS, M)
        L -= step
        s += step*B
    print( (ANS[1][0]*A + ANS[2][0])%M)

if __name__ == '__main__':
    main()

Solution (Python):



def productMatrix(N, A, B):
    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] += A[i][k]*B[k][j]
    return Ret
def modMatrix(N, A, Q): 
    for i in range(N):
        for j in range(N):
            A[i][j] %= Q

    return

def powOfMatrix(N, X, n, Q): 

    Ret = [[1,0,0],[0,1,0],[0,0,1]]
    power = '{:060b}'.format(n)[::-1] 
    for p in power:
        if p == "1":
            Ret = productMatrix(N,Ret, X)

            modMatrix(N, Ret, Q)
        X = productMatrix(N,X,X)

        modMatrix(N, X, Q)

    return Ret

def main():
    L, A, B, M = map( int, input().split())
    s = A
    ANS = [[1,0,0],[0,1,0],[0,0,1]]
    for i in range(1, 37):

        if s >= pow(10,i):
            continue
        P = [[pow(10, i, M), 0, 0], [1,1,0], [0,B,1]]


        step = (pow(10,i)-s+B-1)//B
        if L <= step:
            ANS = productMatrix(3,ANS,powOfMatrix(3,P,L,M))

            modMatrix(3,ANS,M)
            break

        ANS = productMatrix(3,ANS, powOfMatrix(3,P, step,M))


        modMatrix(3,ANS, M)
        L -= step
        s += step*B
    print( (ANS[1][0]*A + ANS[2][0])%M)

        
        
if __name__ == '__main__':
    main()

Solution (Python):
from functools import reduce

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

MOD = M


class ModInt:
    def __init__(self, x):
        self.x = x % MOD

    def __str__(self):
        return str(self.x)

    __repr__ = __str__

    def __add__(self, other):
        return (
            ModInt(self.x + other.x) if isinstance(other, ModInt) else
            ModInt(self.x + other)
        )

    def __sub__(self, other):
        return (
            ModInt(self.x - other.x) if isinstance(other, ModInt) else
            ModInt(self.x - other)
        )

    def __mul__(self, other):
        return (
            ModInt(self.x * other.x) if isinstance(other, ModInt) else
            ModInt(self.x * other)
        )

    def __truediv__(self, other):
        return (
            ModInt(
                self.x * pow(other.x, MOD - 2, MOD)
            ) if isinstance(other, ModInt) else
            ModInt(self.x * pow(other, MOD - 2, MOD))
        )

    def __pow__(self, other):
        return (
            ModInt(
                pow(self.x, other.x, MOD)
            ) if isinstance(other, ModInt) else
            ModInt(pow(self.x, other, MOD))
        )

    def __radd__(self, other):
        return ModInt(other + self.x)

    def __rsub__(self, other):
        return ModInt(other - self.x)

    def __rmul__(self, other):
        return ModInt(other * self.x)

    def __rtruediv__(self, other):
        return ModInt(other * pow(self.x, MOD - 2, MOD))

    def __rpow__(self, other):
        return ModInt(pow(other, self.x, MOD))


def iceil(a, b):
    return (a + b - 1) // b





T = {}
x = A
while x <= A + B * (L - 1):
    k = len(str(x))
    r = iceil(
        min(
            10**k,
            A + B * (L - 1) + 1
        ) - x,
        B
    )
    T[k] = r
    x += r * B


def mut_mul(X, Y):
    return [
        [
            sum(X[i][k] * Y[k][j] for k in range(len(Y)))
            for j in range(len(Y[0]))
        ]
        for i in range(len(X))
    ]


def mut_pow(X, n):
    if n == 0:
        return [
            [1, 0, 0],
            [0, 1, 0],
            [0, 0, 1]
        ]
    else:
        Y = mut_pow(X, n // 2)
        return (
            mut_mul(Y, Y) if n % 2 == 0 else
            mut_mul(mut_mul(Y, Y), X)
        )


def f(p, q, k, r):
    Y = mut_pow(
        [
            [ModInt(10**k), ModInt(1), ModInt(0)],
            [ModInt(0), ModInt(1), ModInt(1)],
            [ModInt(0), ModInt(0), ModInt(1)]
        ],
        r
    )
    Z = mut_mul(
        Y,
        [
            [p],
            [q],
            [ModInt(B)]
        ]
    )
    return (Z[0][0], Z[1][0])


ans, _ = reduce(
    lambda acc, kr: f(acc[0], acc[1], kr[0], kr[1]),
    sorted(T.items()),
    (ModInt(0), ModInt(A))
)

print(ans)

Solution (Python):
def mpw(n,m):
    ret=1
    while m>0:
        if m&1:
            ret*=n
            ret%=M
        n*=n
        n%=M;
        m>>=1
    return ret;
def pw(n,m):
    ret=1
    while m>0:
        if m&1:
            ret*=n
        n*=n
        m>>=1
    return ret;
def delta(i,j):
    if i==j:
        return 1
    return 0
def mtxmpw(A,m):
    ret=[[delta(i,j) for i in range(len(A))] for j in range(len(A))]
    while m>0:
        if m&1:
            ret=mul(ret,A)
        A=mul(A,A)
        m>>=1
    return ret
def mul(A,B):
    ret=[[0 for i in range(len(B[0]))] for j in range(len(A))]
    for i in range(len(A)):
        for j in range(len(B[0])):
            for k in range(len(B)):
                ret[i][j]+=A[i][k]*B[k][j]
                ret[i][j]%=M
    return ret

L,A,B,M=map(int,input().split())
ans=0
d=0
for k in range(18,0,-1):
    l=max((pw(10,(k-1))-A+B-1)//B,0)
    r=min((pw(10,k)-A-1)//B,L-1)
    if l<=r:
        MT=[[1,mpw(10,k),0],[0,mpw(10,k),(M-B*mpw(10,k)%M)%M],[0,0,mpw(10,k)]]
        v=[[(A+B*r)%M],[(A+B*(r-1)%M+M)%M],[1]]
        u=mul(mtxmpw(MT,r-l),v)
        ans+=mpw(10,d)*u[0][0]%M;
        ans%=M;
        d+=(r-l+1)*k
print(ans)

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

def keta(k):
    if A>=10**k:
        return 0,0,0
    begin=10**(k-1)
    end=10**k

    if begin>A+B*(L-1):
        return 0,0,0

    sh=A+max(-1,(begin-1-A)//B)*B+B
    ma=min(A+max(0,(end-1-A)//B)*B,A+B*(L-1))
    kou=(ma-sh)//B+1
    
    return sh,ma,kou

from functools import lru_cache

@lru_cache(maxsize=None)
def f(n,r): 
    if n==1:
        return 1
    if n%2==0:
        return (f(n//2,r) + pow(r,n//2,mod)*f(n//2,r))%mod
    return (r*f(n-1,r)+1)%mod

@lru_cache(maxsize=None)
def g(n,r): 
    if n==1:
        return 0
    if n%2==0:
        return g(n//2,r)+pow(r,n//2,mod)*(g(n//2,r)+n//2*f(n//2,r))
    return r*g(n-1,r)+r*f(n-1,r)

keta_count=1
ANS=0
for i in range(18,0,-1):
    sh,ma,kou = keta(i)
    
    if kou<=0:
        continue
    

    

    ANS=(keta_count*(ma*f(kou,10**i)-B*g(kou,10**i))+ANS)%mod
    
    keta_count = keta_count * pow(10,kou * i,mod)

print(ANS)

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

def keta(k):
    if A>=10**k:
        return 0,0,0
    begin=10**(k-1)
    end=10**k

    if begin>A+B*(L-1):
        return 0,0,0

    sh=A+max(-1,(begin-1-A)//B)*B+B
    ma=min(A+max(0,(end-1-A)//B)*B,A+B*(L-1))
    kou=(ma-sh)//B+1
    
    return sh,ma,kou

from functools import lru_cache

@lru_cache(maxsize=None)
def f(n,r): 
    if n==1:
        return 1
    if n%2==0:
        return (f(n//2,r) + pow(r,n//2,mod)*f(n//2,r))%mod
    return (r*f(n-1,r)+1)%mod

@lru_cache(maxsize=None)
def g(n,r): 
    if n==1:
        return 0
    if n%2==0:
        return g(n//2,r)+pow(r,n//2,mod)*(g(n//2,r)+n//2*f(n//2,r))
    return r*g(n-1,r)+r*f(n-1,r)

keta_count=1
ANS=0
for i in range(18,0,-1):
    sh,ma,kou = keta(i)
    
    if kou<=0:
        continue
    

    

    ANS=(keta_count*(ma*f(kou,10**i)-B*g(kou,10**i))+ANS)%mod
    
    keta_count = keta_count * pow(10,kou * i,mod)

print(ANS)

Solution (Python):
import sys
sys.setrecursionlimit(pow(10, 6))


import numpy as np

def pow_matmul(r, n, m):
  if n == 1:
    return r
  if n % 2 == 0:
    _r = pow_matmul(r, n//2, m)
    return np.mod(np.dot(_r, _r), m)
  else:
    _r = pow_matmul(r, n//2, m)
    _r = np.mod(np.dot(_r, _r), m)
    return np.mod(np.dot(r, _r), m)

def main():
  l, a, b, m = map(int, input().split())
  c = [min(max(0, (10**i - a - 1) // b + 1), l) for i in range(19)]
  ans = 0
  r = np.ones((3,3))
  s = np.array([0, a%m, 1])

  for d in range(1, 19):
    dnum = c[d] - c[d-1]
    if dnum==0:
      continue
    r_s = np.array([[pow(10, d, m), 1, 0], [0, 1, (b%m)], [0, 0, 1]])
    r = pow_matmul(r_s, dnum, m)
    s = np.mod(np.dot(r, s), m)
    if c[d] > l:
      break
  print(s[0])


if __name__ == '__main__':
  main()

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

k = 18

c = [0] * k

t = a + b * (L - 1)

p = 9
temp = 0
for i in range(k):
    if a <= p:
        c[i] = (min(p, t) - a) // b + 1 - temp
        temp += c[i]
    if t <= p:
        break
    p = (p + 1) * 10 - 1

def calc_mat(P, Q):
    ret = [[0,0,0],[0,0,0],[0,0,0]]
    for i in range(3):
        for j in range(3):
            for k in range(3):
                ret[i][j] += P[i][k] * Q[k][j]
            ret[i][j] %= m
    return ret


s = a
x = 0
for i in range(k):
    if c[i] == 0:
        continue
    mat = [[10 ** (i + 1), 0, 0], [1,1,0], [0,b,1]]
    temp_mat = [[1,0,0],[0,1,0],[0,0,1]]
    d = list(reversed(bin(c[i])[2:]))
    p = 1
    for z in d:
        if z == '1':
            temp_mat = calc_mat(temp_mat, mat)
        mat = calc_mat(mat, mat)
    
    temp_x = (temp_mat[0][0] * x + temp_mat[1][0] * s + temp_mat[2][0]) % m
    s = (temp_mat[0][1] * x + temp_mat[1][1] * s + temp_mat[2][1]) % m
    x = temp_x
    

print(x)

Solution (Python):
import sys
sys.setrecursionlimit(1000000)
import numpy as np
input = sys.stdin.readline

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

C = []
for i in range(19):
  C.append(min(max(0,(10**i - A - 1) //B + 1), L))

def R_mul(A, B, M):
  return np.mod((np.dot(A,B)), M)

def R_pow(R, n, M):
  if n == 1:
    return R
  
  R2 = R_mul(R, R, M)
  R2 = R_pow(R2, n//2, M)
  
  if n % 2 == 0:
    return R2
  else:
    return R_mul(R2, R, M)

                 
ans = 0
R = np.ones((3,3))
s = np.array([0, A%M, 1])

for d in range(1,19):
  r = C[d] - C[d-1]
  if (0 == r):
    continue
  R0 = np.array([[pow(10,d,M),1,0],[0,1,(B%M)],[0,0,1]])
  R = R_pow(R0, r, M)
  s = R_mul(R, s, M)
  
  if (C[d] > L):
    break

print(s[0])

Solution (Python):
import sys


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)


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

    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

    return s[0][0]

if __name__ == '__main__':
    print(main())

Solution (Python):
import math
import numpy as np
from decimal import *
L,A,B,M=map(int, input().split())

def matpow(X,n,mod):
  E = np.eye(N=X.shape[0],dtype=int)
  Y = X%mod
  if n == 1:
    return Y
  else:
    return np.dot(matpow(np.dot(Y,Y)%mod,n//2,mod), (Y if n % 2 else E))%mod
  
a = 0
for d in range(1,19):
  left = math.ceil(Decimal(10**(d-1)-A)/Decimal(B))
  right = math.ceil(Decimal(10**d-A)/Decimal(B))-1
  if right < 0:
    continue
  if left > L-1:
    break
  if left > right:
    continue
  if left < 0:
    left = 0
  if right > L-1:
    right = L-1
  mat = np.matrix([[10**d,0,0],[1,1,0],[0,1,1]])
  result = np.dot(np.array([a,A+B*left,B])%M, matpow(mat,right-left+1,M))
  a = result[0,0]%M
  
print(a)

Solution (Python):
import sys
sys.setrecursionlimit(100000000)

def mult(A, B, mod): 
    D = [[(A[0][0] * B[0][0]) % mod, 0, 0], [0, 1, 0], [0, 0, 1]]
    D[1][0] = (A[1][0] * B[0][0] + A[1][1] * B[1][0]) % mod
    D[2][0] = (A[2][0]*B[0][0]+A[2][1]*B[1][0]+A[2][2]*B[2][0]) % mod
    D[2][1] = (A[2][1]*B[1][1]+A[2][2]*B[2][1]) % mod
    return D

def double(C, mod):
    D = [[pow(C[0][0], 2, mod), 0, 0], [(C[1][0] * C[0][0] + C[1][0]) % mod, 1, 0], [(C[2][0] * C[0][0] + C[2][1] * C[1][0] + C[2][0]) % mod, (C[2][1] * 2) % mod, 1]]
    return D

def Cpow(C, n, mod):
    if n == 1: return C
    else:
        D = double(C, mod)
        if n % 2 == 0:
            return Cpow(D, n // 2, mod)
        else:
            return mult(Cpow(D, n // 2, mod), C, mod)

def solve():
    L, a, b, M = map(int, input().split())
    X = 0
    s = a
    nd = -1
    ne = 0
    for d in range(1, 19):
        ne = min((pow(10, d) - a - 1) // b, L-1)
        if ne < 0: continue
        elif ne == nd: continue
        C = [[pow(10, d, M), 0, 0], [1, 1, 0], [0, b % M, 1]]
        D = Cpow(C, ne - nd, M)
        X = (X * D[0][0] + D[1][0] * s + D[2][0]) % M
        s += D[2][1]
        if ne == L - 1: break
        nd = ne
    print(X)
    return

if __name__ == "__main__":
    solve()

Solution (Python):
L, A, B, M = map(int, input().split())
ans = 0
d = 0
def mul(X, Y):
    Z = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    for i in range(3):
        for j in range(3):
            r = 0
            for k in range(3):
                r += X[i][k] * Y[k][j]
            Z[i][j] = r % M
    return Z

x = 0; a = A
while 1:
    kd0 = max(min((10**d - A + B - 1) // B, L), 0)
    kd1 = max(min((10**(d+1) - A + B - 1) // B, L), 0)

    Q = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    P = [[10**(d+1), 0, 0], [1, 1, 0], [0, B, 1]]
    k = kd1 - kd0
    while k:
        if k & 1:
            Q = mul(P, Q)
        P = mul(P, P)
        k >>= 1
    x, a = (x*Q[0][0]+a*Q[1][0]+Q[2][0]) % M, (x*Q[0][1]+a*Q[1][1]+Q[2][1]) % M
    if kd1 == L:
        break
    d += 1
print(x)

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

def mat_pow(x, n) :
    y = [[0] * 3 for _ in range(3)]
        
    for i in range(3) :
        y[i][i] = 1
        
    while n > 0 :
        if n & 1 :
            y = mat_mul(x, y)
        x = mat_mul(x, x)
        n >>= 1
        
    return y
    
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

n0 = -(-(10 ** 0 - A) // B)
n1 = -(-(10 ** 1 - A) // B)

ret = [[0], [A], [1]]

for d in range(1, 19) :
    mat = [[10**d, 1, 0], [0, 1, B], [0, 0, 1]]
    if n0 < 0 and 0 < n1 :
        n0 = 0
        
    if 0 <= n0 < n1 :
        ret = mat_mul(mat_pow(mat, n1 - n0), ret)
    n0, n1 = n1, min(-(-(10 ** (d + 1) - A) // B), L)

print(ret[0][0])

Solution (Python):
import math
import sys
sys.setrecursionlimit(10000)

def nasu(A, B, L, D, M):
  if L == 1:
    return A
  ans = nasuA(A, L, D, M)
  ans = (ans + nasuB(1, L - 1, D, M) * B * D) % M
  return ans % M

def nasuB(B, L, D, M):
  if L == 1:
    return B
  ans = powB(1, L, D, M) % M
  return ans % M

def powB(B, L, D, M):
  if L == 1:
    return B
  k = 0
  t = D
  while T[k + 1] <= L:
    k += 1
    t = t * t % M
    if k + 1 == len(T):
      break
  ans = powB2(B, T[k], D, M) % M
  if L != T[k]:
    ans += nasu(T[k] + 1, 1, L - T[k], D, M) * t % M 
    
    
  return ans % M
 
def powB2(B, L, D, M):
  if L == 1:
    return 1
  ans = powB2(B * 2, L // 2, D * D % M, M) * (D + 1) % M
  ans = ans + nasuA(1, L // 2, D * D % M, M) * B * D % M
  return ans

def nasuA(A, L, D, M):
  if L == 1:
    return A
  ans = powA(A, L, D, M) % M
  return ans % M

def powA(A, L, D, M):
  if L == 1:
    return A
  k = 0
  t = D
  while T[k + 1] <= L:
    k += 1
    t = t * t % M
    if k + 1 == len(T):
      break
  ans = powA2(A, T[k], D, M)
  if L != T[k]:
    ans +=  powA(A, L - T[k], D, M) * t % M
  return ans % M
 
def powA2(A, L, D, M):
  if L == 1:
    return A
  return powA2(A, L // 2, D * D % M, M) * (D + 1) % M

def powmod(a, n, M):
  ans = 1
  while n:
    if n & 1:
      ans = (ans * a) % M
    a = a * a % M
    n >>= 1
  return ans

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

N = math.ceil(math.log(A + 1, 10))
k = pow(10, N)
D = [[0, 0] for _ in range(20)]
while L > 0:
  n = min(L, (k - 1 - A) // B + 1)
  D[N - 1][0] = A
  D[N - 1][1] = n
  L -= n
  N += 1
  k *= 10
  A += n * B
T = [1]
while T[-1] < 10 ** 19:
  T.append(T[-1] * 2)
BI = (((B // M) + 1) * M - B) % M
ans = 0
for i in range(20):
  l = D[i][1]
  a = D[i][0] + (l - 1) * B
  if l == 0:
    continue
  ans = ans * powmod(pow(10, i + 1), l, M) % M
  ans = (ans + nasu(a, BI, l, pow(10, i + 1), M)) % M
print(ans)

Solution (Python):
def mmult(A, B):
    global mod
    n, m, l = len(A), len(B), len(B[0])
    ret = [[0]*l for _ in range(n)]
    for i in range(n):
        for j in range(m):
            for k in range(l):
                ret[i][k] = (ret[i][k]+A[i][j]*B[j][k])%mod
    return ret
def mpow(A, n):
    if n == 0: return [[1 if i == j else 0 for j in range(len(A))] for i in range(len(A))]
    if n % 2: return mmult(mpow(A, n-1), A)
    return mpow(mmult(A, A), n//2)

L, A, B, mod = map(int, input().split())
Y = [max((10**i-A+B-1)//B,0) for i in range(20)]
s = 0
for i in range(19):
    a, b, n = A+B*Y[i], B, max(min(L, Y[i+1])-Y[i], 0)
    M = [[10**(i+1), 0, 0], [1, 1, 0], [0, 1, 1]]
    s = mmult([[s,a,b]], mpow(M,n))[0][0]

print(s)

Solution (Python):
import numpy as np


def exp_np(a, e, m):
    r = np.identity(a.shape[0], dtype=np.int64)
    t = a.copy()
    while e:
        if e % 2:
            r = np.dot(r, t) % m
        t = np.dot(t, t) % m
        e >>= 1
    return r


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

first_d = len(str(a))
tail = a + b * (l - 1)

lo = -1
state = np.array([[0, a % m, 1]])

for d in range(first_d, 19):
    hi = min((10 ** d - 1 - a) // b, l - 1)
    cnt = hi - lo
    coe = np.array([
        [pow(10, d, m), 0, 0],
        [1, 1, 0],
        [0, b % m, 1]
    ], dtype=np.int64)

    coe = exp_np(coe, cnt, m)
    state = np.dot(state, coe) % m

    lo = hi
    if hi == l - 1:
        break

print(state[0][0])