🌻

Juliaで同じものを含む円順列・数珠順列の総数を求める

2025/09/07に公開

6件

 はじめに9月は学園祭・文化祭のシーズンです。Xのポストで東大寺学園の文化祭での数学の問題が流れてきました。

https://x.com/tdjsuken2/status/1964116342132572239
その中で，数珠順列の問題があったので，もう一度考えてみたくなりました。

 以前までの到達点円順列に関しては，コードでリストを書いて数えなくても，考え方で数えられることはわかっていました。以前の私のブログでもまとめています。
https://zenn.dev/dannchu/articles/2a27581585a338
その後，数珠順列についても考えてみたのですが，うまくまとまらず，リストを作成して数えていました。
https://zenn.dev/dannchu/articles/4d35b5d2b4c94c
一方で，数珠順列の総数について求めるアルゴリズムについて書いてあるサイトもありました。
https://www.chart.co.jp/subject/sugaku/suken_tsushin/68/68-1.pdf
この記事を書いた山田一男さんによると，重複数珠順列については，場合分けが4つほどあり，そこからはまとめるには至っていませんでした。今回は，この場合分けをJulia言語で実装し，できればシンプルにまとめられないかを考えました。

 4つの場合分けについてこの4つの場合分けは，同じものを含む円順列を作った後に，裏に返して，同じ並びがあれば，1通りと数えるための場合分けです。簡単のために色の球で数珠を作ることを考えます。
＜場合分け＞
各色の球の個数が全て偶数個　例えば🔴🔴🔴🔴🔵🔵🔵🔵⚪️⚪️
奇数個が1色，あとは全て偶数個　例えば🔴🔴🔴🔵🔵🔵🔵⚪️⚪️
奇数個が2色，あとは全て偶数個　例えば🔴🔴🔴🔵🔵🔵⚪️⚪️
上記以外　例えば🔴🔴🔴🔵🔵🔵⚪️⚪️⚪️

 今回の結果結論から言うと，この4つの場合を次のようにまとめました。結構シンプルになったので，満足しています。
!

 Juliaのコードusing Combinatorics  # multinomial
using Primes         # totient, divisors

# 同じものを含む円順列の総数
function enkan(a)
    l = gcd(a)        # a の最大公約数
    N = sum(a)             # 総和
    p = 0
    for k in divisors(l)
        q = div.(a, k)     # 成分ごと整数除算
        p += totient(k) * multinomial(q...)
    end
    return p ÷ N
end

# 同じものを含む数珠順列の総数
function juzu(a)
    N = sum(a)
    t = enkan(a)
    q = div.(a, 2)
    m = count(isodd, a)

    # 反転の寄与
    if m ≤ 2
        t += multinomial(q...)
    end

    return t ÷ 2
end

 いくつかの例@show juzu([3,2,2])
@show juzu([3,3,3])
@show juzu([3,2,1])
@show juzu([4,6])
@show juzu([4,2,1])
@show juzu([6,3])
@show juzu([2,2,2]);

juzu([3, 2, 2]) = 18

juzu([3, 3, 3]) = 94

juzu([3, 2, 1]) = 6

juzu([4, 6]) = 16

juzu([4, 2, 1]) = 9

juzu([6, 3]) = 7

juzu([2, 2, 2]) = 11
これで，数珠順列も怖くないですね！

Discussion

mixedmoss

はじめまして．私も以前全く別の方法で「同じものを含む数珠順列」の公式を作って簡単な証明をつけたのですが，合っているかどうか不安だったのですが，今回の公式が私の得た公式と同じだったので安心しました．記事をアップしていただきありがとうございます．
なお，私のやり方は

に置いてあります．結果は一致したものの（いつもの様に）やり方がかなり異なるので，結果に関しては自信が持てましたが，証明に関してはいまだ確信がありません．

清水団

ありがとうございます。@mixedmoss さんのサイト見てみますね。
Mathematicaのノートブックを見ましたが，まず，円順列は全リストを数え上げているのでしょうか？
そのリストを元に数珠を作るときの対称性で場合分けしているように見えます。

もし，この方法だとすれば，数が多くなってきた時に処理が大変です。（時間がかかります）
私の式のまとめは，リストアップはしません。計算式で円順列や数珠順列の総数だけを求めています。
（リストアップもしたことがあります。）

数珠の場合分けのスタートは同じですが，私はそれらがまとめられると考えました。（それをこのサイトでまとめてみました）着目したのは「奇数個の種類数」です。

Mathematicaに変換したコードも載せておきますね。

(* 同じものを含む円順列の総数 *)
enkan[a_List] := Module[
  {l, N, p = 0, q},
  l = Apply[GCD, a];          (* a の最大公約数 *)
  N = Total[a];              (* 総和 *)

  Do[
    q = Quotient[a, k];       (* 成分ごと整数除算 *)
    p += EulerPhi[k] * Apply[Multinomial, q],
    {k, Divisors[l]}
  ];

  Quotient[p, N]              (* p ÷ N *)
];

(* 同じものを含む数珠順列の総数 *)
juzu[a_List] := Module[
  {N, t, q, m},
  N = Total[a];
  t = enkan[a];
  q = Quotient[a, 2];
  m = Count[a, _Odd];         (* 奇数成分の個数 *)

  (* 反転の寄与 *)
  If[m <= 2,
    t += Apply[Multinomial, q]
  ];

  Quotient[t, 2]              (* t ÷ 2 *)
];

enkan[{2,2,2}]

juzu[{2,2,2}]

mixedmoss

私はある軸に関して対称な円順列に注目しました．（円順列全てを数え上げてはいません.) そしてそれを最初はリストで数え上げました．すると規則性に気付いたので，公式化して簡単に証明し，それを使って場合の数のみを計算するプログラムも作りました．プログラムは，単純なリスト，その図形表示, 発見した公式による計算の大きく分けて3つを作りました．ついでに Mathematica player でも実行できる物も作成しました. 清水団様のプログラムを拝見しましたが，これは「公式による計算」に当たるのでしょうか？

清水団

ありがとうございます。円順列，数珠順列ともに「リストを列挙して数えている」のではなく，「公式による計算」となります。
「私はある軸に関して対称な円順列に注目しました．」とありますが，できれば一度コードを見てみたいです。

（紹介していただいたサイトに載っているかな？）
「リストを列挙して数えている」とは円順列であれば，1列に並べる順列を作り，それを円形にして，（shiftなどで順繰りにずらして）同じ並びをキャンセルするようなコードをはじめは作りました。数珠も同様で，裏返しにして，同じ並びをキャンセルしてました。このコードであれば，並びものリストも得られるのですが，球が多くなると，求めるのに時間がかかることがわかっていました。

mixedmoss

返信ありがとうございます。「リストを列挙して数える方法」は，清水団様のやり方と基本的に

同じだとおもいます．(shiftとreverseを使ってリストを変換し重なるものを消去) 従って数が大きいと確かに実行不可能です．
コードは私のページの [Notebook download]　に，証明については[PDF]に載っています．また [Notebook mobile download]　の方は Manuplate を使用して，Mathematicaをお持ちでない方でも実験できるようにしたファイルです．特に「証明」についてご意見を頂けるとありがたいです．（自分で書いた証明はいつも自信が持てないです.大学入試レベルなら自信満々ですが．．．)
私見ですが，Mathematica を使ってプログラムを書いたときに残念なのは，GeoGebraと異なり，そのプログラムをweb上で走らせられないことです．確かに Wolfram player と Manipulate を使えばある程度自由に Mathematicaを持っていない人でも試すことができますが，wolfram player をダウンロードしてまで見てくださる方は非常に少ないと思います．またwolfram cloud も低速で使い物になりません この点は是非ともwolframに頑張ってほしいですね．

清水団

数が大きくなっても大丈夫なように，コードを少し変更しました。

using Combinatorics  # multinomial
using Primes         # totient, divisors

# 同じものを含む円順列の総数
function enkan(a)
    l = gcd(a)        # a の最大公約数
    N = sum(a)             # 総和
    p = big(0)
    for k in divisors(l)
        q = div.(a, k)     # 成分ごと整数除算
        p += totient(k) * multinomial(big.(q)...)
    end
    return p ÷ N
end

# 同じものを含む数珠順列の総数
function juzu(a)
    N = sum(a)
    t = enkan(a)
    q = div.(a, 2)
    m = count(isodd, a)

    # 反転の寄与
    if m ≤ 2
        t += multinomial(big.(q)...)
    end

    return t ÷ 2
end

# 大きな数値を見やすく表示する関数（適当に改行）
function show_result(n, width=60)
    s = string(n)
    for i in 1:width:length(s)
        println(s[i:min(i+width-1, length(s))])
    end
    return n
end

show_result(juzu([42,57,334,346,765,1729]))

165701354457414537858061853742887351161347286174442275924096
442670012868851272543023772725632155825457280023828709327839
176007020595333782588752768900553593209656916173846649626888
168220439948044573987243541151148795568178652118574717837178
281014816078207634398116903983784887323696238973507378072681
179569490050432774383776641817280954911274024092391993360581
638876752375318842782195340413918711705807350760099495005967
319969936656248564344002442915674126033064071021590829879975
364650597397262100893457774592909391434902424130808272958871
981091370461498096081910033242938394534991700126756977103936
530149678868693448464964129695075622209140382817030292049119
898419112817621745407108176544224263199124579471294039646334
158189271640698499204864570885015111568015855366729874710632
179541285684185139076914312532562671048615763109813469618726
224599584460279560688063780659816715111274295609198100641267
087595892292479342795473354331182735679958142357576398459803
868629588395096382945986117412418079382536273244760156148310
858954779686717272785061091939250166891086298754609225879162
981357278811969757382151065276470287142319997026494854122274
662912878747174468252283600333372608855308397707586501248930
566165586633382271306244658884218149963775149165861300873378
238337904773360924148147576937669747062093270740050834579686
452537014346576772590601758267852200656076867790393043781435
294946692962930401264791329191928005975195794730058033551190
782586414232538570754733257730238238277750930238887217808445
406322922599009781349266086144569549944545718597957336637221
184469103615713188950755670869923099146495288990759512158379
190722308725620193751762063708253681383668867766736801786830
728950860921777357931600426125078512310575555753816194523220
377886512137882187201441805425394534360302092861440000000000