Description：

给一棵内向基环树，求本质不同的染色方案数。

$1\leqslant n\leqslant 10^5$

Solution：

首先考虑如果是树怎么做，如果这个点有很多棵子树，那么我们可以把每种染色方式看成一个染色，但是有一个问题就是子树可能会同构，我们就用树哈希的方式把所有子树找出来，然后对于每一个子树的等价类，相当于是一个物品相同盒子不同盒子可以空乒乓球装箱问题，可以用插板法解决，让后我们就得出了每棵子树的方案数，然后再来考虑环，环和子树的区别就是子树可以随便换但是环只能旋转，于是我们还是把每一种染色方案都看成一个染色然后在环上做 $polya$ 定理就行了。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<queue>

#include<cctype>

#include<cstring>

using namespace std;

inline int rd()

{

	register int res = 0,f = 1;register char c = getchar();

	while(!isdigit(c)){if(c == '-')f = -1;c = getchar();}

	while(isdigit(c))res = (res << 1) + (res << 3) + c - '0',c = getchar();

	return res * f;

}

int gcd(int a,int b){return (b == 0 ? a : gcd(b,a % b));}

int n,m;

#define MAXN 100010

int fa[MAXN];

#define MOD 1000000007

#define MAX 200010

#define N 200000

int fac[MAX],inv[MAX];

int power(int a,int b)

{

	int res = 1;

	for(;b > 0;b = b >> 1,a = 1ll * a * a % MOD)if(b & 1)res = 1ll * res * a % MOD;

	return res;

}

struct edge

{

	int to,nxt;

}e[MAXN << 1];

int edgenum = 0;

int lin[MAXN] = {0};

int ind[MAXN];

void add(int a,int b)

{

	++ind[b];

	e[++edgenum] = (edge){b,lin[a]};lin[a] = edgenum;

	e[++edgenum] = (edge){a,lin[b]};lin[b] = edgenum;

	return;

}

bool incir[MAXN];

int siz[MAXN];

struct hash

{

	#define V1 1001001

	#define V2 19260817

	#define M1 1000000007

	#define M2 998244353

	int v1,v2,f;

	friend bool operator < (hash a,hash b)

	{

		return (a.v1 == b.v1 ? (a.v2 == b.v2 ? a.f < b.f : a.v2 < b.v2) : a.v1 < b.v1);

	}

	friend bool operator == (hash a,hash b)

	{

		return (a.v1 == b.v1 && a.v2 == b.v2);

	}

}f[MAXN];

hash s[MAXN];

int top = 0;

int C(int n,int m)

{

	if(n < m)return 0;

	if(n <= N)return 1ll * fac[n] * inv[m] % MOD * inv[n - m] % MOD;

	else

	{

		int ans = inv[n - m];

		for(int i = m + 1;i <= n;++i)ans = 1ll * ans * i % MOD;

		return ans;

	}

}

int calc(int n,int m)

{

	return C(n + m - 1,m - 1);

}

void dp(int k,int fa)

{

	siz[k] = 1;

	for(int i = lin[k];i != 0;i = e[i].nxt)

	{

		if(e[i].to == fa || incir[e[i].to])continue;

		dp(e[i].to,k);

		siz[k] += siz[e[i].to];

	}

	top = 0;

	for(int i = lin[k];i != 0;i = e[i].nxt)

	{

		if(e[i].to == fa || incir[e[i].to])continue;

		s[++top] = f[e[i].to];

	}

	sort(s + 1,s + 1 + top);

	f[k].v1 = 1;f[k].v2 = 1;

	for(int i = 1;i <= top;++i)

	{

		f[k].v1 = (1ll * f[k].v1 * V1 % M1 + s[i].v1) % M1;

		f[k].v2 = (1ll * f[k].v2 * V2 % M2 + s[i].v2) % M2;

	}

	s[top + 1].v1 = -1;s[top + 1].v2 = -1;

	f[k].f = m % MOD;

	int cnt = 0;//cout << " : " << k << endl;cout << "st" << endl;

	for(int i = 1;i <= top;++i)

	{

		if(!(s[i] == s[i + 1]))

		{

			++cnt;//cout << cnt << " -> " << s[i].f << endl;

			f[k].f = 1ll * f[k].f * calc(cnt,s[i].f) % MOD;

			cnt = 0;

		}

		else ++cnt;

	}//cout << "ed" << endl;

	f[k].v1 = 1ll * f[k].v1 * siz[k] % MOD;

	f[k].v2 = 1ll * f[k].v2 * siz[k] % MOD;

	return;

}

int cir[MAXN];

int nxt[MAXN];

int main()

{

	fac[0] = 1;for(int i = 1;i <= N;++i)fac[i] = 1ll * fac[i - 1] * i % MOD;

	inv[N] = power(fac[N],MOD - 2);for(int i = N - 1;i >= 0;--i)inv[i] = 1ll * inv[i + 1] * (i + 1) % MOD;

	scanf("%d%d",&n,&m);

	for(int i = 1;i <= n;++i)add(i,fa[i] = rd());

	queue<int> q;

	for(int i = 1;i <= n;++i)if(ind[i] == 0)q.push(i);

	while(!q.empty())

	{

		int k = q.front();q.pop();

		for(int i = lin[k];i != 0;i = e[i].nxt)

		{

			if(ind[e[i].to] == 0)continue;

			if((--ind[e[i].to]) == 0)q.push(e[i].to);

		}

	}

	for(int i = 1;i <= n;++i)if(ind[i] != 0)incir[i] = true,++cir[0];

	for(int i = 1;i <= n;++i)if(incir[i])dp(i,0);

	for(int i = 1;i <= n;++i)

	{

		if(incir[i])for(int k = 1,cur = i;k <= cir[0];++k,cur = fa[cur])cir[k] = cur;

	}

	for(int i = 2,j = 0;i <= cir[0];++i)

	{

		while(j && !(f[cir[i]] == f[cir[j + 1]]))j = nxt[j];

		if(f[cir[i]] == f[cir[j + 1]])++j;

		nxt[i] = j;

	}

	int len = nxt[cir[0]];

	while(cir[0] % (cir[0] - len) != 0)len = nxt[len];

	len = cir[0] - len;

	int val = 1;

	for(int i = 1;i <= len;++i)val = 1ll * val * f[cir[i]].f % MOD;

	int tot = cir[0] / len;

	int ans = 0;

	for(int i = 1;i <= tot;++i)

	{

		int c = gcd(i,tot);

		ans = (ans + power(val,c)) % MOD;

	}

	ans = 1ll * ans * power(tot,MOD - 2) % MOD;

	cout << ans << endl;

	return 0;

}