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

  def __init__(self, mat, mod=10**9+7):
    assert(len(mat) and len(mat[0]))
    from copy import deepcopy
    self.MOD = mod
    self.mat = deepcopy(mat)

  def matunit(self, n):
    return Matrix([[(x==y) for x in range(n)] for y in range(n)], mod=self.MOD)

  def __add__(self, other):
    assert(len(self.mat)==len(other.mat) and len(self.mat[0])==len(other.mat[0]))
    h = len(self.mat)
    w = len(self.mat[0])
    ret = [[(self.mat[y][x]+other.mat[y][x])%self.MOD for x in range(w)] for y in range(h)]
    return ret

  def __mul__(self, other):
    assert(len(self.mat[0])==len(other.mat))
    h = len(self.mat)
    w = len(other.mat[0])
    ret = [[0 for _ in range(w)] for _ in range(h)]
    for y in range(h):
      for z in range(len(other.mat)):
        for x in range(w):
          ret[y][x] += self.mat[y][z]*other.mat[z][x]
    for y in range(h):
      for x in range(w):
        ret[y][x] %= self.MOD
    return Matrix(ret, mod=self.MOD)

  def __pow__(self, k):
    assert(len(self.mat)==len(self.mat[0]))
    from copy import deepcopy
    ret = self.matunit(len(self.mat))
    n = Matrix(deepcopy(self.mat), mod=self.MOD)
    while k:
      if k&1:
        ret = ret * n
      n = n * n
      k >>= 1
    return ret

  def __str__(self):
    return "[{}]".format("\n".join(str(row) for row in self.mat))

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

d = [0 for i in range(20)]
for i in range(1, 19):
  l = -1
  r = L
  while r-l>1:
    m = (l+r)//2
    if len(str(A+B*m))>i:
      r = m
    else:
      l = m
  d[i] = r

x = Matrix([[0, A, 1]], mod=M)
for i in range(1, 19):
  dt = d[i]-d[i-1]
  y = Matrix([[10**i, 0, 0], [1, 1, 0], [0, B, 1]], mod=M)
  y = y ** dt
  x *= y
print(x.mat[0][0])

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






def c111(n, l, m):
    if n <= 0:
        return 0
    if n == 1:
        return 1
    if n % 2 == 1:
        return (c111(n - 1, l, m) * pow(10, l, m) + 1) % m
    half = c111(n // 2, l, m)
    return (half * pow(10, (n // 2) * l, m) + half) % m






def c123(n, l, m):
    if n <= 0:
        return 0
    if n == 1:
        return 1
    if n % 2 == 1:
        return (c123(n - 1, l, m) + c111(n, l, m)) % m
    half = c123(n // 2, l, m)
    return (half * pow(10, (n // 2) * l, m) + half + (n // 2) * c111(n // 2, l, m)) % m

fst, lst = a, a + b * (l - 1)

fst_l, lst_l = len(str(fst)), len(str(lst))

res, margin = 0, 0
for keta in reversed(range(fst_l, lst_l + 1)):
    num_l = a + b * ((10 ** (keta - 1) - a + b - 1) // b)
    num_r = a + b * ((10 ** keta - a + b - 1) // b - 1)
    if keta == fst_l:
        num_l = fst
    if keta == lst_l:
        num_r = lst
    if num_l > num_r:
        continue
    sz = (num_r - num_l) // b + 1
    _111 = num_l * c111(sz, keta, m)
    _123 = b * c123(sz - 1, keta, m)
    res += (pow(10, margin, m) * (_111 + _123)) % m
    margin += sz * keta
print(res % m)

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







def get_indices(d):
	i_start = (10**d - A + (B-1)) / B
	i_end = (10**(d+1) - A + (B-1)) / B
	i_start = max(0, i_start)
	i_end = max(0, i_end)
	i_start = min(L, i_start)
	i_end = min(L, i_end)
	return (d+1, i_start, i_end)
	
index_pairs = [get_indices(n) for n in xrange(0, 19)]

def get_chunk_value(num_digs, i_start, i_end):
	if i_start == i_end:
		return (0, 0)
	
	a = A + B*(i_end - 1)
	d = -B
	n = i_end - i_start
	r = 10 ** num_digs
	
	Mt = M*((1-r)**2)
	
	numerator = (1-r)*(a - (a + n*d)*pow(r, n, Mt)) + d*r*(1 - pow(r, n, Mt))
	numerator %= Mt
	assert numerator % (1-r)**2 == 0
	return (num_digs*(i_end-i_start), numerator/(1-r)**2)
	
	


chunk_values = [get_chunk_value(n, i_start, i_end) for n, i_start, i_end 
													in index_pairs]


ans = 0
digits_offset = 0
for num_digs, value in reversed(chunk_values):
	ans += value * pow(10, digits_offset, M)
	digits_offset += num_digs
	
ans %= M
print ans

Solution (Python):
L,A,B,M=map(int,input().split())
C=(L-1)*B+A
mod=M
X,s=0,A
now=A
data=[0]*19
ruiseki=[0]*19
Ami=(A-1)//B
Cmi=C//B
for i in range(19):
    now=10**i-1
    if now<A:
        ruiseki[i]=0
    elif now>=C:
        ruiseki[i]=L
    else:
        ruiseki[i]=(now-A)//B+1
def mul33(A,B):
    global mod
    res=[[0]*3 for i in range(3)]
    for i in range(3):
        for j in range(3):
            res[i][j]=sum([A[i][k]*B[k][j]%mod for k in range(3)])
    return res
def ruijou(vec,Cd):
    rui=dict()
    if Cd==0:
        return [[1,0,0],[0,1,0],[0,0,1]]
    now=vec
    ans=[[1,0,0],[0,1,0],[0,0,1]]
    Cds=bin(Cd)[2:][::-1]
    for i in range(len(Cds)):
        if Cds[i]=="1":
            ans=mul33(ans,now)
        now=mul33(now,now)
    return ans
def mul(X,s,Cd,d):
    global mod
    global B
    vec=[[pow(10,d,mod),0,0],[1,1,0],[0,B,1]]
    vec=ruijou(vec,Cd)
    Xc=X
    sc=s
    X=(Xc*vec[0][0]+sc*vec[1][0]+vec[2][0])%mod
    s=(Xc*vec[0][1]+sc*vec[1][1]+vec[2][1])%mod
    return X,s
for d in range(1,19):
    Cd=ruiseki[d]-ruiseki[d-1]
    X,s=mul(X,s,Cd,d)
print(X)

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






def c111(n, l, m):
    if n <= 1:
        return 1
    if n % 2 == 1:
        return (c111(n - 1, l, m) * pow(10, l, m) + 1) % m
    half = c111(n // 2, l, m)
    return (half * pow(10, (n // 2) * l, m) + half) % m






def c123(n, l, m):
    if n <= 1:
        return 1
    if n % 2 == 1:
        return (c123(n - 1, l, m) + c111(n, l, m)) % m
    half = c123(n // 2, l, m)
    return (half * pow(10, (n // 2) * l, m) + half + (n // 2) * c111(n // 2, l, m)) % m

fst = a
lst = a + b * (l - 1)

fst_l = len(str(fst))
lst_l = len(str(lst))

res = 0
margin = 0
for keta in reversed(range(fst_l, lst_l + 1)):
    num_l = a + b * ((10 ** (keta - 1) - a + b - 1) // b)
    num_r = a + b * ((10 ** keta - a + b - 1) // b - 1)
    if keta == fst_l:
        num_l = fst
    if keta == lst_l:
        num_r = lst
    if num_l > num_r:
        continue
    sz = (num_r - num_l) // b + 1
    _111 = num_l * c111(sz, keta, m)
    _123 = b * c123(sz - 1, keta, m)
    if sz == 1:
        _123 = 0
    res += (pow(10, margin, m) * (_111 + _123)) % m
    margin += sz * keta
print(res % m)

Solution (Python):
import os
import sys

import numpy as np
from numpy.linalg import matrix_power

if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")

sys.setrecursionlimit(2147483647)
INF = float("inf")

L, A, B, M = list(map(int, sys.stdin.readline().split()))
MOD = M


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

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

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

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


size = 0


vec = [ModInt(0), ModInt(A), ModInt(B)]

digits = 1
while size < L:
    
    
    cnt = max(0, (10 ** digits - 1 - A) // B - size + 1)
    cnt = min(cnt, L - size)

    mat = np.array([
        [ModInt(pow(10, digits, MOD)), ModInt(0), ModInt(0)],
        [ModInt(1), ModInt(1), ModInt(0)],
        [ModInt(0), ModInt(1), ModInt(1)],
    ])
    if cnt > 0:
        vec = np.dot(vec, matrix_power(mat, cnt))
    digits += 1
    size += cnt
print(vec[0])

Solution (Python):
import copy
import sys

stdin = sys.stdin

ni = lambda: int(ns())
na = lambda: list(map(int, stdin.readline().split()))
ns = lambda: stdin.readline().rstrip()  

L,A,B,mod = na()

low = 1
high = 10


def matpow(M, v, e, mod):
    A = copy.deepcopy(M)
    w = copy.deepcopy(v)
    while e > 0:
        if e&1:
            w = mulv(A, w, mod)
        A = mul(A, A, mod)
        e >>= 1
    return w


def mulv(M, v, mod):
    n = len(M)
    m = len(v)
    ret = [0] * n
    for i in range(n):
        s = 0
        for j in range(m):
            s += M[i][j] * v[j]
        ret[i] = s % mod
    return ret


def mul(A, B, mod):
    n = len(A)
    m = len(B)
    o = len(B[0])
    ret = [[0] * o for _ in range(n)]
    for i in range(n):
        for j in range(o):
            s = 0
            for k in range(m):
                s += A[i][k] * B[k][j]
            ret[i][j] = s % mod
    return ret








v = [0, A, B]
ra = A

while low < 1e18:
    mat = [[high%mod, 1, 0], [0, 1, 1], [0, 0, 1]]
    step = max(0, min(L, (high-ra+B-1)//B))
    v = matpow(mat, v, step, mod)
    
    ra = ra + B * step
    L -= step

    low *= 10
    high *= 10

print(v[0])

Solution (Python):
l, a, b, m = [int(i) for i in input().split()]

cd = []
for i in range(18):
    cd.append(max(0, (10 ** (i + 1) - 1 - a) // b + 1))
    if cd[-1] >= l:
        cd[-1] = l
        break

cd_sum = 0
cd_2 = []

for i in cd:
    cd_2.append(i-cd_sum)
    cd_sum += i-cd_sum

X = [0, a, 1]
B = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
for i in range(len(cd_2)):
    A = [[(10 ** (i + 1)), 0, 0], [1, 1, 0], [0, b, 1]]
    B = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

    n = cd_2[i]
    while n > 0:
        if n % 2 == 1:
            B = [[(B[0][0] * A[0][0] + B[0][1] * A[1][0] + B[0][2] * A[2][0]) % m,
                  (B[0][0] * A[0][1] + B[0][1] * A[1][1] + B[0][2] * A[2][1]) % m,
                  (B[0][0] * A[0][2] + B[0][1] * A[1][2] + B[0][2] * A[2][2]) % m],
                 [(B[1][0] * A[0][0] + B[1][1] * A[1][0] + B[1][2] * A[2][0]) % m,
                  (B[1][0] * A[0][1] + B[1][1] * A[1][1] + B[1][2] * A[2][1]) % m,
                  (B[1][0] * A[0][2] + B[1][1] * A[1][2] + B[1][2] * A[2][2]) % m],
                  [(B[2][0] * A[0][0] + B[2][1] * A[1][0] + B[2][2] * A[2][0]) % m,
                   (B[2][0] * A[0][1] + B[2][1] * A[1][1] + B[2][2] * A[2][1]) % m,
                   (B[2][0] * A[0][2] + B[2][1] * A[1][2] + B[2][2] * A[2][2]) % m]]

            n -= 1

        else:
            A = [[(A[0][0] * A[0][0] + A[0][1] * A[1][0] + A[0][2] * A[2][0]) % m,
                  (A[0][0] * A[0][1] + A[0][1] * A[1][1] + A[0][2] * A[2][1]) % m,
                  (A[0][0] * A[0][2] + A[0][1] * A[1][2] + A[0][2] * A[2][2]) % m],
                 [(A[1][0] * A[0][0] + A[1][1] * A[1][0] + A[1][2] * A[2][0]) % m,
                  (A[1][0] * A[0][1] + A[1][1] * A[1][1] + A[1][2] * A[2][1]) % m,
                  (A[1][0] * A[0][2] + A[1][1] * A[1][2] + A[1][2] * A[2][2]) % m],
                  [(A[2][0] * A[0][0] + A[2][1] * A[1][0] + A[2][2] * A[2][0]) % m,
                   (A[2][0] * A[0][1] + A[2][1] * A[1][1] + A[2][2] * A[2][1]) % m,
                   (A[2][0] * A[0][2] + A[2][1] * A[1][2] + A[2][2] * A[2][2]) % m]]
            n //= 2

    X[0] = (X[0] * B[0][0] + X[1] * B[1][0] + X[2] * B[2][0]) % m
    X[1] = (X[0] * B[0][1] + X[1] * B[1][1] + X[2] * B[2][1]) % m
    X[2] = (X[0] * B[0][2] + X[1] * B[1][2] + X[2] * B[2][2]) % m

print(X[0] % m)

Solution (Python):
L,A,B,M=map(int,input().split())
K=19
gap=[-1 for k in range(K)]
for k in range(K):
    tmp=10**k-A
    if tmp%B==0:
        gap[k]=-1+tmp//B
    else:
        gap[k]=tmp//B
    gap[k]=min(L-1,gap[k])
    gap[k]=max(-1,gap[k])
num=[0 for k in range(K)]
for k in range(1,K):
    
    num[k]=-gap[k-1]+gap[k]
def matrix(d):
    return [[pow(10,d,M),1,0],[0,1,B],[0,0,1]]
def mattime(P,Q):
    R=[[0,0,0],[0,0,0],[0,0,0]]
    for i in range(3):
        for j in range(3):
            for k in range(3):
                R[i][j]+=P[i][k]*Q[k][j]
                R[i][j]%=M
    return R
def matexp(P,n):
    if n==0:
        return [[1,0,0],[0,1,0],[0,0,1]]
    elif n%2==0:
        return matexp(mattime(P,P),n//2)
    else:
        return mattime(P,matexp(mattime(P,P),n//2))
Mat=[[1,0,0],[0,1,0],[0,0,1]]
for d in range(1,K):
    Mat=mattime(matexp(matrix(d),num[d]),Mat)
print((A*Mat[0][1]+Mat[0][2])%M)

Solution (Python):
L,A,B,M=map(int,input().split())
K=19
gap=[-1 for k in range(K)]
for k in range(K):
    tmp=10**k-A
    if tmp%B==0:
        gap[k]=-1+tmp//B
    else:
        gap[k]=tmp//B
    gap[k]=min(L-1,gap[k])
    gap[k]=max(-1,gap[k])
num=[0 for k in range(K)]
for k in range(1,K):
    
    num[k]=-gap[k-1]+gap[k]
def matrix(d):
    return [[pow(10,d,M),1,0],[0,1,B],[0,0,1]]
def mattime(P,Q):
    R=[[0,0,0],[0,0,0],[0,0,0]]
    for i in range(3):
        for j in range(3):
            for k in range(3):
                R[i][j]+=P[i][k]*Q[k][j]
                R[i][j]%=M
    return R
def matexp(P,n):
    if n==0:
        return [[1,0,0],[0,1,0],[0,0,1]]
    elif n%2==0:
        return matexp(mattime(P,P),n//2)
    else:
        return mattime(P,matexp(mattime(P,P),n//2))
Mat=[[1,0,0],[0,1,0],[0,0,1]]
for d in range(1,K):
    Mat=mattime(matexp(matrix(d),num[d]),Mat)
print((A*Mat[0][1]+Mat[0][2])%M)

Solution (Python):
def f_takahashi_basics_in_education_and_learning(L, A, B, M, DIGIT_MAX=18):
    import numpy

    
    def matrix_power(a, n):
        a = numpy.core.asanyarray(a)
        fmatmul = numpy.core.dot

        z = result = None
        while n > 0:
            z = a % M if z is None else fmatmul(z, z) % M
            n, bit = divmod(n, 2)
            if bit:
                result = z if result is None else fmatmul(result, z) % M
        return result

    
    num_pow_d = [0] * (DIGIT_MAX + 1)
    for d in range(DIGIT_MAX + 1):
        if 10**d - 1 >= A:
            
            num_pow_d[d] = min(((10**d - 1 - A) // B) + 1, L)  

    ans = numpy.array([0, A, 1])
    for d in range(1, DIGIT_MAX + 1):
        c_d = num_pow_d[d] - num_pow_d[d - 1]
        matrix = numpy.array([[(10**d) % M, 0, 0], [1, 1, 0], [0, B, 1]], dtype='object')
        m = numpy.eye(3, dtype='object') if c_d == 0 else matrix_power(matrix, c_d)
        ans = ans.dot(m) % M
    return ans[0]

L, A, B, M = [int(i) for i in input().split()]
print(f_takahashi_basics_in_education_and_learning(L, A, B, M))

Solution (Python):
l, a, b, MOD = [int(item) for item in input().split()]
vec = [0, a%MOD, 1]

def mat_mat_mul3(m1, m2):
    m3 = [[0] * 3 for _ in range(3)]
    for i in range(3):
        for j in range(3):
            for k in range(3):
                m3[i][j] += m1[i][k] * m2[k][j]
                m3[i][j] %= MOD
    return m3

def mat_vec_mul3(m1, v1):
    v2 = [0] * 3
    for i in range(3):
        for k in range(3):
            v2[i] += v1[k] * m1[i][k]
            v2[i] %= MOD
    return v2

A_d_pow = [[None] * 66 for _ in range(20)]
for d in range(20):
    for i in range(66):
        if i == 0:
            A = [[1, 0, 0], 
                 [0, 1, 0], 
                 [0, 0, 1]]
            A_d_pow[d][i] = A
            continue
        if i == 1:
            A = [[10**d % MOD, 1, 0], 
                 [0, 1, b % MOD], 
                 [0, 0, 1]]
            A_d_pow[d][i] = A
            continue
        A_d_pow[d][i] = mat_mat_mul3(A_d_pow[d][i-1], A_d_pow[d][i-1])

multimes = [0] * 20
for i in range(1,20):
    if 10**i > a:
        cnt = (10**i - 1 - a) // b + 1
        if cnt >= l:
            multimes[i] = l
            break
        else:
            multimes[i] = cnt
            l -= cnt
            a += cnt * b

for i, num in enumerate(multimes):
    if num == 0:
        continue
    for j in range(66):
        if num & 1 << j:
            vec = mat_vec_mul3(A_d_pow[i][j+1], vec)

print(vec[0])

Solution (Python):
n,a,b,MOD=map(int,input().split())
def matmul(A,B):
	global MOD
	N,M,L=len(A),len(B),len(B[0])
	ret=[[0 for j in range(L)] for i in range(N)]
	for i in range(N):
		for j in range(M):
			for k in range(L):
				ret[i][k] += A[i][j]*B[j][k]%MOD
				ret[i][k] %= MOD
	return ret
ans,c=[[0,a%MOD,1]],0
for i in range(1,19):
	if not a<10**i:
		continue
	m=[[10**i%MOD,0,0],[1,1,0],[0,b%MOD,1]]
	q = -c
	c = min(((10**i-1-a)//b)+1,n)
	q += c
	while q:
		if q&1:
			ans = matmul(ans,m)
		m = matmul(m,m)
		q >>= 1
print(ans[0][0])

Solution (Python):
def f_takahashi_basics_in_education_and_learning(L, A, B, M, DIGIT_MAX=18):
    """numpyの行列演算だと L, A, B が 10**18 程度のときうまく動作しない"""
    def identity_matrix(n):
        """n×n単位行列"""
        return [[int(i == j) for i in range(n)] for j in range(n)]

    def multiply(a, b):
        """行列a, bの積"""
        row = len(a)
        col = len(b[0])
        t = len(b)
        ret = [[0] * col for _ in range(row)]
        for i in range(row):
            for j in range(col):
                for k in range(t):
                    ret[i][j] += a[i][k] * b[k][j]
        for i in range(row):
            for j in range(col):
                ret[i][j] %= M
        return ret

    def power(a, k):
        """行列aのk乗"""
        n = len(a)
        ret = identity_matrix(n)
        while k:
            if k & 1:
                ret = multiply(ret, a)
            a = multiply(a, a)
            k >>= 1
        return ret
 
    
    num_pow_d = [0] * (DIGIT_MAX + 1)
    for d in range(DIGIT_MAX + 1):
        if 10**d - 1 < A:
            continue
        
        num_pow_d[d] = min(((10**d - 1 - A) // B) + 1, L)  
    ans = 0
    array = [[0, A, 1]]
    for d in range(1, DIGIT_MAX + 1):
        c_d = num_pow_d[d] - num_pow_d[d - 1]
        matrix = [[10**d, 0, 0], [1, 1, 0], [0, B, 1]]
        array = multiply(array, power(matrix, c_d))
    return array[0][0]

L, A, B, M = [int(i) for i in input().split()]
print(f_takahashi_basics_in_education_and_learning(L, A, B, M))

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

ten = 0

def pow_(x, t):
  if t == 0:
    return 1
  if t % 2 == 1:
    return (pow_(x, t - 1) * x) % M
  else:
    y = pow_(x, t // 2)
    return (y * y) % M

def f(l):
  if l == 0:
    return 0
  if l % 2 == 1:
    pl = l - 1
    x = f(pl)
    return x * ten + 1
  else:
    pl = l // 2
    x = f(pl)
    return x * pow_(ten, pl) + x

def g(l):
  if l == 0:
    return 0
  if l % 2 == 1:
    pl = l - 1
    x = g(pl)
    return x * ten + B * pl
  else:
    pl = l // 2
    x = g(pl)
    return x * pow_(ten, pl) + x + B * pl * f(pl)

last = A + B * (L - 1)
ans = 0
ten = 10
for _ in range(18):
  l = ten // 10
  r = ten - 1
  
  if last < l:
    ten *= 10
    continue
  if A > r:
    ten *= 10
    continue
  na = 0 
  nl = 0 

  if A >= l:
    na = A
  else:
    na = (l - A + B - 1) // B * B + A
    na = min(na, last)
  
  nlast = 0
  if last <= r:
    nlast = last
  else:
    nlast = (r - A) // B * B + A
  nl = (nlast - na) // B + 1

  ans *= pow_(ten, nl) % M
  ans += (na * f(nl)) % M
  ans += g(nl) % M
  ten *= 10

print(ans % M)

Solution (Python):
def mat_dot(A, B, mod):
    assert len(A[0]) == len(B), 'invalid_size'
    
    L = len(A)
    M = len(A[0])
    N = len(B[0])
    
    res = [[0]*N for _ in range(L)]
    
    for i in range(L):
        for j in range(N):
            a = 0
            for k in range(M):
                a = (a+A[i][k]*B[k][j]) % mod
            res[i][j] = a
            
    return res

def mat_pow(A, x, mod):
    N = len(A)
    res = [[0]*N for _ in range(N)]
    for i in range(N):
        res[i][i] = 1
    for i in range(x.bit_length()):
        if 2**i & x:
            res = mat_dot(res, A, mod)
        A = mat_dot(A, A, mod)
    return res
      
L, A, B, M = map(int, input().split())

pre = -1
S = []
for k in range(50):
    if len(str(A)) > k+1:
        continue
    j = -((A - 10**(k+1))//B) - 1
    j = min(j, L-1)
    j = j - pre
    pre += j
    if j:
        S.append((j, k+1))

Ans = [[0], [A], [1]]
for p, k in S:
    cmat = [[pow(10, k, M), 1, 0], [0, 1, B], [0, 0, 1]]
    
    Ans = mat_dot(mat_pow(cmat, p, M), Ans, M)
print(Ans[0][0])

Solution (Python):



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

def matpow(A, p, mod): 
    if p % 2 == 1:
        return matmul(A, matpow(A, p - 1, mod), mod)
    elif p > 0:
        b = matpow(A, p // 2, mod)
        return matmul(b, b, mod)
    else:
        return [[1 if i == j else 0 for j in range(len(A))] for i in range(len(A))]



l,a,d,m = [int(i) for i in input().split()]
r = [(10**(c+1)-1-a)//d for c in range(18)]
MOD=m

ans = 0
flag = 0

for c, rc in enumerate(r):
    if rc < 0: continue
    C=[[pow(10,c+1,MOD),1,0], [0,1,d], [0,0,1]]
    if c==0:
        C = matpow(C,rc+1,MOD)
        ansvec = [[ans],[a],[1]]
    else:
        lc = max(r[c-1]+1,0)

        C = matpow(C,min(l-1,rc)-lc+1,MOD)
        ansvec = [[ans],[a+d*(lc)],[1]]
    ans = matmul(C,ansvec,MOD)[0][0]

    if rc > l:
        break

print(ans)

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

ten = 0

def pow_(x, t):
  if t == 0:
    return 1
  if t % 2 == 1:
    return (pow_(x, t - 1) * x) % M
  else:
    y = pow_(x, t // 2)
    return (y * y) % M

def f(l):
  if l == 0:
    return 0
  if l % 2 == 1:
    pl = l - 1
    x = f(pl)
    return x * ten + 1
  else:
    pl = l // 2
    x = f(pl)
    return x * pow_(ten, pl) + x

def g(l):
  if l == 0:
    return 0
  if l % 2 == 1:
    pl = l - 1
    x = g(pl)
    return x * ten + B * pl
  else:
    pl = l // 2
    x = g(pl)
    return x * pow_(ten, pl) + x + B * pl * f(pl)

last = A + B * (L - 1)
ans = 0
ten = 10
for _ in range(18):
  l = ten // 10
  r = ten - 1
  
  if last < l:
    ten *= 10
    continue
  if A > r:
    ten *= 10
    continue
  na = 0 
  nl = 0 

  if A >= l:
    na = A
  else:
    na = (l - A + B - 1) // B * B + A
    na = min(na, last)
  
  nlast = 0
  if last <= r:
    nlast = last
  else:
    nlast = (r - A) // B * B + A
  nl = (nlast - na) // B + 1

  ans *= pow_(ten, nl) % M
  ans += (na * f(nl)) % M
  ans += g(nl) % M
  ten *= 10

print(ans % M)

Solution (Python):
def dot(a,b,m):
    c = [[0,0,0],[0,0,0],[0,0,0]]
    for i in range(3):
        for j in range(3):
            for k in range(3):
                c[i][j] += a[i][k]*b[k][j]
                c[i][j] %= m
    return c
    
def powm(a,n,m):
    if n == 0:
        return [[1,0,0],[0,1,0],[0,0,1]]
    for i in range(3):
        for j in range(3):
            a[i][j] %= m
    r = powm(a,n//2,m)
    if n%2==0:
        return dot(r,r,m)
    else:
        return dot(dot(r,r,m),a,m)

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

a = [[1,0,0],[0,1,0],[0,0,1]]
for d in range(1,19):
    C = min(max(0,(10**d-A-1)//B+1),L)-min(max(0,(10**(d-1)-A-1)//B+1),L)
    a = dot(powm([[pow(10,d,M),1,0],[0,1,B%M],[0,0,1]],C,M),a,M)
ans = (a[0][1]*A%M+a[0][2])%M
print(ans)

Solution (Python):
def unitmat(n):
    return [[int(i==j) for i in range(n)] for j in range(n)]
def matadd(A,B):
    h=len(A); w=len(A[0])
    C=[[A[i][j]+B[i][j] for j in range(w)] for i in range(h)]
    return C
def matmul(A,B):
    h=len(A)
    w=len(B[0])
    l=len(B)
    C=[[0]*w for _ in range(h)]
    for i in range(h):
        for j in range(w):
            for k in range(l):
                C[i][j]+=A[i][k]*B[k][j]
    for i in range(h):
        for j in range(w):
            C[i][j]%=mod
    return C
def matpow(A,k):
    n=len(A)
    B=unitmat(n)
    while k:
        if k&1:
            B=matmul(B,A)
        A=matmul(A,A)
        k>>=1
    return B

l,a,b,mod=map(int,input().split())
X=[[0],[a],[1]]
maxdigit=18
Idx=[0]*(maxdigit+1)
for i in range(maxdigit+1):
    if 10**i-1<a:
        continue
    Idx[i]=min((10**i-1-a)//b+1,l)
for i in range(1,maxdigit+1):
    f=[[10**i,1,0],[0,1,b],[0,0,1]]
    X=matmul(matpow(f,Idx[i]-Idx[i-1]),X)
print(X[0][0])

Solution (Python):
def unitmat(n):
    return [[int(i==j) for i in range(n)] for j in range(n)]
def matadd(A,B):
    h=len(A); w=len(A[0])
    C=[[A[i][j]+B[i][j] for j in range(w)] for i in range(h)]
    return C
def matmul(A,B):
    h=len(A)
    w=len(B[0])
    l=len(B)
    C=[[0]*w for _ in range(h)]
    for i in range(h):
        for j in range(w):
            for k in range(l):
                C[i][j]+=A[i][k]*B[k][j]
    for i in range(h):
        for j in range(w):
            C[i][j]%=mod
    return C
def matpow(A,k):
    n=len(A)
    B=unitmat(n)
    while k:
        if k&1:
            B=matmul(B,A)
        A=matmul(A,A)
        k>>=1
    return B

l,a,b,mod=map(int,input().split())
X=[[0],[a],[1]]
maxdigit=18
Idx=[0]*(maxdigit+1)
for i in range(maxdigit+1):
    if 10**i-1<a:
        continue
    Idx[i]=min((10**i-1-a)//b+1,l)
for i in range(1,maxdigit+1):
    f=[[10**i,1,0],[0,1,b],[0,0,1]]
    X=matmul(matpow(f,Idx[i]-Idx[i-1]),X)
print(X[0][0])

Solution (Python):
def mat_dot(one, two, mod):
	return [[sum([(one[i][k]*two[k][j])%mod for k in range(len(two))])%mod for j in range(len(two[0]))] for i in range(len(one))]

def mat_pow(mat, exp, mod):
	size = len(mat)
	res = [[0 for _ in range(size)] for _ in range(size)]
	for i in range(size):
		res[i][i] = 1
	cnt = 0
	while (1<<cnt) <= exp:
		if (exp>>cnt)&1:
			res = mat_dot(res, mat, mod)
		mat = mat_dot(mat, mat, mod)
		cnt += 1
	return res

l, a, b, m = map(int, input().split())
c = a + (l-1)*b
n = len(str(c))
dgt = [0 for _ in range(n+1)]
for i in range(1, n):
	dgt[i] = max(0, (10**i -a+b-1) // b) - max(0, (10**(i-1) -a+b-1) // b)
dgt[n] = l-sum(dgt)
d_sum = [0 for _ in range(n+1)]
for i in range(n-1, 0, -1):
	d_sum[i] = d_sum[i+1] + dgt[i+1]*(i+1)
fr = [0 for _ in range(n+1)]
fr[1] = a
for i in range(2, n+1):
	fr[i] = fr[i-1] + b*dgt[i-1]
l = [[0, a%m, 1]]
for d in range(1, n+1):
	k = [[(10**d)%m, 0, 0], [1, 1, 0], [0, b%m, 1]]
	j = mat_pow(k, dgt[d], m)
	l = mat_dot(l, j, m)

print(l[0][0]%m)

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

T = []
for w in range(0, len(str(A+B*(L-1)))):
  l0 = ((10**w-1)-A)/B+1
  l1 = ((10**(w+1)-1)-A)/B+1
  T += [(max(0, min(l1, L)-max(l0, 0)), A+B*max(0, l0), w+1)]



def f(L, w):
  if L==1:
    return 1
  t = f(L/2, 2*w)*(pow(10, w, M) + 1)%M
  if L%2!=0:
    t = (t*pow(10, w, M) + 1)%M
  return t

def g(L, B, w):
  if L==1:
    return 0
  t = (g(L/2, B*2, w*2)*(pow(10, w, M) + 1) + B*f(L/2, w*2))%M
  if L%2!=0:
    t = (t*pow(10, w, M) + B*(L-1))%M
  return t

ans = 0
w = 0
for t in T[::-1]:
  if t[0]>0:
    ans += (t[1]*f(t[0], t[2]) + g(t[0], B, t[2]))*pow(10, w, M)
    ans %= M
  w += t[2]*t[0]
print ans

Solution (Python):


MOD = None

def matunit(n):
    return [ [ int(y == x) for x in range(n) ] for y in range(n) ]

def matadd(a, b):
    assert len(a) == len(b)
    assert len(a[0]) == len(b[0])
    h = len(a)
    w = len(a[0])
    c = [ [ a[y][x] + b[y][x] for x in range(w) ] for y in range(h) ]
    return c

def matmul(a, b):
    assert len(a[0]) == len(b)
    h = len(a)
    k = len(b)
    w = len(b[0])
    c = [ [ 0 for _ in range(w) ] for _ in range(h) ]
    for y in range(h):
        for z in range(k):
            for x in range(w):
                c[y][x] += a[y][z] * b[z][x]
    for y in range(h):
        for x in range(w):
            c[y][x] %= MOD
    return c

def matpow(a, k):
    assert len(a) == len(a[0])
    n = len(a)
    b = matunit(n)
    while k:
        if k & 1:
            b = matmul(b, a)
        a = matmul(a, a)
        k >>= 1
    return b

def binsearch(l, r, pred): 
    assert l < r
    l -= 1
    while r - l > 1:
        m = (l + r) // 2
        if pred(m):
            r = m
        else:
            l = m
    return r

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

    global MOD
    MOD = m

    x = [ [ 0 ], [ a ], [ 1 ] ]  
    DIGITS = 20
    d = [ binsearch(0, l, lambda i: len(str(a + b * i)) >= d) for d in range(DIGITS) ]
    for i in range(DIGITS - 1):
        l = d[i]
        r = d[i + 1]
        f = [ [ 10 ** i, 1, 0 ], [ 0, 1, b ], [ 0, 0, 1 ] ]
        x = matmul(matpow(f, r - l), x)
    print(x[0][0] % m)


if __name__ == '__main__':
    main()

Solution (Python):
def S1(r,T,M):
  
  
  
  if r == 1:
    return T%M
  x = pow(r,T,M*(r-1)) - 1
  return (x//(r-1))%M

def S2(r,T,M):
  
  if r == 1:
    x = T*(T-1)//2
    return x%M
  Mr = M*(r-1)
  x = (T-1) * pow(r,T,Mr)
  x -= S1(r,T,Mr)
  x += 1
  x //= (r-1)
  return x

def S(A,B,L,M):
  result = 0
  d = len(str(A)) 
  start = A
  r = 10**d
  while True:
    n1 = ((r-1)-A)//B+1
    items = min(n1,L) 
    last = A + (items-1)*B
    x = last * S1(r,items,M)
    x -= B * S2(r,items,M)
    result *= pow(r,items,M)
    result += x
    result %= M
    
    r *= 10
    d += 1 
    L -= items
    A += B*items
    if L <= 0:
      break
  return result%M

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

Solution (Python):
def solve(L, A, B, M):
    ret, digits = 0, 0
    for d in range(18, 0, -1):
        lo = 10 ** (d - 1)
        hi = 10 ** d - 1
        if hi < A or A + B * (L - 1) < lo:
            continue
        imin = 0 if lo <= A else (lo - A + B - 1) // B
        init = A + B * imin
        if init > hi:
            continue
        imax = L - 1 if A + B * (L - 1) <= hi else imin + (hi - init) // B
        n = imax - imin + 1
        p = 10 ** d % M
        a = matpow([[1, 0, 1, 0], [0, 1, 0, 1], [0, 0, p, p], [0, 0, 0, p]], n, M)
        
        
        
        sub = (A + B * imax) % M * a[1][3] % M + M - (B * a[0][3] % M)
        ret += sub % M * pow10(digits, M) % M
        digits += d * (imax - imin + 1)
    return ret % M


def pow10(p, mod):
    if p % 2 == 1:
        return 10 * pow10(p - 1, mod) % mod
    elif p > 0:
        sub = pow10(p // 2, mod)
        return sub * sub % mod
    else:
        return 1


def matpow(a, p, mod):
    if p % 2 == 1:
        return matmul(a, matpow(a, p - 1, mod), mod)
    elif p > 0:
        b = matpow(a, p // 2, mod)
        return matmul(b, b, mod)
    else:
        n = len(a)
        return [[1 if i == j else 0 for j in range(n)] for i in range(n)]


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


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