#include <iostream>#include <sstream>#include <iomanip>#include <utility>#include <vector>#include <list>#include <string>#include <stack>#include <queue>#include <deque>#include <set>#include <map>#include <algorithm>#include <functional>#include <numeric>#include <bitset>#include <complex>#include <cstdio>#include <cstring>#include <cmath>#include <cstdlib>#include <ctime>#include <climits>using namespace std;#define lp for(;;)#define repf(i,a,b) for (int i=(a);i<(b);++i)#define rrepf(i,a,b) for (int i=(a)-1;i>=(b);--i)#define rep(i,n) repf(i,0,n)#define rrep(i,n) rrepf(i,n,0)#define ft(i,a,b) for (int i=(a);i<=(b);++i)#define fdt(i,a,b) for (int i=(a);i>=(b);--i)#define for_nonempty_subsets(subset,set) for (int subset=set;subset;subset=(subset-1)&(set))#define for_in_charset(i,charset) for (cstr i=(charset);*i;i++)#define whl while#define rtn return#define fl(x,y) memset((x),char(y),sizeof(x))#define clr(x) fl(x,char(0))#define cpy(x,y) memcpy(x,y,sizeof(x))#define sf scanf#define pf printf#define vec vector#define pr pair#define que queue#define prq priority_queue#define itr iterator#define x first#define y second#define pb push_back#define mp make_pair#define ins insert#define ers erase#define lb lower_bound#define ub upper_bound#define rnk order_of_key#define sel find_by_order#define sz(x) (int((x).size()))#define all(x) (x).begin(),(x).end()#define srt(x) sort(all(x))#define uniq(x) srt(x),(x).erase(unique(all(x)),(x).end())#define rev(x) reverse(all(x))#define shf(x) random_shuffle(all(x))#define nxtp(x) next_permutation(all(x))typedef long long int lli;typedef double db;typedef const char* cstr;typedef string str;typedef vec<int> vi;typedef vec<vi> vvi;typedef vec<lli> vl;typedef vec<vl> vvl;typedef vec<bool> vb;typedef vec<vb> vvb;typedef vec<char> vc;typedef vec<vc> vvc;typedef vec<str> vs;typedef pr<int,int> pii;typedef pr<lli,lli> pll;typedef pr<db,db> pdd;typedef vec<pii> vpii;typedef vec<pll> vpll;typedef vec<pdd> vpdd;typedef map<int,int> mii;typedef map<str,int> msi;typedef map<char,int> mci;typedef set<int> si;typedef set<str> ss;typedef que<int> qi;int oo=(~0u)>>1;lli ooll=(~0ull)>>1;db inf=1e+10;db eps=1e-10;db gam=0.5772156649015328606;db pi=acos(-1.0);int dx[]={1,0,-1,0,1,-1,-1,1,0};int dy[]={0,1,0,-1,1,1,-1,-1,0};int MOD=1000000007;template<typename type>inline bool cmax(type& a,const type& b){rtn a<b?a=b,true:false;}template<typename type>inline bool cmin(type& a,const type& b){rtn b<a?a=b,true:false;}template<typename type>inline type sqr(const type& x){rtn x*x;}template<typename type>inline type mod(const type& x){rtn x%MOD;}inline int sgn(const db& x){rtn (x>+eps)-(x<-eps);}inline int dbcmp(const db& a,const db& b){rtn sgn(a-b);}template<typename type>inline pr<type,type> operator-(const pr<type,type>& x){rtn mp(-x.x,-x.y);}template<typename type>inline pr<type,type> operator+(const pr<type,type>& a,const pr<type,type>& b){rtn mp(a.x+b.x,a.y+b.y);}template<typename type>inline pr<type,type> operator-(const pr<type,type>& a,const pr<type,type>& b){rtn mp(a.x-b.x,a.y-b.y);}template<typename type>inline pr<type,type> operator*(const pr<type,type>& a,const type& b){rtn mp(a.x*b,a.y*b);}template<typename type>inline pr<type,type> operator/(const pr<type,type>& a,const type& b){rtn mp(a.x/b,a.y/b);}template<typename type>inline pr<type,type>& operator-=(pr<type,type>& a,const pr<type,type>& b){rtn a=a-b;}template<typename type>inline pr<type,type>& operator+=(pr<type,type>& a,const pr<type,type>& b){rtn a=a+b;}template<typename type>inline pr<type,type>& operator*=(pr<type,type>& a,const type& b){rtn a=a*b;}template<typename type>inline pr<type,type>& operator/=(pr<type,type>& a,const type& b){rtn a=a/b;}template<typename type>inline type cross(const pr<type,type>& a,const pr<type,type>& b){rtn a.x*b.y-a.y*b.x;}template<typename type>inline type dot(const pr<type,type>& a,const pr<type,type>& b){rtn a.x*b.x+a.y*b.y;}template<typename type>inline type gcd(type a,type b){if(b)whl((a%=b)&&(b%=a));rtn a+b;}template<typename type>inline type lcm(type a,type b){rtn a*b/gcd(a,b);}inline lli bin_pow(lli x,lli y){lli z=1;whl(y){if(y&1)z=mod(z*x);x=mod(sqr(x)),y>>=1;}rtn z;}template<typename istream,typename first_type,typename second_type>inline istream& operator>>(istream& cin,pr<first_type,second_type>& x){rtn cin>>x.x>>x.y;}template<typename ostream,typename first_type,typename second_type>inline ostream& operator<<(ostream& cout,const pr<first_type,second_type>& x){rtn cout<<x.x<<" "<<x.y;}template<typename istream,typename type>inline istream& operator>>(istream& cin,vec<type>& x){rep(i,sz(x))cin>>x[i];rtn cin;}template<typename ostream,typename type>inline ostream& operator<<(ostream& cout,const vec<type>& x){rep(i,sz(x))cout<<x[i]<<(i+1==sz(x)?"":" ");rtn cout;}inline ostream& pdb(int prcs,db x){rtn cout<<setprecision(prcs)<<fixed<<(sgn(x)?(x):0);}template<typename type>inline void bit_inc(vec<type>& st,int x,type inc){whl(x<sz(st))st[x]+=inc,x|=x+1;}template<typename type>inline type bit_sum(const vec<type>& st,int x){type s=0;whl(x>=0)s+=st[x],x=(x&(x+1))-1;rtn s;}template<typename type>inline type bit_kth(const vec<type>& st,int k){int x=0,y=0,z=0;whl((1<<(++y))<=sz(st));fdt(i,y-1,0){if((x+=1<<i)>sz(st)||z+st[x-1]>k)x-=1<<i;else z+=st[x-1];}rtn x;}inline void make_set(vi& st){rep(i,sz(st))st[i]=i;}inline int find_set(vi& st,int x){int y=x,z;whl(y!=st[y])y=st[y];whl(x!=st[x])z=st[x],st[x]=y,x=z;rtn y;}inline bool union_set(vi& st,int a,int b){a=find_set(st,a),b=find_set(st,b);rtn a!=b?st[a]=b,true:false;}inline void make_set(vpii& st){rep(i,sz(st))st[i]=mp(i,1);}inline int find_set(vpii& st,int x){int y=x,z;whl(y!=st[y].x)y=st[y].x;whl(x!=st[x].x)z=st[x].x,st[x].x=y,x=z;rtn y;}inline bool union_set(vpii& st,int a,int b){a=find_set(st,a),b=find_set(st,b);rtn a!=b?(st[a].y>st[b].y?st[a].x=b,st[a].y+=st[b].y:st[b].x=a,st[b].y+=st[a].y),true:false;}template<typename type>inline void merge(type& a,type& b){if(sz(a)<sz(b))swap(a,b);whl(sz(b))a.ins(*b.begin()),b.ers(b.begin());}#pragma comment(linker, "/STACK:31457280")const int MAXN=100000;lli inv[MAXN];int n,m;int par[MAXN];vi adj[MAXN];int tv;int vis[MAXN];int onc[MAXN];vi cir;vec<pr<pll,lli> > info;vl f;int rpt;lli rptf;lli C(lli n,lli m){if (n<m) rtn 0;else{lli ans=inv[m];ft(i,n-m+1,n) ans=(ans*i)%MOD;rtn ans;}}pr<pll,lli> calc(int u){vec<pr<pll,lli> > lst;rep(i,sz(adj[u])){int v=adj[u][i];if (onc[v]!=tv) lst.pb(calc(v));}srt(lst);pr<pll,lli> ans=mp(mp(lli(1),lli(1)),1);rep(i,sz(lst)){ans.x.y=(ans.x.y*bin_pow(2,lst[i].x.x)+lst[i].x.y)%MOD;ans.x.x+=lst[i].x.x;}ans.x.y=(ans.x.y*2+0)%MOD;ans.x.x+=1;rep(i,sz(lst)){int b=i,e=i+1;whl(e<sz(lst)&&lst[e]==lst[b]) e++;lli ttl=e-b;lli cnt=(lli(m-1)*lst[b].y)%MOD;ans.y=(ans.y*C(cnt+ttl-1,ttl))%MOD;i=e-1;}rtn ans;}int match(){vi pi(sz(info));{pi[0]=-1;int j=-1;repf(i,1,sz(info)){whl(j!=-1&&info[j+1]!=info[i]) j=pi[j];if (info[j+1]==info[i]) j++;pi[i]=(j!=-1&&i+1<sz(info)&&info[j+1]==info[i+1])?pi[j]:j;}}{int j=-1;repf(i,1,sz(info)*2){whl(j!=-1&&info[j+1]!=info[i%sz(info)]) j=pi[j];if (info[j+1]==info[i%sz(info)]) j++;if (j==sz(info)-1) rtn i-j;}}}int main(){lli fac=1;repf(i,1,MAXN) fac=(fac*i)%MOD;inv[MAXN-1]=bin_pow(fac,MOD-2);rrep(i,MAXN-1) inv[i]=(inv[i+1]*(i+1))%MOD;whl(~sf("%d%d",&n,&m)){rep(i,n) sf("%d",&par[i]),--par[i],adj[par[i]].pb(i);int u;++tv;vis[u=0]=tv;whl(cmax(vis[u=par[u]],tv));onc[u]=tv,cir.pb(u);whl(cmax(onc[u=par[u]],tv)) cir.pb(u);rep(i,sz(cir)) info.pb(calc(cir[i]));f=vl(sz(cir)+1);f[1]=0;lli g=m;ft(i,2,sz(cir)){g=(g*(m-1))%MOD;f[i]=(g-f[i-1])%MOD;}rpt=match();rptf=1;rep(i,rpt) rptf=(rptf*info[i].y)%MOD;lli ans=0;int tim=sz(cir)/rpt;ft(i,1,tim){int j=gcd(tim,i);ans=(ans+bin_pow(rptf,j)*f[rpt*j])%MOD;}ans=(ans*bin_pow(tim,MOD-2))%MOD;if (ans<0) ans+=MOD;pf("%dn",int(ans));rep(i,n) adj[i].clear();cir.clear();info.clear();}}
