Kawazoe(@riverplus)さんのCodeIqの問題です。
問題とテストケースはこちらにお世話になりました。解答もあります。
->https://blog.goo.ne.jp/r-de-r/e/f2b057442e37b97f008680d1eb1c087c
こちらにも情報があります。
->https://togetter.com/li/1105901
辺がnの正方形にm個の駒を置く場合、はぐれルークがk個となる配置は、
一辺nの正方形のはぐれルークとなるk個の配置のみの組み合わせ数に、
一辺n-kの正方形のはぐれルークがないm-k個の駒の配置の場合の数をかけたもの。
一方一辺nの正方形のはぐれルークにならないh個の駒の配置の場合の数は、
一辺nの正方形にh個の駒を置くすべての組み合わせから、
はぐれルークがある配置の場合の数を除いたもの。
したがって再帰の形になります。
Juliaで書きました。
binomial(n,r)はnCrです。
一辺n0の正方形にはぐれルークをmi個配置する組み合わせの数が
n0Cmi*n0Pmi=binomial(n0,mi)*div(factorial(n0),factorial(n0-mi))
一辺nの正方形にはぐれルークがないようにh個の駒を並べる場合の数がrec(n,h)
n0-sqrt(m-mi)>=miのsqrt(m-mi)
はm-mi個の駒を縦横に並べられる最小の正方形の一辺
（投稿時私の環境ではideoneでJavaしか表示されず実行できません）
#ver 1.53, 1.73,1.10.3
function solve(n0,m0,a)
function rec(n,h)
if h == 1 return 0
elseif h == n^2 || h == 0
return 1
else
sum=binomial(n^2,h)
i =1
while h >= i && n - sqrt(h-i) >= i
sum -= binomial(n,i) * div(factorial(n),factorial(n-i)) * rec(n-i,h-i)
i += 1
end
return sum
end
end
msum = 0
for m in 1 : m0
tmp = 0
mi = 1
while m >= mi && n0 - sqrt(m-mi) >= mi
tmp += binomial(n0,mi) * div(factorial(n0),factorial(n0-mi)) * rec(n0-mi,m-mi) * mi
mi += 1
end
msum += tmp / binomial(n0^2,m)
end
ans=Int(floor(msum * 1000))
s = a == ans ? "ok " : "no "
println(s,ans)
end
solve(2, 2, 1666)
solve(2, 3, 1666)
solve(3, 3, 2642)
solve(3, 4, 2928)
solve(4, 3, 3228)
solve(4, 4, 3967)
solve(4, 5, 4428)
solve(4, 7, 4797)
solve(4, 16,4857)
Prologで別解を書きました。
一辺nの正方形にm個の駒をならべ、はぐれ駒がk個の場合の数をf(n,m,k)としますと、
f(n,m,k)=f(n-1,m-1,k-1)n^2/kとなります。
(辺を一つ足してk個めの駒を置く配置がn^2ありますが、同じ状態がk個ずつ生じる)
そこで
f(n,m,k)=f(n-k,m-k,0)*n^2*(n-1)^2*...*(n-k+1)^2/(k(k-1)*..*1)
=f(n-k,m-k,0)*(n!)^2/((n-k)!)^2/k!
rec/4で上記の式を扱い、loop3/5では
f(n,m,0)=binomial(n^2,m)-(f(n,m,1)+...+f(n,m,m))
の部分を扱っています。
loop1/5がmのループ、loop2/5がkのループです。
ただし上と比べて早いわけではありません。
（投稿時私の環境ではideoneでJavaしか表示されず実行できません）
%swi-Prolog ver 7.4.2,7.6.4,9.2.5
%:-initialization(start).%(Online Prolog Compiler.)
% start.
ncr(_,0,1).
ncr(N,1,R):-R is N.
ncr(N,M,R):-N1 is N-1,M1 is M-1,ncr(N1,M1,R1),R is R1*N div M,!.
fact(0,1):-!.
fact(N,R):-fact(N,1,R).
fact(1,R,R):-!.
fact(N,R1,R):-R2 is N*R1,N1 is N-1,fact(N1,R2,R).
/*
2024.8.3解説に合うように変更
loop3(N,M,I,R,R):-(M<I;N<I),!.
loop3(N,M,I,R1,R):-
rec(N,M,I,R2),R3 is R1+R2,I1 is I+1,loop3(N,M,I1,R3,R).
rec(N,M,0,R):-N1 is N^2,ncr(N1,M,X),loop3(N,M,1,0,R1),R is X-R1,!.
rec(N,M,MI,R):-N1 is N-MI,M1 is M-MI,fact(N,X),fact(N1,Y),fact(MI,Z),rec(N1,M1,0,R1),R is R1*X^2/(Y^2*Z).
*/
loop3(N,M,I,R1,R):-(M<I;N<I),!,N1 is N^2,ncr(N1,M,X),R is X-R1.
loop3(N,M,I,R1,R):-
rec(N,M,I,R2),R3 is R1+R2,I1 is I+1,loop3(N,M,I1,R3,R).
rec(N,M,MI,R):-N1 is N-MI,M1 is M-MI,fact(N,X),fact(N1,Y),fact(MI,Z),loop3(N1,M1,1,0,R1),R is R1*X^2/(Y^2*Z).
loop2(N,M,MI,R,R):-(M<MI;N<MI),!.
loop2(N0,M,MI,R1,R):-
rec(N0,M,MI,R2), R3 is R1+R2*MI,MI1 is MI+1,loop2(N0,M,MI1,R3,R).
loop1(_,M0,M,R,R):-M=:=M0+1.
loop1(N0,M0,M,R1,R):-
loop2(N0,M,1,0,R2),N1 is N0^2,ncr(N1,M,X),
R3 is R1+R2/X*1000,M1 is M+1,loop1(N0,M0,M1,R3,R).
solve(N0,M0,R):-loop1(N0,M0,1,0,R1),R is floor(R1), !.
start:-f(N,K,A),(solve(N,K,R)->(R==A->Str=" ok ";Str=" no ");
write(fail)),write(" "),write(Str),writeln(R),fail.
start.
f(2, 2, 1666).
f(2, 3, 1666).
f(3, 3, 2642).
f(3, 4, 2928).
f(4, 3, 3228).
f(4, 4, 3967).
f(4, 5, 4428).
f(4, 7, 4797).
f(4, 16,4857).
大した躓きもなく解答出来て悦に入っていましたが、下書きを始めたころに出題者の
「うまく書くと H(10^5, 10^5) ぐらいでも余裕で１秒切れるよ」
というtweetを見つけてしまいました。
->https://x.com/riverplus/status/866625907020095490
recという関数を多用していますので、ここをうまく纏められればよいのですが、
ネストが深すぎて、たいして強くもない私の脳力ではこんがらかって無理です。
そもそもbinomial(ncr)をまともに計算している時点で大きな数は扱えません。
F(n,m)単独よりG(n,m)のほうが簡単になるのかとも思いましたが、
うまくいきません。
と敗北宣言をしていたのですが、しつこく挑戦していたら、投稿前日に解けました。
解けてみると結構正解の周りをさまよっていた気がします。
f(n,m,k)=f(n-1,m-1,k-1)*n^2/k,
f(n-1,m-1,0)=binomial(n-1)^2,m-1)-Σ(1<=k<=m-1)f(n-1,m-1,k)なので
F(n,m)*binomial(n^2,m)=Σ(1<=k<>m)f(n,m,k)*k
Σ(1<=k<=m)f(n,m,k)*k (2024.10.6 <>を訂正)
=n^2*(f(n-1,m-1,0)+Σ(1<=k<=m-1)f(n-1,m-1,k))
=n^2*binomial(n-1)^2,m-1)
F(n,m)=n^2*binomial(n-1)^2,m-1)/binomial(n^2,m)
このままでは大きい数は扱えませんので、
F(n,m)=F(n,m-1)*vという形に変換します。
#ver 1.53, 1.73,1.10.3
function solve(n,m,a)
sum = tmp = 1
for i in 2 : m
# sum+=n^2*binomial((n-1)^2,i-1)/binomial(n^2,i)
tmp *= i * ((n-1)^2-i+2) / ((i-1) * (n^2-i+1))
sum += tmp
end
ans=Int(floor(sum * 1000))
s = a == ans ? "ok " : "no "
println(s,ans)
end
solve(2, 2, 1666)
solve(2, 3, 1666)
solve(3, 3, 2642)
solve(3, 4, 2928)
solve(4, 3, 3228)
solve(4, 4, 3967)
solve(4, 5, 4428)
solve(4, 7, 4797)
solve(4, 16,4857)
#solve(10^5,10^5,Inf)