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: Wrong Answer

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 = mul(n%2, (Y+mul(n, power(10, n*i))))
    return Y


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))
    #bb = 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: Wrong Answer

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 = mul(n%2, (Y+mul(n, power(10, n*i))))
    return Y


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))
    bb = mul(b, div(r-n*P + (n-1)*P*r, mul(1-r, 1-r)))
    bb2 = 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: Wrong Answer

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):
    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)
    return (
            (
                (X+tn*(G(tn+1, i)-1) % mod) * power(10, tn*i) % mod
                +X
                )
            + (n%2)*(n*i)*power(10, n*i) % mod
            ) % mod

#NNN = (10**5)
#inverse = [0, 1]
#for i in range( 2, NNN + 1 ):
#    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )


def f(a, b, r, n, i):
    if n == 0:
        return 0
    P = power(r,n)
    aa = mul(a, G(n, i))
    bb = mul(b, G2(n-1, i))
    #print("aa")
    #print(a, aa, bb)
    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
    r = power(10, i)
    n = ks[i-1]
    a = sl - ti2*b
    if n:
        R = (R + (f(a, b, r, 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: Wrong Answer

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 G(n, i):
    if n == 1:
        return 1
    tn = n//2
    return (G(tn, i) * (1 + power(10, tn*i)) + (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)
    return (
            (
                (X+tn*(G(tn+1, i)-1)) * power(10, tn*i)
                +X
                )
            + (n%2)*(n*i)*power(10, n*i)
            ) % mod

#NNN = (10**5)
#inverse = [0, 1]
#for i in range( 2, NNN + 1 ):
#    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )

def div(a, b):
    #return mul(a, power(b, mod-2))
    #return (a//b) % mod
    return a * modinv(b) % mod

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

sl = A + (L-1) * B
#print(sl)
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
    r = power(10, i)
    n = ks[i-1]
    a = sl - ti2*b
    if n:
        R = (R + (f(a, b, r, 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: Wrong Answer

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

#NNN = (10**5)
#inverse = [0, 1]
#for i in range( 2, NNN + 1 ):
#    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )

def div(a, b):
    #return mul(a, power(b, mod-2))
    #return (a//b) % mod
    return a * modinv(b) % mod

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

sl = A + (L-1) * B
#print(sl)
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
    r = power(10, i)
    n = ks[i-1]
    a = sl - ti2*b
    if n:
        R = (R + (f(a, b, r, n))*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: Wrong Answer

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
    # 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: Wrong Answer

Solution code:
import sys
from io import StringIO
import unittest
import math
class TestClass(unittest.TestCase):
    def assertIO(self, input, output):
        stdout, stdin = sys.stdout, sys.stdin
        sys.stdout, sys.stdin = StringIO(), StringIO(input)
        resolve()
        sys.stdout.seek(0)
        out = sys.stdout.read()[:-1]
        sys.stdout, sys.stdin = stdout, stdin
        self.assertEqual(out, output)
    def test_入力例_1(self):
        input = """5 3 4 10007"""
        output = """5563"""
        self.assertIO(input, output)
    def test_入力例_2(self):
        input = """4 8 1 1000000"""
        output = """891011"""
        self.assertIO(input, output)
    def test_入力例_3(self):
        input = """107 10000000000007 1000000000000007 998244353"""
        output = """39122908"""
        self.assertIO(input, output)
def resolve():
    (L,A,B,M)=map(int,input().split())
    st=(L-1)*B+A
    beforemod=1
    modlist=[0,1]
    mod=0
    for i in range(19):
        modlist.append(modlist[-1]*10%M)
    for i in range(L):
        mod=((beforemod*(st%M))%M+mod)%M

        beforemod=(modlist[math.floor(math.log10(st)+1)+1]*beforemod)%M        

        st-=B

    print(mod)





if __name__ == "__main__":
    unittest.main()

Submission status: Wrong Answer

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,20):
        left = max(0,math.ceil((p-a)/b))
        if left >= l:
            break
        p *= 10
        right = min(l,math.ceil((p-a)/b))
        if right <= 0:
            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: Wrong Answer

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,20):
        left = max(0,math.ceil((p-a)/b))
        if left >= l:
            break
        p *= 10
        right = min(l,math.ceil((p-a)/b))
        if right <= left:
            continue
        mat = mat_pow([[p,1,0],[0,1,b],[0,0,1]],right-left,m)
        ans *= pow(10,d*(right-left),m)
        ans += (mat[0][1]*(a+b*left)%m+mat[0][2]%m)
        ans %= m
    print(ans)
    return

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

Submission status: Wrong Answer

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,20):
        left = max(0,math.ceil((p-a)/b))
        if left >= l:
            break
        p *= 10
        right = min(l,math.ceil((p-a)/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)
        ans %= m
    print(ans)
    return

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

Submission status: Wrong Answer

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,math.ceil((p-a)/b))
        if left >= l:
            break
        p *= 10
        right = min(l,math.ceil((p-a)/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)
        ans %= m
    print(ans)
    return

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

Submission status: Wrong Answer

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)
        if L < (pow(10,i)-1-s)//B+1:
            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, (pow(10,i)-1-s)//B + 1,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 -= (pow(10,i)-1-s)//B + 1
        s = ((pow(10,i)-1-s)//B+1)*B
    print( (ANS[1][0]*A + ANS[2][0])%M)
#    print( int((ANS[1,0]*A + ANS[2,0])%M))
        
        
if __name__ == '__main__':
    main()