mod = 998244353
omega = pow(3,119,mod)
rev_omega = pow(omega,mod-2,mod)

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

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

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

inv_list = [inv[i] for i in range(N+1)]

mod = 998244353

import numpy as np

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]

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)

a = [0]*(N+1)
for c in cycle:
    a[c-1] += 1

f = [0 for i in range(N+1)]
for i in range(1,N+1):
    if not a[i]:
        continue
    for j in range(1,N//i+1):
        f[i*j] += -a[i] * inv_list[j]
        f[i*j] %= mod

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

f = exp(f,N+1)

nf = np.array([0]*(N+1),np.int64)
nf[n:] = f[:N+1-n]
nf = (pow(-1,n)*nf) % mod
f = list(map(int,nf))

g = [(f[N-i] * g1[N-i]) % mod for i in range(N+1)]
e_x = [g2[i] for i in range(N+1)]

g = np.array(g,np.int64)
e_x = np.array(e_x,np.int64)

g_e_x = list(map(int,convolve2(g,e_x,N+1)))

f = [(g_e_x[N-i] * g2[i]) % mod for i in range(N+1)]

for i in range(1,N+1,2):
    f[i] = (-f[i]) % mod

f = [f[i+n] for i in range(N-n+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
poly_Stirling = Res
ans = 0
for j in range(N-n+1):
    ans += poly_Stirling[n+j] * f[j]
    ans %= mod

print((ans*(-1)**N)%mod)