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

本記事は[素数大富豪 Advent Calendar 2018](https://adventar.org/calendars/2963)の5日目の記事かつ[Wolfram Mathematica Advent Calendar 2018](https://qiita.com/advent-calendar/2018/mathematica)の5日目の記事です。

「[素数大富豪]()」は、数学イベント「[MATH POWER](http://mathpower.sugakubunka.com/)」で大会が行われるなど、巷で人気のトランプゲームです。

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

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


## ユーティリティー集

### 素数判定（定義）

```
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偶数の抽出](#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偶数の抽出](#qk偶数の抽出)
[使い方3: 抽出・除外](#抽出除外)

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

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

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

[使い方1: 3の倍数の抽出](#3の倍数の抽出)
[使い方2: 3で割り切れない数の抽出](#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](https://gist.github.com/butchi/160eb11a5a81c2d42bfa752c06026c0c)にも上げましたのでご活用ください。


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