CodeIQ:「ロンリー・ルーク」問題の解答

Kawazoe(@riverplus)さんのCodeIqの問題です。
問題とテストケースはこちらにお世話になりました。解答もあります。
~~->https://blog.goo.ne.jp/r-de-r/e/f2b057442e37b97f008680d1eb1c087c~~
(リンク切れになっていました2025.5.2)
URLが変わっていましたが見つかりました。(2025.5.3追加)
->https://sangaku0418.hatenablog.com/entry/2017/05/18/100100

ーーー2025.5.2追加
n,kを自然数とし、一辺nマスの正方形の盤に重ならないようにk個の駒を置くとき、
自身から見て縦横に駒がないときはぐれルークと言い、はぐれルークの数の期待値を
F(n,k)とし、1<=k<=mのすべてのkに対するF(n,K)の和をG(n,m)とするとき
10^3とG(n,m)の積の整数部分を求める問題です。
但し　1≦n≦4, 1≦m≦n^2　です。
ーーー

こちらにも情報があります。
->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しか表示されず実行できません）

```Julia
#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しか表示されず実行できません）

```prolog
%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という形に変換します。

```Julia
#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)
```