AGC006

死んだ。
まためっちゃ面白いセットだった。sugim48すごすぎ。

A:流石にやるだけ
B:いきなり難しいconstructive. 実験したらrotate(ans+x-2,ans+x,ans+N) が見つかったのでこれを出力した。D考えてるときに真ん中でa<b<cってなってるなら絶対bが残ることに気づいたのでわざわざこんなことしなくてもよい.
C:期待値なのでそのまま計算していいことには気づいたのに、差分取るのに気づけなかった・・・ Air Pollution という問題がAOJにあるんですけど、それと同じ発想。なんで気づけなかったかなあ・・・
D:全く何もわからなかったんですが、にぶたんと言われたら解ける。こういうのなくしたいなあ(あまりにぶたんして嬉しい形に見えなかったけど、中央値→x以下が2個 みたいな言い換えを考えるとbool値にすると嬉しいことがわかる)
E:二個離れたところのswap+その二個と間の一個が逆向きになる という操作なんだけど、二つの列の転倒数のparityを見ればOKだった.まあなんかflipの場所を変えるのは出来そうだったので多分そうだろうと思ったけどそうなのか。
F:まだなんもわかってない

Cはともかく、D解けないのはダメなので気をつけないと

FunctionalEquation (SRM456hard)

おもしろかった.

問題:整数1≤C≤16 と, x[i],y[i]が与えられる.
f:Z->Zで,f(2f(x)-x+1)=f(x)+C を満たすという条件のもとでsigma |f(x[i])-y[i]| を最小化せよ.
x,yの長さは<=10000.値は10^9.

まずxに2f(x)-x+1を再代入しするとf(x+2C)=f(x)+2Cが得られる.
逆に,
f(x+2C)=f(x)+2C for any x かつ
f(2f(x)-x+1)=f(x)+C for any x in [0,2C)
が成り立てば元の式も成り立つことがわかるので,これを満たすようなものを考える.
f(0)~f(2C-1)を決めれば自動的にすべての値が決まる.
あるf(x)を決めると,xを代入してf(2f(x)-x+1)の値も決まる.そこから更に何か決まるかというと,できることがxに2f(x)-x+1を代入するしかないが,これはさっきやったf(x+2C)=f(x)+2Cになるのでやる意味がない.

冷静に考えると,xと2f(x)-x+1は(mod 2Cで)偶奇が異なるので, 偶数xに対して,f(x)=a と決めると,元の式から,y=2a-x+1として,f(y)=a+Cが決まる.それ以外(mod 2Cで)は決まらない.
偶数xがどの奇数yとペアを組むかはa mod Cに依存することがわかるので,これを固定する.この時,f(x)=aのとき,f(y)=b(=a+C)とすると,f(x)=a+?Cのときはf(y)=b-?Cになることがわかる,?には任意の整数が入れられる.この時関連する部分(つまりx[i]%2Cがxかy)での値を最小化するには,だいたい中央値っぽいところを?Cとして取れば良いので簡単に計算できる.

さいごにbitDPしておわり.

modが16*2で32になって焦るが,xと2f(x)-x+1の偶奇が違う, というところに気づくのが難しかった.
これ非数オリ勢だとはじめのFEパートどのくらいで気づくんだろう

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;
bool B(int x,int y){return (x>>y)&1;}
typedef long long ll;
const int MN=10000;
ll inf=1e18;
ll xs[MN],ys[MN];
ll cost[16][16];
ll dp[17][1<<16];
class FunctionalEquation{
	public:
	long long minAbsSum(int C, int N, int xzero, int xprod, int xadd, int xmod, int yzero, int yprod, int yadd, int ymod){
		xs[0]=xzero,ys[0]=yzero;
		rep(i,N-1) xs[i+1]=(xs[i]*xprod+xadd)%xmod,ys[i+1]=(ys[i]*yprod+yadd)%ymod;
		int M=C*2;
		rep(t,C){
			//f(x)=a+?*C,f(y)=b-?*C
			int x=t*2;
			rep(a,C){
				int yo=2*a-x+1;
				int y=(yo%M+M)%M;
				int b=a+C+(y-yo);
				vector<ll> vc;
				rep(i,N){	// |?*C - tmp|
					if(xs[i]%M==x){
						ll tmp=ys[i]+x-xs[i]-a;
						vc.pb(tmp);
					}
					if(xs[i]%M==y){
						ll tmp=-(ys[i]+y-xs[i]-b);
						vc.pb(tmp);
					}
				}
				sort(all(vc));
				int K=vc.size();
				if(K==0){
					cost[x/2][y/2]=0;
					continue;
				}
				ll med;
				if(K%2==0) med=(vc[K/2-1]+vc[K/2])/2;
				else med=vc[K/2];
				ll qo=med/C;
				ll mn=inf;
				for(ll q=qo-2;q<=qo+2;q++){
					ll sum=0;
					rep(i,K) sum+=abs(q*C-vc[i]);
					chmin(mn,sum);
				}
				cost[x/2][y/2]=mn;
//				printf("f(%d)=%d+?C,  f(%d)=%d+?C,   cost=%lld\n",x,a,y,b,mn);
			}
		}
//		rep(i,C) rep(j,C) printf("cost[%d][%d]=%lld\n",i,j,cost[i][j]);
		rep(i,C+1) rep(j,1<<C) dp[i][j]=inf;
		dp[0][0]=0;
		rep(i,C) rep(b,1<<C) if(dp[i][b]!=inf){
			rep(j,C) if(!B(b,j)){
				chmin(dp[i+1][b|(1<<j)],dp[i][b]+cost[i][j]);
			}
		}
		return dp[C][(1<<C)-1];
	}
};

yukicoder No.284 門松と魔法(2)

発想は良問.実装が大変.

問題:数列が与えられる,部分列Bであって,B[i]!=B[i+2] かつ (B[i]<B[i+1]>B[i+2] または B[i]>B[i+1]<B[i+2]) を任意のiに対して満たすようなものの長さの最大値を求めよ.

まずはじめB[i]!=B[i+2]を完全に見落としてたためO(N^2)DPから自然に生えるsegtreeを書いて提出したらWAだったので,問題文を読むと違うことがわかった.

このsegtreeでは最後がbで終わるような列 というのの最大値をsegで持っていたけど,最後から二番目aの情報も持って置かなければいけないことがわかる.でも,条件は,今xを更新しようとしてるとするとx!=a なので,aとしてはでかい方からたかだか2個しか持たなくて良い.
なので,seg[b]に,(mx,a)の組(aは異なる)を2つ持っておけばできる.(mergeはちょっと面倒だけど)

ただ,これを持とうとするとちゃんとqueryの際にどっから最大値が来たか,というのを求めないといけないので,idも持っておく.
しかも更新のためには2番目の最大値も求めなければいけないが,これは冷静に考えるともう一度1番目のとこを除いてクエリを投げれば良い.

実装がしんどかった(2番目がまだないとか,一箇所しか候補がないなら2番目はなしとか,細かい部分が大変).

はじめoperator+でmapを使ってたらTLEしたので,bubble sortにしたら通った.

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;
template<class S,class T> ostream& operator<<(ostream& o,const pair<S,T> &p){return o<<"("<<p.fs<<","<<p.sc<<")";}
template<class T> ostream& operator<<(ostream& o,const vector<T> &vc){o<<"sz = "<<vc.size()<<endl<<"[";for(const T& v:vc) o<<v<<",";o<<"]";return o;}

struct D{
	int segid[2],id[2],mx[2];			//id[i] -> segid[i] -> nxt
	D(){
		segid[0]=segid[1]=-1;
		id[0]=id[1]=-1;
		mx[0]=mx[1]=-1;
	}
	D(int x){
		segid[0]=x;
		segid[1]=-1;
		id[0]=-1;
		id[1]=-1;
		mx[0]=1;
		mx[1]=-1;
	}
	D(int a,int b,int c,int d,int e,int f){
		id[0]=a,segid[0]=b,mx[0]=c,id[1]=d,segid[1]=e,mx[1]=f;
	}
	const static D e;
/*	D operator+(const D& r) const {
		typedef pair<int,int> P;
		typedef pair<P,int> PP;
		map<int,P> mp;
		rep(i,2) if(mx[i]>=0) chmax(mp[id[i]],P(mx[i],segid[i]));
		rep(i,2) if(r.mx[i]>=0) chmax(mp[r.id[i]],P(r.mx[i],r.segid[i]));
		vector<PP> vp;
		for(auto it:mp) vp.pb(PP(it.sc,it.fs));
		sort(all(vp),greater<PP>());
		int K=vp.size();
		D ret;
		rep(i,min(2,K)) ret.id[i]=vp[i].sc,ret.mx[i]=vp[i].fs.fs,ret.segid[i]=vp[i].fs.sc;
		return ret;
	}*/
	D operator+(const D& r) const {
		int tid[4],tsegid[4],tmx[4],K=0;
		rep(i,2) if(mx[i]>=0) tid[K]=id[i],tsegid[K]=segid[i],tmx[K]=mx[i],K++;
		rep(i,2) if(r.mx[i]>=0) tid[K]=r.id[i],tsegid[K]=r.segid[i],tmx[K]=r.mx[i],K++;
		rep(i,K-1) rep(j,K-1-i) if(tmx[j]<tmx[j+1]) swap(tmx[j],tmx[j+1]),swap(tsegid[j],tsegid[j+1]),swap(tid[j],tid[j+1]);

		int I=0;
		D ret;
		rep(i,K){
			bool ok=1;
			rep(j,i){
				if(tid[j]==tid[i]){
					ok=0;
					break;
				}
			}
			if(!ok) continue;
			ret.id[I]=tid[i],ret.segid[I]=tsegid[i],ret.mx[I]=tmx[i],I++;
			if(I==2) break;
		}
		return ret;
	}
	void inc(){
		rep(i,2) if(mx[i]>=0) mx[i]++;
	}
	friend ostream& operator<<(ostream& o,const D& d){
		o<<"d=   mx="<<d.mx[0]<<" at ("<<d.id[0]<<","<<d.segid[0]<<")"<<endl;
		o<<"     mx="<<d.mx[1]<<" at ("<<d.id[1]<<","<<d.segid[1]<<")"<<endl;
		return o;
	}
};
const D D::e = D();

template<class D>
struct segtree{
	static const int N=1<<17;

	D e=D::e;
	vector<D> seg;
	segtree():seg(N*2,e){}
	void update(int k,D val){
		k+=N;
		seg[k]=val;
		k/=2;
		while(k){
			seg[k]=seg[k*2]+seg[k*2+1];
			k/=2;
		}
	}
	D calc(int a,int b,int l=0,int r=N,int k=1){
		if(b<=a||b<=l||r<=a) return e;
		if(a<=l&&r<=b) return seg[k];
		return calc(a,b,l,(l+r)/2,k*2)+calc(a,b,(l+r)/2,r,k*2+1);
	}
};
segtree<D> up,down;


int N;
int x[100000];
vector<int> xs;
int main(){
	cin>>N;
	rep(i,N){
		scanf("%d",x+i);
		xs.pb(x[i]);
	}
	sort(all(xs));
	xs.erase(unique(xs.begin(),xs.end()),xs.end());
	int K=xs.size();
	rep(i,N) x[i]=lower_bound(all(xs),x[i])-xs.begin();
	rep(i,N){
//		printf("i=%d\n",i);
		int X=x[i];
		D dval,uval;
		{//go up
			D dmx=down.calc(0,X);
			int mx=-2,arg;
			rep(j,2){
				if(dmx.id[j]!=X){
					mx=dmx.mx[j];
					arg=dmx.segid[j];
					break;
				}
			}
			if(mx>=0){
				D dmx2 = down.calc(0,arg) + down.calc(arg+1,X);
				int mx2=-2,arg2;
				rep(j,2){
					if(dmx2.id[j]!=X){
						mx2=dmx2.mx[j];
						arg2=dmx2.segid[j];
						break;
					}
				}
				dval = D(arg,X,mx,arg2,X,mx2);
			}else{
				dval = D(arg,X,mx,-1,-1,-1);
			}
			dval.inc();
			dval = dval + D(X);
		}

		{//go down
			D umx=up.calc(X+1,K);
			int mx=-2,arg;
			rep(j,2){
				if(umx.id[j]!=X){
					mx=umx.mx[j];
					arg=umx.segid[j];
					break;
				}
			}
			if(mx>=0){
				D umx2 = up.calc(X+1,arg) + up.calc(arg+1,K);
				int mx2=-2,arg2;
				rep(j,2){
					if(umx2.id[j]!=X){
						mx2=umx2.mx[j];
						arg2=umx2.segid[j];
						break;
					}
				}
				uval = D(arg,X,mx,arg2,X,mx2);
			}else{
				uval = D(arg,X,mx,-1,-1,-1);
			}
			uval.inc();
			uval = uval + D(X);
		}
		// cout<<"dval"<<endl<<dval<<endl;
		// cout<<"uval"<<endl<<uval<<endl;
		up.update(X,dval);
		down.update(X,uval);
	}
	int ans=max(up.calc(0,K).mx[0],down.calc(0,K).mx[0]);
	if(ans<=2) ans=0;
	cout<<ans<<endl;
}

天下一2016本選

ほんとうにつらい・・・
C,Dをやって,Eでcinを使ってTLEをした後,FをスルーしてAとBに向かってしまい,結局BはN≤10までしか解かないコードを,Aは996頂点0辺の1彩色可能なグラフをsubmitして終了してしまった.
F,得意系だったのでスルーしてしまったのは本当に残念.次は絶対日和らないようにします・・・
F:

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;
template<class S,class T> ostream& operator<<(ostream& o,const pair<S,T> &p){return o<<"("<<p.fs<<","<<p.sc<<")";}
template<class T> ostream& operator<<(ostream& o,const vector<T> &vc){o<<"sz = "<<vc.size()<<endl<<"[";for(const T& v:vc) o<<v<<",";o<<"]";return o;}
typedef long long ll;
int N,G,P[36];
ll dp[37][37][1297];
ll C[37][37];
int main(){
	rep(i,37){
		C[i][0]=C[i][i]=1;
		for(int j=1;j<i;j++) C[i][j]=C[i-1][j-1]+C[i-1][j];
	}
	cin>>N>>G;
	rep(i,N) cin>>P[i];
	int sum=0;
	rep(i,N) sum+=P[i];
//	int S=N*G;
	int S=sum-N*G;
	if(S<=0){
		cout<<sum<<endl;
		return 0;
	}
	double ans=0;
	rep(last,N){
		rep(i,N+1) rep(j,N+1) rep(s,1297) dp[i][j][s]=0;
		dp[0][0][0]=1;
		rep(i,N){
			rep(j,N){
				rep(s,1297){
					dp[i+1][j][s]+=dp[i][j][s];
					if(i!=last && s+P[i]<=1296) dp[i+1][j+1][s+P[i]]+=dp[i][j][s];
				}
			}
		}
		double ex=0;
		for(int j=0;j<N;j++){
			ex+=(double)N/(N-j);
			for(int s=max(0,S-P[last]);s<S;s++){
/*				if(dp[N][j][s]){
					printf("dp[N][%d][%d]=%I64d\n",j,s,dp[N][j][s]);
					show(ex*G+(sum-s-P[last]));
					show((double)dp[N][j][s]/C[N][j+1]/(j+1));
				}*/
				ans+=(ex*G+(sum-s-P[last])) * ((double)dp[N][j][s]/C[N][j+1]/(j+1));
			}
		}
	}
	printf("%.12f\n",ans);
}

JAG2016夏合宿

に参加しました

day1

補習に行って、なぜか自分のPC上の仮想環境(VMware)が立ち上がらなくなっていることに気づく。そこでいろいろやってたら遅刻した(ごめんなさい)。sugimさんがいないので自己紹介でsugimを名乗りかける。談話室でボドゲをやる。The Bossが初見だったけど面白かった。寝る。

day2

5:30くらいにクッソまぶしい朝の光が窓から刺さってきたので起きたら本当にパソコンが壊れていて焦る(windowsが起動しない)。とりあえず朝食を解いてふて寝した。コンテスト直前にカスタマーサービスに連絡したらどうやらwindows 10 anniversary updateの不具合だったらしくBIOSの設定をちょっと変えたらなおった。こわすぎ。

コンテスト1: 2完差でpoに負けて2位。やばすぎ。ていうかやるだけしか解けなかった・・・。うちのチームだとある程度の考察とある程度の実装みたいなのをなぎ倒すのが得意なのでこういうセットが一番相対順位上がるのかなあ(ダメ)

万豚記に行くのを談話室で待ってたらもう出発してたのでこの日は行かなかった。

day3

ちゃんとカーテンを閉めておいたのでゆっくり起きれた。移動してヒカリエに行く。景色が良かった。レゴブロックで遊べる場所があって、人々がレゴブロックが入っている穴をレゴブロックで蓋をして遊んでいた(最終的に完成していた)。

めっちゃくつろげる空間で模擬地区。
わふれるかが寝坊して遅刻していた。

コンテスト2:Eでジャッジがバグっていたのとwafrelkaの遅刻で初動がうまくいかず、7完4位だった。取り組めば解けそうな問題が多くて良いセットだったと思うので残念(最後F,I,Kの3並列だったし流石にやばい)

AGC004があった。writerはsugim48(理由:4は偶数)だった。A~Dを倒した後E,Fが全くわからず48位。ほんとにクッソ良問セットでびっくりした。

コンテスト後談話室で駄弁る。「エセ芸術家」というゲームをやったがクッソ面白かった。全体的にみんなエセ芸術家の推理能力を過信しすぎていて、適当なのを書いてエセ芸術家が誰かわからなくなることが多かった。
皆もあててみてください(ほとんど無理だと思います) 記事の最後に答えを載せておきます。

day4

ギリギリだったがルームメイトと朝飯を食っていったら本当にギリギリになった、すみません(提案したのは僕です)。
今回はSleepの残り二人は参加していなかったため、snukeを誘ってチームを組んだ。チーム名はigma425(いつもの)

コンテスト3: Open Cupの過去問。ICPC系面倒と考察をsnukeがやっていたため、残り時間を使って面倒DPの実装をしていたらコンテストが終わった。10完でその場では一位だった(が、virtualでやばい人がいっぱいいた 18位)

終わった後解散だったが、snukeとsugimさんと万豚記へ向かう。おいしかった。ずっとクソなぞなぞを生成していた。

まとめ

楽しかったです。JAGのセットも割と好きだったし、準備もお疲れ様でした。

お絵かきの答え

食べ物: たいやき
動物: アリクイ
ことわざ: 石の上にも三年
たてもの: 橋
機械: 掃除機
競技プログラマ: わふれるか

CF AIM Tech Round 3 (Div. 1) (708)

A:場合分け はじめのaじゃないとこからaじゃないとこまでだけどaaa->aazに注意
B:場合分け 0がない時と1がない時とそれ以外で,それ以外だと0と1の数がわかるので全体と合ってることを確認した後は適切に0001111からずらす(計算して).
C:全方向木DP.dp[v]=v以下でN/2以下のサイズの部分木の内でもっとも大きい物 をもつとvをcentroidにしたければvの子uのうちN/2よりでかいものからできるだけでかく(つまりdp[u])を切り取ってvとつなげばいいのでこれはdpを計算すればわかる あとは全方向にすれば終わり.
コンテスト後にすぬけが言ってたように,そもそも重心を根として木DPすれば,切るのは重心をまたいだ部分なので楽に実装できるっぽい.
D:min_cost_flow.
元々a->bにf流れていてmaxがcっていうのを,

		if(f<=c){
			sum+=f;
			add_edge(a,T,f,0);
			add_edge(S,b,f,0);
			add_edge(b,a,f,1);
			add_edge(a,b,c-f,1);
			add_edge(a,b,INF,2);
		}else{
			ans+=f-c;
			sum+=c;
			add_edge(a,T,c,0);
			add_edge(S,b,c,0);
			add_edge(b,a,c,1);
			add_edge(a,b,f-c,0);
			add_edge(a,b,INF,2);
		}

こうする.(add_edgeはfrom,to,capacity,cost)
SとTっていう新たな頂点を追加して,

• f≤c の時

すでにf流してある というのは, a->TとS->bにf流さなければいけない という設定にして,f* ≤ f の時, (fを減らす) コスト1でb->aに流し戻せる(ただしfまで)
f* > f の時, (fを増やす) cまで(つまり新たなグラフでc-fまで)はコスト1で増やせて,それ以降はfとcを共に増やす必要が有るためコスト2になる.

• f>c の時

まず絶対にf-cのコストは掛かるため,これを答えに足しあわせておく.そして今c流してるという設定にしておく.前述のf-cのコストの範囲で,c流す~f流すまでは自由にできる(cを減らし,fを増やすことで)ので,f-c分はコスト0で流せる.さらに追加分はコスト2.fを減らすのはコスト1でcまで減らせる.
sumでS->Tに流さねばならない分を計算していて,これでS->Tにsum流す必要があって,それとは別に0からN-1に自由に流して良くて,これはN-1から0へ辺を貼ればOK.その後min_cost_flow(S,T,sum)を呼べばおわり.

E:dp[i][l][r]=i段目が[l,r]でまだ連結な期待値 みたいなのを持つと更新をうまくやればO(W^2 * H)でできる.更新をうまくやるのにlsum[i][l]=sigma dp[i][x][l] みたいなのを持つんだけど,実はdpはいらずにlsum,rsumだけ持って更新ができて,この更新もO(1)でできるためO(HW)になり終了 っぽいけど面倒だし実装してない.

コンテスト:
C開いたらモロ全方向木DPで笑ってしまった(ついこないだ書いたばっかり) やる→Bで1WA→Aやる→E考える→順位表見る→Dが結構通されている→D通す みたいな感じだった.
Bが落ちてキレた.

A:Submission #20123051 - Codeforces
B:Submission #20134838 - Codeforces
C:Submission #20116938 - Codeforces
D:Submission #20131523 - Codeforces

TreeDistance (TCO14 Round 3B d1 hard)

問題:N頂点の木Tが与えられます.「辺を一つ除いて,新たに付け加えて木にする」 という操作をK回以下出来るとき何種類の木が出来るでしょう?N,K≤50

まず, TからK回以下で木Xに出来る⇔TとXの両方に含まれる辺がN-1-K本以上(→は自明,左はいつかのこどふぉEに構成する問題が出た(
CF345 - sigma425のブログ).

Tの辺集合をEとして,Eのsubset Sに対して,SがTとXの両方に含まれる ようなものの個数を考えると,Sの連結成分の個数をC個,その大きさをn1,n2,..nCとすると,N^(C-2) * Πni になる.
これは非自明なんだけど,Prufer Sequenceと木との対応を知っていれば,そこからCayleyの式N^(N-2)が出てくることとの一般化として,これがわかる(sequenceに出てくる回数+1が次数で,連結成分ciとつなぐときはつなぎ方がni通りあるため,結局因数分解すると(n1+n2+...)^(C-2) * (+1回分の係数) となる.)

なので後は木DPで,dp[v][x][n]=v以下でx本を明確には固定せずに,今連結成分の大きさがn とやって,あとは包除.

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;
typedef long long ll;
ll mod=1e9+7;
void add(ll &x,ll y){
	x+=y;
	x%=mod;
}
ll C[51][51];
vector<int> G[50];
int N;
int sz[50];
ll dp[50][51][51];
ll tmp[51][51][51];
void dfs(int v){
	sz[v]=1;
	for(int u:G[v]) dfs(u),sz[v]+=sz[u];
	rep(k,G[v].size()+1) rep(i,51) rep(j,51) tmp[k][i][j]=0;
	tmp[0][0][1]=1;
	int S=1;
	int K=G[v].size();
	rep(i,K){
		int u=G[v][i];
		rep(x,S) rep(n,S+1) if(tmp[i][x][n]){
			rep(ax,sz[u]+1) rep(an,sz[u]+1) add(tmp[i+1][x+ax][n+an],tmp[i][x][n]*dp[u][ax][an]);
		}
		S+=sz[u];
	}
	rep(x,S) rep1(n,S){
		dp[v][x][n]=tmp[K][x][n];
		add(dp[v][x+1][0],tmp[K][x][n]*n);
	}
//	rep(x,S+1) rep(n,S+1) printf("dp[%d][%d][%d]=%lld\n",v,x,n,dp[v][x][n]);
}
ll Co(int x,int y){
	if(x==y) return 1;
	return C[x][y];
}
ll ex(ll x,int p){
	if(p==-1){
		ll a=1;
		rep(i,x+1){
			if(a%N==0){
				a/=N;
				break;
			}
			a+=mod;
		}
		return a;
	}
	ll a=1;
	while(p){
		if(p%2) a=a*x%mod;
		x=x*x%mod;
		p/=2;
	}
	return a;
}
class TreeDistance{
	public:
	int countTrees(vector <int> p, int K){
		rep(i,51){
			C[i][0]=C[i][i]=1;
			for(int j=1;j<i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
		}
		N=p.size()+1;
		rep(i,N-1) G[p[i]].pb(i+1);
		dfs(0);
		K=max(0,N-1-K);
		ll ans=0;
		for(int i=K;i<=N-1;i++){
			ll ad=dp[0][N-i][0]*ex(N,N-i-2)%mod*Co(i-1,K-1)%mod;
			if((i-K)%2==1) ad=mod-ad;
			add(ans,ad);
		}
		return ans;
	}
};