素数大富豪研究のためのMathematicaユーティリティー関数

本記事は素数大富豪 Advent Calendar 2018の5日目の記事かつWolfram Mathematica Advent Calendar 2018の5日目の記事です。

「素数大富豪」は、数学イベント「MATH POWER」で大会が行われるなど、巷で人気のトランプゲームです。

私は第1回MATH POWER杯でベスト4に入る実績を持っているのですが、最近はもりしーさんをはじめとした北海道勢が強すぎて足元にも及ばないので、自分の持ち味のプログラミングを武器とした研究で勝たねばと思っているところです。

以下、関数定義と使い方をざーっと書き下します。

ユーティリティー集

素数判定（定義）

superscript[a_, b_] := If[b == 1, a, Superscript[a, b]]
centerDot[li_] := If[Length[li] == 1, First[li], CenterDot @@ li]
j[n_] := If[PrimeQ[n], ToString[PrimePi[n]] <> "番目の素数", 
  centerDot[superscript @@@ FactorInteger[n]]]

使い方

倍数（定義）

引数が check の倍数かどうかを判定する純関数です。
3 の倍数かつ 5 の倍数か、など、複数の倍数をチェックすることもできます。

今回、関数にオプションを付けられる OptionsPattern と OptionValue を初めて使ってみましたが、正直、この後も出てくる origin オプションのデフォルト値は True か False のどちらにする方がいいか揺れているので、今後改変する可能性ありです。

divisible[check_Integer, OptionsPattern[origin -> True]] := 
 Divisible[#, check] && If[OptionValue[origin], True, # != check] &
divisible[checkLi_List, OptionsPattern[origin -> True]] := 
 AnyTrue[Table[divisible[check, origin -> OptionValue[origin]][#], {check, checkLi}], TrueQ] &

使い方

QK偶数・QK奇数（定義）

下一桁だけを見て合成数かどうかを判定できるのが2の倍数と5の倍数（下一桁が 0, 2, 4, 5, 6, 8）の場合であることから、素数大富豪においては「5は偶数」と言われています。

言い方が難しいので、今後2の倍数または5の倍数のことを「QK偶数」、そうでない数を「QK奇数」と定義することにします。

qkEvenQ := divisible[{2, 5}]
qkOddQ := ! qkEvenQ[#] && # != 1 &

使い方1: QK偶数の抽出
使い方2: 抽出・除外

倍数チェック（定義）

1001チェックを通過するかどうかなどを調べる関数を定義します。
この関数が一番研究に役立つと思います。

チェックに使われる主な数は
$91 (= 7 \times 13$ )
$969 (= 3 \times 17 \times 19$ )
$1001 (= 3 \times 17 \times 19$ )
$2001 (= 3 \times 23 \times 29$ )
です。

checkQ[checkLi_, OptionsPattern[origin -> False]] := ! AnyTrue[Table[
  divisible[check, origin -> OptionValue[origin]][#],
  {check, checkLi}], TrueQ] &

pass91checkQ := checkQ[{7, 13}]
pass969checkQ := checkQ[{3, 17, 19}]
pass1001checkQ := checkQ[{7, 11, 13}]
pass2001checkQ := checkQ[{3, 23, 29}]

使い方1: 倍数チェック
使い方2: 倍数チェック（応用）

select, reject（定義）

Select 関数を拡張して複数の条件の指定を可能とする select 関数を定義し、また条件に当てはまるものを除外する reject 関数も定義します。

selectAnyTrue と rejectAnyTrue は条件のどれかに当てはまるものだけ抽出／除外する関数です。

select[li_, crit_] := Select[li, crit]
selectAllTrue[li_, critLi_List] := select[li, AllTrue[Table[crit[#], {crit, critLi}], TrueQ] &];
selectAnyTrue[li_, critLi_List] := select[li, AnyTrue[Table[crit[#], {crit, critLi}], TrueQ] &];
select[li_, critLi_List] := selectAllTrue[li, critLi]

reject[li_, crit_] := select[li, ! crit[#] &]
rejectAllTrue[li_, critLi_List] := reject[li, AllTrue[Table[crit[#], {crit, critLi}], TrueQ] &];
rejectAnyTrue[li_, critLi_List] := reject[li, AnyTrue[Table[crit[#], {crit, critLi}], TrueQ] &];
reject[li, critLi_List] := rejectAllTrue[li, critLi]

素数の抽出・除外（定義）

selectPrime[li_] := select[li, PrimeQ]
rejectPrime[li_] := reject[li, PrimeQ]

使い方1: 素数の抽出
使い方2: 合成数の抽出

偶数の抽出・除外（定義）

selectEven[li_] := select[li, EvenQ]
rejectEven[li_] := reject[li, EvenQ]

使い方1: 偶数の抽出
使い方2: 奇数の抽出

倍数の抽出・除外（定義）

selectDivisible[li_, checkLi_List, OptionsPattern[origin -> False]] := select[li, divisible[checkLi, origin -> OptionValue[origin]]]
selectDivisible[li_, check_Integer, OptionsPattern[origin -> False]] := selectDivisible[li, {check}, origin -> OptionValue[origin]]

rejectDivisible[li_, checkLi_List, OptionsPattern[origin -> False]] := reject[li, divisible[checkLi, origin -> OptionValue[origin]]]
rejectDivisible[li_, check_Integer, OptionsPattern[origin -> False]] := rejectDivisible[li, {check}, origin -> OptionValue[origin]]

使い方1: 倍数の抽出
使い方2: QK偶数の抽出
使い方3: 抽出・除外

3の倍数の抽出・除外（定義）

3の倍数抽出はよく使うので専用の関数を用意しました。

selectM3[li_] := selectDivisible[li, 3]
rejectM3[li_] := rejectDivisible[li, 3]

使い方1: 3の倍数の抽出
使い方2: 3で割り切れない数の抽出

確率（定義）

最後に、抽出された数が全体の数の何%かどうかを表示する関数です。
デフォルトでは割合をそのまま値として出力しますが、
disp -> "account" で 分子 / 分母 の表記、
disp -> "numerical" で値を数値化（デフォルトで 1/5 のところを 0.2 と表示）、
disp -> "percent" でパーセント表記
とオプションを選べます。

percent[x_] := ToString[N[100 x]] <> "%"
ratio[li_, crit_, OptionsPattern[displayStyle -> "none"]] := With[
  {numerator = Length[select[li, crit]], denominator = Length[li]},
  With[
    {r = numerator/denominator, disp = OptionValue[displayStyle]}, 
    Piecewise[{
      {HoldForm[numerator/denominator], disp == "account"},
      {N[r], disp == "numerical"},
      {percent[N[r]], disp == "percent"},
      {r, True}
    }]
  ]
]

使い方

使い方

素数判定

In[1]:= 1001 // j
Out[1]= 7・11・13

In[2]:= j /@ Range[20]
Out[2]= {1, "1番目の素数", "2番目の素数", 2², "3番目の素数",
  2・3, "4番目の素数", 2³, 3², 2・5,
  "5番目の素数", 2²・3, "6番目の素数", 2・7, 3・5,
  2⁴, "7番目の素数", 2・3², "8番目の素数", 2²・5
}

素数の抽出

1から20までの素数

In[3]:= selectPrime[Range[20]]
Out[3]= {2, 3, 5, 7, 11, 13, 17, 19}

合成数の抽出

1から20までの合成数
※ 1を含むので注意

In[4]:= rejectPrime[Range[20]]
Out[4]= {1, 4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20}

倍数

17は17の倍数か

In[5]:= divisible[17][17]
Out[5]= True

17は17を含まない17の倍数か

In[6]:= divisible[17, origin -> False][17]
Out[6]= False

300は3の倍数かつ5の倍数か

In[7]:= divisible[{3, 5}][300]
Out[7]= True

偶数の抽出

In[8]:= selectEven[Range[20]]
Out[8]= {2, 4, 6, 8, 10, 12, 14, 16, 18, 20}

奇数の抽出

In[9]:= rejectEven[Range[20]]
Out[9]= {1, 3, 5, 7, 9, 11, 13, 15, 17, 19}

倍数の抽出

1から100までのうち13の倍数（13そのものは含まない）
※ 13そのものを含むときは origin->True オプションを与える

In[10]:= selectDivisible[Range[100], 13]
Out[10]= {26, 39, 52, 65, 78, 91}

QK偶数の抽出

In[11]:= selectDivisible[Range[20], {2, 5}]
Out[11]= {4, 6, 8, 10, 12, 14, 15, 16, 18, 20}
In[12]:= select[Range[20], qkEvenQ]
Out[12]= {2, 4, 5, 6, 8, 10, 12, 14, 15, 16, 18, 20}

※ 正式なQK偶数は後者

抽出・除外

QK奇数の抽出

In[13]:= select[Range[50], qkOddQ]
Out[13]= {3, 7, 9, 11, 13, 17, 19, 21, 23, 27, 29, 31, 33, 37, 39, 41, 43, 47, 49}

10で割り切れない数（10を含む）

In[14]:= rejectDivisible[Range[30], 10]
Out[14]= {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 21, 22, 23, 24, 25, 26, 27, 28, 29}

2の倍数でも5の倍数でもない数（2と5を含む）

In[15]:= rejectDivisible[Range[50], {2, 5}]
Out[15]= {1, 2, 3, 5, 7, 9, 11, 13, 17, 19, 21, 23, 27, 29, 31, 33, 37, 39, 41, 43, 47, 49}

3の倍数の抽出

In[16]:= selectM3[Range[20]]
Out[16]= {6, 9, 12, 15, 18}

3で割り切れない数の抽出

In[17]:= rejectM3[Range[20]]
Out[17]= {1, 2, 3, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20}

倍数チェック

91チェックを通過する数

In[18]:= select[Range[90, 110], pass91checkQ]
Out[18]= {90, 92, 93, 94, 95, 96, 97, 99, 100, 101, 102, 103, 106, 107, 108, 109, 110}

969チェックを通過する数

In[19]:= select[Range[950, 970], pass969checkQ]
Out[19]= {953, 955, 956, 958, 959, 961, 962, 964, 965, 967, 968, 970}

1001チェックを通過する数

In[20]:= select[Range[1000, 1010], pass1001checkQ]
Out[20]= {1000, 1002, 1003, 1004, 1005, 1006, 1007, 1009, 1010}

2001チェックを通過する数

In[21]:= select[Range[2000, 2010], pass2001checkQ]
Out[21]= {2000, 2002, 2003, 2005, 2006, 2008, 2009}

1001チェックを通過しない数

In[22]:= reject[Range[1000, 1100], pass1001checkQ]
Out[22]= {1001, 1008, 1012, 1014, 1015, 1022, 1023, 1027, 1029,
1034, 1036, 1040, 1043, 1045, 1050, 1053, 1056, 1057, 1064, 1066,
1067, 1071, 1078, 1079, 1085, 1089, 1092, 1099, 1100}

倍数チェック（応用）

3の倍数となる3桁以下のQK奇数

In[23]:= select[Range[999], {qkOddQ, divisible[3]}]
Out[23]= {3, 9, 21, 27, 33, 39, 51, 57, 63, 69, 81, 87, 93, 99, 111,
117, 123, 129, 141, 147, 153, 159, 171, 177, 183, 189, 201, 207, 213,
219, 231, 237, 243, 249, 261, 267, 273, 279, 291, 297, 303, 309, 321,
327, 333, 339, 351, 357, 363, 369, 381, 387, 393, 399, 411, 417, 423,
429, 441, 447, 453, 459, 471, 477, 483, 489, 501, 507, 513, 519, 531,
537, 543, 549, 561, 567, 573, 579, 591, 597, 603, 609, 621, 627, 633,
639, 651, 657, 663, 669, 681, 687, 693, 699, 711, 717, 723, 729, 741,
747, 753, 759, 771, 777, 783, 789, 801, 807, 813, 819, 831, 837, 843,
849, 861, 867, 873, 879, 891, 897, 903, 909, 921, 927, 933, 939, 951,
957, 963, 969, 981, 987, 993, 999}

確率

4桁のQK奇数が素数である確率

In[24]:= ratio[Range[1000, 9999], PrimeQ, displayStyle -> "percent"]
Out[24]= "11.7889%"

1001チェックと2001チェックを通過する4桁のQK奇数が素数である割合

In[25]:= ratio[
  select[Range[1000, 9999],{qkOddQ, pass1001checkQ, pass2001checkQ}],
  PrimeQ,
  displayStyle -> "account"
]

Out[25]= 1061/1594

終わりに

今回の記事は基本的にMathematica使いじゃないと有用ではありませんが、このユーティリティー関数をベースにプログラムを組めば、素数大富豪大会で勝ち抜く戦略を効率的に研究できるようになると思います。
Mathematicaならではの関数定義もいくらかありますが、Pythonなど他の言語でも同じようなユーティリティー関数が出てくれば素数大富豪研究がもっと広まるんじゃないかと思います。

ユーティリティー関数はGistにも上げましたのでご活用ください。

明日の素数大富豪 Advent Calendarはtatyam_primeさんの「素数判定アプリを作ってみたお話」です。一体どんなプログラムが飛び交うのでしょうか！