同じものを含む円順列の総数（julia版）

はじめに

8/10にtwitterで次のような数学の問題が流れてきました。

https://twitter.com/TubeSoling/status/1424748510923018251

場合の数の問題ですが，リツートでc++で解いているものがありました。

https://twitter.com/OJun46225932/status/1424967384570335235

面白いなあ。。。と思って，juliaでもできそうだなと思って，やってみました。

20210810.jl

using Combinatorics
seq = [1,1,1,4,4,5]
p=union(permutations(seq))
p[[8]]

1-element Vector{Vector{Int64}}:
[1, 1, 4, 5, 1, 4]

ものすごいシンプルで数学っぽいです！

ちなみにpの中身はこんな感じです。

p

60-element Vector{Vector{Int64}}:
[1, 1, 1, 4, 4, 5]
[1, 1, 1, 4, 5, 4]
[1, 1, 1, 5, 4, 4]
[1, 1, 4, 1, 4, 5]
[1, 1, 4, 1, 5, 4]
[1, 1, 4, 4, 1, 5]
[1, 1, 4, 4, 5, 1]
[1, 1, 4, 5, 1, 4]
[1, 1, 4, 5, 4, 1]
[1, 1, 5, 1, 4, 4]
[1, 1, 5, 4, 1, 4]
[1, 1, 5, 4, 4, 1]
[1, 4, 1, 1, 4, 5]
⋮
[4, 5, 1, 4, 1, 1]
[4, 5, 4, 1, 1, 1]
[5, 1, 1, 1, 4, 4]
[5, 1, 1, 4, 1, 4]
[5, 1, 1, 4, 4, 1]
[5, 1, 4, 1, 1, 4]
[5, 1, 4, 1, 4, 1]
[5, 1, 4, 4, 1, 1]
[5, 4, 1, 1, 1, 4]
[5, 4, 1, 1, 4, 1]
[5, 4, 1, 4, 1, 1]
[5, 4, 4, 1, 1, 1]

以前，同じものを含む円順列の総数の求め方をまとめていたので，「これはjuliaでできるのではないか？」と思ってチャレンジです。

同じものを含む円順列の総数

以前，同じものを含む円順列の総数はMathematicaで関数化しました。

https://qiita.com/dannchu/items/d004e96e8f5114fe59de

一般化

友田 勝久さんがまとめているものを参考にしました。
https://tomodak.com/report/junretsu2012.pdf

$S_1,\,S_2,\,\cdots,\,S_m$の$m$種類の球がそれぞれ$n_1,\,n_2,\,\cdots,\, n_m$個の球が作る円順列の総数を $f(n_1,\,\cdots,\,n_m)$と表すと，

$$f(n_1,\,\cdots,\,n_m)=\frac1N\sum_{k|l}\phi(k)\frac{\left(\frac{N}k\right)!}{\left(\frac{n_1}k\right)!\left(\frac{n_2}k\right)!\cdots\left(\frac{n_m}k\right)!}$$

• $l=\gcd(n_1,\,\cdots,\,n_m)$，$N=n_1+n_2+\cdots+n_m$です。

• $\sum_{k|l}$は$l$の正の約数である$k$で和をとります。

• $\phi(k)$はオイラー(Euler)のトーシェント関数で， $k$以下で$k$と互いに素な自然数の個数を返します。

コード化（その1）

最初，3つのパッケージProjectEulerUtil.jl Combinatorics.jl Primes.jlを使うことにしました。数学のパッケージ揃ってます！

enkan1.jl

using ProjectEulerUtil  #proper_divisors(l) 約数のリストを求める。
using Combinatorics #multinomial(a)多項係数を求める
using Primes    #totient(i)　オイラーのトーシェント関数

function enkan1(a)
    l=gcd(a)    #リスト（配列）aの最大公約数を求める。
    N=sum(a)    #aの総和
    A=union(push!(proper_divisors(l),l))    #lの約数のリスト。ちょっと使いにくい。
        #lの約数でl自身は除かれてしまうので，push!でlをリストに加えている。
        #しかし，l=1の時は除かれない。l=1の時は[1,1]となってしまうので，unionで重複を除いている。
    p=0
    for k in A
        q=map(x -> x÷k,a) 
        p +=totient(k)*multinomial(q...)
    end
    println(p÷N)
end

コード化（その2）

proper_divisorsがちょっと使いにくかったので，約数のリスト求めるdivisorsという関数を作ることにしました。（どこかのパッケージにあるような気もします。。。）

enkan.jl

using Combinatorics #multinomial(a)多項係数を求める
using Primes    #totient(i)　オイラーのトーシェント関数

function divisors(n)    #約数のリストを求める関数
    X=[]
    for i=1:n
        if n % i==0
            X =push!(X,i)
        end
    end
    X
end

function enkan(a)
    l=gcd(a)    #リスト（配列）aの最大公約数を求める。
    N=sum(a)    #aの総和
    A=divisors(l)
    p=0
    for k in A
        q=map(x -> x÷k,a) 
        p +=totient(k)*multinomial(q...)
    end
    p÷N
end

これで完成です！いろいろ確かめてみます。

赤球4個，白球6個を円形に並べる総数

enkan([4 6])

22

赤球2個，青球3個，白球4個を円形に並べる総数

enkan([2 3 4])

140

赤球3個，青球3個，白球3個を円形に並べる総数

enkan([3 3 3])

188

赤球4個，青球4個，白球4個，黒球4個を円形に並べる総数

enkan([4 4 4 4])

3941598

大丈夫そうです！