🤖

同じものを含む円順列の総数(julia版)

2021/08/11に公開約3,100字

はじめに

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_mm種類の球がそれぞれ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(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

大丈夫そうです!

Discussion

ログインするとコメントできます