Xor-sequences
You are given n integers a1,  a2,  ...,  an.A sequence of integers x1,  x2,  ...,  xk is called a "xor-sequence" if for every 1  ≤  i  ≤  k - 1 the number of ones in the binary representation of the number xi  xi  +  1's is a multiple of 3 and  for all 1 ≤ i ≤ k. The symbol  is used for the binary exclusive or operation.How many "xor-sequences" of length k exist? Output the answer modulo 109 + 7.Note if a = [1, 1] and k = 1 then the answer is 2, because you should consider the ones from a as different.
The first line contains two integers n and k (1 ≤ n ≤ 100, 1 ≤ k ≤ 1018) — the number of given integers and the length of the "xor-sequences".The second line contains n integers ai (0 ≤ ai ≤ 1018).
Print the only integer c — the number of "xor-sequences" of length k modulo 109 + 7.
mod = 10 ** 9 + 7

def popcount(n):
    c = (n & 6148914691236517205) + (n >> 1 & 6148914691236517205)
    c = (c & 3689348814741910323) + (c >> 2 & 3689348814741910323)
    c = (c & 1085102592571150095) + (c >> 4 & 1085102592571150095)
    c = (c & 71777214294589695) + (c >> 8 & 71777214294589695)
    c = (c & 281470681808895) + (c >> 16 & 281470681808895)
    c = (c & 4294967295) + (c >> 32 & 4294967295)
    return c

def matmul(A, B):
    res = [[0] * len(B[0]) for _ in [None] * len(A)]
    for (i, resi) in enumerate(res):
        for (k, aik) in enumerate(A[i]):
            for (j, bkj) in enumerate(B[k]):
                resi[j] += aik * bkj
                resi[j] %= mod
    return res

def matpow(A, p):
    if p % 2:
        return matmul(A, matpow(A, p - 1))
    elif p > 0:
        b = matpow(A, p // 2)
        return matmul(b, b)
    else:
        return [[int(i == j) for j in range(len(A))] for i in range(len(A))]
(n, k) = map(int, input().split())
a = list(map(int, input().split()))
mat = [[0] * n for i in range(n)]
for i in range(n):
    for j in range(n):
        if popcount(a[i] ^ a[j]) % 3 == 0:
            mat[i][j] = 1
matp = matpow(mat, k - 1)
ansmat = matmul(matp, [[1] for i in range(n)])
ans = 0
for i in range(n):
    ans += ansmat[i][0]
print(ans % mod)
import sys
input = sys.stdin.readline

def matmul(a, b, n):
    c = [[0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            x = 0
            for k in range(n):
                x += a[i][k] * b[k][j] % mod
                x %= mod
            c[i][j] = x
    return c
(n, k) = map(int, input().split())
mod = pow(10, 9) + 7
a = list(map(int, input().split()))
b = [set() for _ in range(n)]
for i in range(n):
    a0 = bin(a[i])
    for j in range(1, len(a0)):
        if a0[-j] == '1':
            b[i].add(j)
x = [[0] * n for _ in range(n)]
for i in range(n):
    x[i][i] = 1
y = [-1] * (n * n)
y = [[-1] * n for _ in range(n)]
for i in range(n):
    y[i][i] = 1
for i in range(n):
    bi = b[i]
    for j in range(i + 1, n):
        bj = b[j]
        c = 0
        for l in bj:
            if not l in bi:
                c += 1
        for l in bi:
            if not l in bj:
                c += 1
        s = 0 if c % 3 else 1
        y[i][j] = s
        y[j][i] = s
k -= 1
k = list(bin(k)[2:])
while k:
    if k.pop() == '1':
        x = matmul(x, y, n)
    y = matmul(y, y, n)
ans = 0
for i in x:
    ans += sum(i) % mod
    ans %= mod
sys.stdout.write(str(ans))
import sys, math
input = sys.stdin.readline
INF = int(1000000000.0) + 7
mod = INF

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

def many(a, m):
    if m == 1:
        return a
    if m % 2 == 0:
        k = many(a, m // 2)
        return matrixmult(k, k)
    else:
        k = many(a, m // 2)
        return matrixmult(matrixmult(k, k), a)

def solve():
    (n, k) = map(int, input().split())
    data = list(map(int, input().split()))
    a = [1] * n
    base = [[0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if bin(data[i] ^ data[j])[2:].count('1') % 3 == 0:
                base[i][j] = 1
            else:
                base[i][j] = 0
    if k == 1:
        print(n)
    else:
        result = many(base, k - 1)
        ans = 0
        for i in range(n):
            ans += sum(result[i])
        print(ans % mod)
t = 1
while t:
    t -= 1
    solve()