Problem:


Score : 600 points

Problem Statement
Q tuples of integers (L_i, A_i, B_i, M_i). For each tuple, answer the following question. -->

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



Sample input:
5 3 4 10007


Sample output:
5563

Submission status: Accepted

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

Submission status: Accepted

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

# c111 -> 1001001001みたいなやつを求める
# n : 塊の数
# l : 一つの塊(00..001)の長さ
# m : mod
# ex) n=3, l=2 -> 10101
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

# c123 -> 1002003004みたいなやつを求める
# n : 塊の数
# l : 一つの塊(00..001)の長さ
# m : mod
# ex) n=3, l=2 -> 10203
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)

Submission status: Accepted

Solution code:
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
# 全部つなげたのを X とすると、
# [今の X, 次に加算する数, 公差]
vec = [ModInt(0), ModInt(A), ModInt(B)]
# 加算する数の桁数
digits = 1
while size < L:
    # 桁数が digits である要素の数
    # digits 以下の要素の数の合計から今までに加算した数を引く
    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])

Submission status: Accepted

Solution code:
import copy
import sys

stdin = sys.stdin

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

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


# x = x * high + val
# val += B
# (high 1 0)
# (0 1 1)
# (0 0 1)

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)
    # print(low, high, step, ra + B*step, v)
    ra = ra + B * step
    L -= step

    low *= 10
    high *= 10

print(v[0])

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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):
    #k-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)

Submission status: Accepted

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

    # https://github.com/numpy/numpy/blob/master/numpy/linalg/linalg.py L560から
    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

    # 数列に対して、d桁の要素の個数を求める
    num_pow_d = [0] * (DIGIT_MAX + 1)
    for d in range(DIGIT_MAX + 1):
        if 10**d - 1 >= A:
            # dに対してA + B*n <= 10**d - 1 となる最大のn(L以下しか取りえない)
            num_pow_d[d] = min(((10**d - 1 - A) // B) + 1, L)  # +1はnが0スタートのため

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

Submission status: Accepted

Solution code:
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])

Submission status: Accepted

Solution code:
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])

Submission status: Accepted

Solution code:
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
 
    # 数列に対して、d桁の要素の個数を求める
    num_pow_d = [0] * (DIGIT_MAX + 1)
    for d in range(DIGIT_MAX + 1):
        if 10**d - 1 < A:
            continue
        # dに対してA + B*n <= 10**d - 1 となる最大のn(L以下しか取りえない)
        num_pow_d[d] = min(((10**d - 1 - A) // B) + 1, L)  # +1はnが0スタートのため
    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))

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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])

Submission status: Accepted

Solution code:
# coding: utf-8
# Your code here!

def matmul(A, B, mod): # A,B: 行列
    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): #A^p
    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
#print(r)
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)
#        print(lc)
        C = matpow(C,min(l-1,rc)-lc+1,MOD)
        ansvec = [[ans],[a+d*(lc)],[1]]
    ans = matmul(C,ansvec,MOD)[0][0]
#    print(ans)
    if rc > l:
        break

print(ans)

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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])

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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)]
#for t in T:
#  print t

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

Submission status: Accepted

Solution code:
#!/usr/bin/env python3

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): # [l, r)
    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 ] ]  # (answer, s_i, 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()

Submission status: Accepted

Solution code:
def S1(r,T,M):
  # 1+r+...+r^{T-1}
  # (r^T-1)/(r-1)
  # mod M(r-1)で計算する
  if r == 1:
    return T%M
  x = pow(r,T,M*(r-1)) - 1
  return (x//(r-1))%M

def S2(r,T,M):
  # sum ir^i
  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)

Submission status: Accepted

Solution code:
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)
        # a[0][3] = sum p^i for i in [0, n-1]
        # a[1][3] = sum i * p^i for i in [0, n-1]
        # sum (A + B * (imax - i)) * p^i for i in [0, n-1]
        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))