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:
#!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))
        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()

Submission status: Accepted

Solution code:
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)#n//2を二回行う
    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])

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.readline
#from numpy import matrix, eye, dot
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): #N×N行列のmod
    for i in range(N):
        for j in range(N):
            A[i][j] %= Q
    return

def powOfMatrix(N, X, n, Q): #N×N行列のn乗
    Ret = [[1,0,0],[0,1,0],[0,0,1]]
    power = '{:060b}'.format(n)[::-1] #log2(pow(10,18)) < 60
    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()

Submission status: Accepted

Solution code:
#import sys
#input = sys.stdin.readline
#from numpy import matrix, eye, dot
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): #N×N行列のmod
    for i in range(N):
        for j in range(N):
            A[i][j] %= Q
    return

def powOfMatrix(N, X, n, Q): #N×N行列のn乗
    Ret = [[1,0,0],[0,1,0],[0,0,1]]
    power = '{:060b}'.format(n)[::-1] #log2(pow(10,18)) < 60
    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()

Submission status: Accepted

Solution code:
#import sys
#input = sys.stdin.readline
#from numpy import matrix, eye, dot
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): #N×N行列のmod
    for i in range(N):
        for j in range(N):
            A[i][j] %= Q
#            A[i,j] %= Q
    return

def powOfMatrix(N, X, n, Q): #N×N行列のn乗
#    Ret = eye(N)
    Ret = [[1,0,0],[0,1,0],[0,0,1]]
    power = '{:060b}'.format(n)[::-1] #log2(pow(10,18)) < 60
    for p in power:
        if p == "1":
            Ret = productMatrix(N,Ret, X)
#            Ret = dot(Ret, X)
            modMatrix(N, Ret, Q)
        X = productMatrix(N,X,X)
#        X = dot(X,X)
        modMatrix(N, X, Q)
#        print(Ret)
    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):
#        print(i, ANS)
        if s >= pow(10,i):
            continue
        P = [[pow(10, i, M), 0, 0], [1,1,0], [0,B,1]]
#        P = matrix([[pow(10, i, M), 0, 0], [1,1,0], [0,B,1]])
#        print((pow(10,i)-1-s)//B+1)
        step = (pow(10,i)-s+B-1)//B
        if L <= step:
            ANS = productMatrix(3,ANS,powOfMatrix(3,P,L,M))
#            ANS = dot(ANS,powOfMatrix(3,P,L,M))
            modMatrix(3,ANS,M)
            break
#        print(powOfMatrix(3,P, (pow(10,i)-1-s)//B + 1,M))
        ANS = productMatrix(3,ANS, powOfMatrix(3,P, step,M))
#        ANS = dot(ANS, powOfMatrix(3,P, (pow(10,i)-1-s)//B + 1,M))
#        print( (ANS[1][0]*A + ANS[2][0])%M, (pow(10,i)-1-s)//B + 1)
        modMatrix(3,ANS, M)
        L -= step
        s += step*B
    print( (ANS[1][0]*A + ANS[2][0])%M)
#    print( int((ANS[1,0]*A + ANS[2,0])%M))
        
        
if __name__ == '__main__':
    main()

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

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

def keta(k):# 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): # 1+r+...+r**(n-1)
    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): # 0+r+2*r**2+...+(n-1)*r**(n-1)
    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
    

    #print(i,sh,ma,kou)

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

print(ANS)

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(pow(10, 6))
# input = sys.stdin.readline

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

Submission status: Accepted

Solution code:
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)
    #print(temp_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, s)

print(x)

Submission status: Accepted

Solution code:
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)
  #print(d, R, s)
  if (C[d] > L):
    break

print(s[0])

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(100000000)

def mult(A, B, mod): #行列A, Bの内積を求める。
    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()

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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 
    #ans += powB(1, L - T[k], D, M) * t % M
    #ans += powA(-T[k], 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)

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

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

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

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

Submission status: Accepted

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

# first i such that A + i*B has d digits.
# A + i*B >= 10^d
# i*B >= 10^d - A
# i >= 
# i = (10^d - A + (B-1))/B

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)
	
	

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

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