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())
MOD=M

#n*n行列の積を繰り返し二乗法で求める
#単位行列
def imat(n):
  ret=[[1 if i==j else 0 for j in range(n)] for i in range(n)]
  return ret

#行列の積A*B
def prod_mat(amat,bmat):
  res_mat = [[sum([amat[i][j]*bmat[j][k] for j in range(len(bmat))]) for k in range(len(bmat[0]))] for i in range(len(amat))]
  return res_mat  

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

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

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

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

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

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

Submission status: Accepted

Solution code:
import sys
import time

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


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

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

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

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

    print(res[0])


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
import numpy as np

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

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

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

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

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

Submission status: Accepted

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

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

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

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

    return ansline*a+ansconst*b

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

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

Submission status: Accepted

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

ten = 1
def f(x):

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

def g(x):

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

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

print(ans)

Submission status: Accepted

Solution code:

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

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

"""

#行列A,Bの積(不可能だとエラー)
def matrix_mul(A,B):

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

    for ai in range(len(A)):

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

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

    return ans

#行列Aのx乗(当然正方行列じゃないとだめ)
def matrix_pow(A,x):

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

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

        x//=2
    return ans

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

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

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

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

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

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

    #print (i,k,now)
print (now[0][0] % mod)

Submission status: Accepted

Solution code:
import sys
import itertools
import queue
input = sys.stdin.readline
sys.setrecursionlimit(100000)
# mod = 10 ** 9 + 7

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


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


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


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

    return res


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

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

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

    print(X)


if __name__ == "__main__":
    main()

Submission status: Accepted

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

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

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

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

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

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


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
import numpy as np

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

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

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

Submission status: Accepted

Solution code:
from numpy.linalg import matrix_power
from numpy import dot, array

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

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

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

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


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

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

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

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

    size += cnt
    if size >= L:
        break

print(ans[0])

Submission status: Accepted

Solution code:
from numpy.linalg import matrix_power
from numpy import dot

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

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

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

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


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

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

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

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

    size += cnt
    if size >= L:
        break

print(ans[0])

Submission status: Accepted

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

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

# mat = [[1, 0, 0],                                                                                                                                                                                                                                            
#        [0, 1, 0],                                                                                                                                                                                                                                            
#        [0, 0, 1]]                                                                                                                                                                                                                                            
# vec = (0, a, 1)                                                                                                                                                                                                                                              

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

# print((mat[0][0] * vec[0] + mat[0][1] * vec[1] + mat[0][2] * vec[2]) % mod)                                                                                                                                                                                  
print(mat[0][0])

Submission status: Accepted

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

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

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

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


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

SL = A + (L-1) * B
ts = 0
ks = [0]
for i in range(len(str(SL))+1):
    # ks[i] = Count(10**i <= Si <10**(i+1))
    ks.append(max(min((10**(i+1)-1),SL)-(A-B), 0) // B - ts)
    ts += ks[-1]
ti = 0
a = SL
R = 0
# i = 2, n = 3, a = 17, B=2
# 17 * 1500 * 130000
# r += Sum_{0<=j<=n-1}((a-Bj)*pow(10**i, j)) * pow(10, ti)
# j = 0 -> (a+Bj)*pow(10**i, j) = a
for i in range(len(str(SL))+1,0,-1):
    n = ks[i]
    R = (R + f(a, B, n, i)*power(10, ti)) % mod
    ti += i*n
    a -= n*B
print(R)

Submission status: Accepted

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

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

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

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


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

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

    ti += i*n
    ti2 += n
print(R)
#$print(ks)

Submission status: Accepted

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


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



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


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

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

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

    return y


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


print(ret[0][0])

Submission status: Accepted

Solution code:
from operator import mul

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

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

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

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

print(As[0][0])

Submission status: Accepted

Solution code:
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

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

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

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

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

print(S)

Submission status: Accepted

Solution code:
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

from functools import lru_cache

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

@lru_cache(None)
def F(r,n,MOD):
    # 1+r+...+r^{n-1} をダブリングで求める
    if n == 0:
        return 0
    q = n//2
    x = F(r,q,MOD)
    # 1+r^qをかける
    x *= 1 + pow(r,q,MOD); x %= MOD
    if n&1:
        x *= r; x += 1; x %= MOD
    return x

@lru_cache(None)
def G(r,n,MOD):
    # sum ir^i をダブリングで求める
    if n == 0:
        return 0
    q = n//2
    x = G(r,q,MOD)
    # 1+r^qをかける
    x *= 1 + pow(r,q,MOD); x %= MOD
    # q(r^q+...+r^{2q-1})を加える
    x += F(r,q,MOD) * pow(r,q,MOD) * q % MOD
    if n&1:
        x *= r; x += F(r,n,MOD); x -= 1; x %= MOD
    return x

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

print(S)

Submission status: Accepted

Solution code:
#!usr/bin/env python3
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

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

#B
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

#C
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

#D
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

#E
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

#F
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)) #left := a+bi >= pを満たす最小の整数iだがceilを使うとなぜかバグった、そのためfloorで-leftを求めてからleftを求めた
        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

#Solve
if __name__ == "__main__":
    F()