I found an interesting problem on https://codeforces.com/problemset/problem/1105/C:
Ayoub had an array $a$ of integers of size $n$ and this array had two interesting properties:   All the integers in the array were between $l$ and $r$ (inclusive).  The sum of all the elements was divisible by $3$. 

Unfortunately, Ayoub has lost his array, but he remembers the size of the array $n$ and the numbers $l$ and $r$, so he asked you to find the number of ways to restore the array. 

Since the answer could be very large, print it modulo $10^9 + 7$ (i.e. the remainder when dividing by $10^9 + 7$). In case there are no satisfying arrays (Ayoub has a wrong memory), print $0$.


-----Input-----

The first and only line contains three integers $n$, $l$ and $r$ ($1 \le n \le 2 \cdot 10^5 , 1 \le l \le r \le 10^9$) — the size of the lost array and the range of numbers in the array.


-----Output-----

Print the remainder when dividing by $10^9 + 7$ the number of ways to restore the array.


-----Examples-----
Input
2 1 3

Output
3

Input
3 2 2

Output
1

Input
9 9 99

Output
711426616



-----Note-----

In the first example, the possible arrays are : $[1,2], [2,1], [3, 3]$.

In the second example, the only possible array is $[2, 2, 2]$.

I tried it in Python, but could not do it. Can you solve it?
from sys import stdin, stdout
prime = 10**9 + 7

def unitMatrix(size):
    out = [[0]*size for i in range(size)]
    for i in range(size):
        out[i][i] = 1
    return out

def matrixMult(pre, post):
    rows = len(pre)
    mid = len(post)
    columns = len(post[0])
    out = [[0]*columns for i in range(rows)]
    for i in range(rows):
        for j in range(columns):
            for k in range(mid):
                out[i][j] += pre[i][k]*post[k][j]
                out[i][j] = out[i][j]%prime
    return out

def matrixExp(base, power):
    if power == 0:
        return unitMatrix(len(base))
    else:
        if power%2 == 0:
            half = matrixExp(base, power/2)
            return matrixMult(half, half)
        else:
            half = matrixExp(base, (power-1)/2)
            return matrixMult(matrixMult(half, half), base)

def main():
    n, l, r = [int(i) for i in stdin.readline().split()]

    vals = [[1], [0], [0]]
    mods = [(r-l+1)//3 ,(r-l+1)//3 ,(r-l+1)//3]
    remainder = (r-l+1)%3
    if remainder >= 1:
        mods[l%3] += 1
    if remainder == 2:
        mods[(l+1)%3] += 1
        
    transforms = [[mods[0], mods[2], mods[1]], [mods[1], mods[0], mods[2]], [mods[2], mods[1], mods[0]]]
    power = matrixExp(transforms, n)
    out = matrixMult(power, vals)
    print(out[0][0])

main()