def cmb(n, r, mod):#コンビネーションの高速計算　
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

mod = 998244353
N = 2*10**5
g1 = [1]*(N+1) # 元テーブル
g2 = [1]*(N+1) #逆元テーブル
inverse = [1]*(N+1) #逆元テーブル計算用テーブル

for i in range( 2, N + 1 ):
    g1[i]=( ( g1[i-1] * i ) % mod )
    inverse[i]=( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2[i]=( (g2[i-1] * inverse[i]) % mod )
inverse[0]=0

mod = 998244353

import numpy as np

N = 2**18
inv = [1]*(N+1) #逆元テーブル計算用テーブル

for i in range( 2, N + 1 ):
    inv[i]=( ( -inv[mod % i] * (mod//i) ) % mod )
inv[0]=0

inv = np.array(inv,np.int64)

def convolve(f,g,limit):
    fft_len =1
    while 2*fft_len<len(f)+len(g)-1:
        fft_len *= 2
    fft_len *=2
    Ff = np.fft.rfft(f,fft_len)
    Fg = np.fft.rfft(g,fft_len)

    Fh = Ff * Fg

    h=np.fft.irfft(Fh,fft_len)

    h= np.rint(h).astype(np.int64)

    return h[:min(len(f)+len(g)-1,limit)]

def convolve2(f,g,limit,p=998244353):

    f1,f2=np.divmod(f,1<<15)
    g1,g2=np.divmod(g,1<<15)

    a = convolve(f1,g1,limit)%p
    c = convolve(f2,g2,limit)%p
    b = (convolve(f1+f2,g1+g2,limit)-(a+c))%p

    h = (a<<30) + (b<<15) +c
    return h[:limit] % p

def inverse(f,limit):
    g = np.array([pow(int(f[0]),mod-2,mod)],np.int64)
    f = list(f)
    n = (len(f)-1).bit_length()
    F = f + [0]*(2**n-len(f))

    f=np.array(F,np.int64)

    for i in range(1,n+1):
        h = convolve2(g,f[:2**i],2**i)
        h = (-h) % mod
        h[0] = (h[0] + 2) %mod
        g = convolve2(g,h,2**i)
    return g[:limit]

def integral(f,limit):
    F = np.zeros(len(f),np.int64)
    F[1:] = f[:len(f)-1] * inv[1:len(f)]
    return (F % mod)[:limit]

def diff(f,limit):
    arange = np.array([i for i in range(len(f))],np.int64)
    res = (f * arange) % mod
    res = np.resize(res[1:],limit)
    res[-1] = 0
    return res[:limit]

def log(f,limit):
    res = convolve2(diff(f,limit),inverse(f,limit),limit)
    return integral(res,limit)

def exp(f,limit):
    l = len(f)
    L = 1<<((l-1).bit_length())
    n = L.bit_length()-1
    f = np.resize(f,L)
    f[L:] = 0

    res = np.array([1],np.int64)
    for i in range(1,n+1):
        res = np.resize(res,2**i)
        res[2**(i-1):] = 0
        g = log(res,2**i)
        h = (f[:2**i]-g[:2**i]) % mod
        h[0] = (h[0] + 1) % mod
        res = convolve2(res,h,2**i)
    return res[:limit]

def pow_poly(f,k,limit):
    l = len(f)
    L = 1<<((l-1).bit_length())
    n = L.bit_length()-1
    f = np.resize(f,L)
    f[L:] = 0

    g = (k * log(f,limit)) % mod
    h = exp(g,limit)

    return h[:limit]

def main():
    N,M = map(int,input().split())
    P = list(map(int,input().split()))

    P = [P[i]-1 for i in range(N)]
    cycle = []
    used = [False]*N
    for i in range(N):
        if not used[i]:
            used[i] = True
            c = 1
            pos = i
            while not used[P[pos]]:
                pos = P[pos]
                used[pos] = True
                c += 1
            cycle.append(c)

    n = len(cycle)

    binom_poly = np.array([1],np.int64)
    for i in range(n):
        c = [cmb(cycle[i],j,mod)*pow(-1,j,mod) for j in range(1,cycle[i]+1)]
        c[0] = (c[0] + 1) % mod
        c = np.array(c,np.int64)
        c %= mod
        binom_poly = convolve2(binom_poly,c,len(binom_poly)+cycle[i]-1)

    g = [g2[i+1] for i in range(N+1)]
    g = np.array(g,np.int64)
    res = pow_poly(g,M,N+1)
    res = (res * g2[M]) % mod
    res = list(res)
    Res = [0 for i in range(N+1)]
    for i in range(N-M+1):
        Res[i+M] = (res[i] * g1[i+M]) % mod
    res = Res

    ans = 0
    for j in range(N-n+1):
        ans += res[n+j] * binom_poly[j]
        ans %= mod
    print((ans*(-1)**N)%mod)

if __name__=="__main__":
    main()