d, l, r, k = map(int, input().split()) mod = 10**9 + 7 ans = 1 fac = [1]*(2**(d-1)+1) for i in range(1, 2**(d-1)+1): fac[i] = (fac[i-1] * i) % mod depth1, depth2 = -1, -1 for depth in range(d): i_l, i_r = 2**depth, 2**(depth+1)-1 nodes = 2**depth if i_l <= l <= i_r: depth1 = depth nodes -= 1 if i_l <= r <= i_r: depth2 = depth nodes -= 1 ans = (ans * fac[nodes]) % mod if depth1 != depth2: lca_depth = depth1 - (k - (depth2 - depth1)) // 2 need_edge = (depth1 - lca_depth) * 2 + (depth2 - depth1) # print(depth1, depth2, lca_depth, need_edge) if (k - (depth2 - depth1)) % 2 == 1 or lca_depth < 0 or depth1 < lca_depth or k != need_edge: ans = 0 else: ans = ans * pow(2, depth2 - lca_depth, mod) // 2 * pow(2, depth1, mod) % mod else: if k % 2 == 1 or (d-1)*2 < k: ans = 0 else: nodes = 2 ** depth1 unit = 2 ** (k//2) n = nodes // unit ans = ans * (pow(unit//2, 2, mod) * n) * 2 % mod print(ans % mod)