Everyone knows what the Fibonacci sequence is. This sequence can be defined by the recurrence relation: 

We'll define a new number sequence *A*<sub class="lower-index">*i*</sub>(*k*) by the formula: 

In this problem, your task is to calculate the following sum: *A*<sub class="lower-index">1</sub>(*k*) + *A*<sub class="lower-index">2</sub>(*k*) + ... + *A*<sub class="lower-index">*n*</sub>(*k*). The answer can be very large, so print it modulo 1000000007 (10<sup class="upper-index">9</sup> + 7).
The first line contains two space-separated integers *n*, *k* (1 ≤ *n* ≤ 10<sup class="upper-index">17</sup>; 1 ≤ *k* ≤ 40).
def add_mat(a, b, mod):
	return [[(a[i][j] + b[i][j]) % mod for j in range(len(a[0]))] for i in range(len(a))]


def mul_mat(a, b, mod):
	res = []
	for i in range(len(a)):
		res.append([])
		for j in range(len(b[0])):
			res[i].append(sum(a[i][_] * b[_][j] for _ in range(len(a[0]))) % mod)

	return res


def pow_mat(b, n, mod):
	if n == 1:
		return b

	tmp = pow_mat(b, n//2, mod)
	res = mul_mat(tmp, tmp, mod)
	if n % 2 == 1:
		res = mul_mat(res, b, mod)

	return res


def generate_coefs(n, a0, a1=1, mod=None):
	"""
	Generate coefficients in the expansion of (a0+a1*x)^n.
	"""
	binom = [1]
	for i in range(n):
		binom.append((n-i) * binom[-1] // (i+1))

	res = []
	for i in range(n+1):
		res.append(binom[i] * a1**i * a0**(n-i))
		if mod is not None:
			res[i] %= mod

	return res


def calculate(n, k, mod):
	"""
	Calculate the sum: (A_1(k) + A_2(k) + ... + A_n(k)) % mod
	with F_1 = 1, F_2 = 2, F_i = F_{i-1} + F_{i-2} (i > 2);
	     A_i(k) = F_i × i^k (i >= 1).
	"""
	if n <= 1:
		return n

	ker_size = 2*(k+1)+1
	A = [None] * ker_size
	for i in range(ker_size):
		A[i] = [0] * ker_size

	A[ker_size-1][ker_size-1] = 1

	for x in range(k+1):
		A[k+1+x][x] = 1
		coefs_1 = generate_coefs(x, 1, 1)
		coefs_2 = generate_coefs(x, 2, 1)
		# print(x, coefs_1, coefs_2)

		for i in range(x+1):
			A[i][k+1+x] = coefs_2[i]
			A[k+1+i][k+1+x] = coefs_1[i]

			if x == k:
				A[i][k+2+x] = coefs_2[i]
				A[k+1+i][k+2+x] = coefs_1[i]

	# print(*A, sep='\n')
	b0 = [[1, 0, 0, 0, 1, 1, 1, 1, 1]]
	# fn1 = f(n-1)
	# b(n-1)  = [fn1*(n-1)**0..k, fn*n**0..k, s]
	# => b(n) = [fn*n**0..k, (fn1+fn)*(n+1)**0..k, s+...]
	bn = mul_mat(b0, pow_mat(A, n, mod), mod)
	# print(bn)

	return bn[0][-1]


if __name__ == '__main__':
	n, k = map(int, input().split())
	print(calculate(n, k, 10**9 + 7))
