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 if lca_depth != depth1 else 1) * 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)