import sys
MOD = 10**9 + 7

def main():
    import sys
    sys.setrecursionlimit(1 << 25)
    p, n, k, b = map(int, sys.stdin.readline().split())
    a = list(map(int, sys.stdin.readline().split()))
    
    # 预处理 pow_k[x] = x^k mod p
    pow_k = [0] * p
    for x in range(p):
        pow_k[x] = pow(x, k, p)
    
    # 处理 a_i=0 的情况
    m = 0
    for ai in a:
        if ai == 0:
            m += 1
    # 这些变量的贡献总是0，所以总共有 p^m 种情况
    ans_p_m = pow(p, m, MOD)
    
    # 剩下的变量
    a_nonzero = [ai for ai in a if ai != 0]
    n_nonzero = len(a_nonzero)
    if n_nonzero == 0:
        if b == 0:
            print(ans_p_m % MOD)
        else:
            print(0)
        return
    
    # 现在处理剩下的n_nonzero个变量
    # 预处理每个变量i的贡献频率
    freqs = []
    for ai in a_nonzero:
        freq = [0] * p
        for x in range(p):
            v = (ai * pow_k[x]) % p
            freq[v] += 1
        freqs.append(freq)
    
    # 现在，我们需要计算这些freq的多项式乘积
    # 初始化dp
    dp = [0] * p
    dp[0] = 1
    
    # 迭代计算每个变量的贡献
    for freq in freqs:
        new_dp = [0] * p
        for j in range(p):
            if dp[j] == 0:
                continue
            for v in range(p):
                if freq[v] == 0:
                    continue
                new_v = (j + v) % p
                new_dp[new_v] = (new_dp[new_v] + dp[j] * freq[v]) % MOD
        dp = new_dp
    
    # 最终答案是 dp[b] * ans_p_m mod MOD
    total = (dp[b] * ans_p_m) % MOD
    print(total)

if __name__ == "__main__":
    main()