数論・素数に関するメソッドを実装する

数論・素数に関するメソッドを C# で一通り実装してまとめました。
通常の数値計算プログラミング用にも、また競技プログラミング用のライブラリとしても使えます。

機能は次の通りです。

Euclid 互除法

  • 最大公約数 (GCD)
  • 最小公倍数 (LCM)
  • ax + by = 1 の解 (Euclid 互除法の拡張)

素数

下の O(x) は計算量のオーダーを表します。

  • 素因数分解
    • O(√n)
  • 約数の列挙
    • O(√n)
  • 素数判定
    • O(√n)
  • n 以下の素数の列挙
    • およそ O(n)
  • m 以上 M 以下の素数の列挙
    • およそ O(√M + (M – m))

ソースコードは以下の通りです。
ここではシンプルな実装とし、例外処理も全体的に省略しています。
さらに無駄な演算処理を除くように改良することで、もう少し速くすることはできるはずです。

using System;
namespace AlgorithmLib.Numerics
{
public static class Euclidean
{
public static int Gcd(int a, int b) { for (int r; (r = a % b) > 0; a = b, b = r) ; return b; }
public static int Lcm(int a, int b) => a / Gcd(a, b) * b;
public static long Gcd(long a, long b) { for (long r; (r = a % b) > 0; a = b, b = r) ; return b; }
public static long Lcm(long a, long b) => a / Gcd(a, b) * b;
// ax + by = 1 の解
// 前提: a と b は互いに素。
// ax + by = GCD(a, b) の解を求める場合、予め GCD(a, b) で割ってからこの関数を利用します。
public static (long x, long y) ExtendedEuclid(long a, long b)
{
if (b == 0) throw new ArgumentOutOfRangeException(nameof(b));
if (b == 1) return (1, 1 - a);
var q = Math.DivRem(a, b, out var r);
var (x, y) = ExtendedEuclid(b, r);
return (y, x - q * y);
}
}
}
view raw Euclidean.cs hosted with ❤ by GitHub
using System;
using System.Collections.Generic;
namespace AlgorithmLib.Numerics
{
public static class Primes
{
/// <summary>
/// 指定された自然数を素因数分解します。O(√n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>素因数のコレクション。n = 1 の場合は空の配列。</returns>
public static long[] Factorize(long n)
{
var r = new List<long>();
for (long x = 2; x * x <= n && n > 1; ++x)
while (n % x == 0)
{
r.Add(x);
n /= x;
}
// √n を超える素因数はたかだか 1 個であり、その次数は 1。
if (n > 1) r.Add(n);
return r.ToArray();
}
/// <summary>
/// 指定された自然数の約数をすべて求めます。O(√n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>約数のコレクション。</returns>
public static long[] Divisors(long n)
{
var r = new List<long>();
for (long x = 1; x * x <= n; ++x)
if (n % x == 0) r.Add(x);
for (int i = r.Count - 1; i >= 0; --i)
{
var v = n / r[i];
if (v > r[i]) r.Add(v);
}
return r.ToArray();
}
/// <summary>
/// 指定された自然数が素数かどうかを判定します。O(√n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>素数である場合に限り true。</returns>
public static bool IsPrime(long n)
{
for (long x = 2; x * x <= n; ++x)
if (n % x == 0) return false;
return n > 1;
}
/// <summary>
/// 指定された自然数以下の素数をすべて求めます。およそ O(n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>素数のコレクション。</returns>
public static long[] GetPrimes(long n)
{
var b = new bool[n + 1];
for (long p = 2; p * p <= n; ++p)
if (!b[p])
for (long x = p * p; x <= n; x += p)
b[x] = true;
var r = new List<long>();
for (long x = 2; x <= n; ++x) if (!b[x]) r.Add(x);
return r.ToArray();
}
/// <summary>
/// 指定された範囲内の素数をすべて求めます。およそ O(√M + (M - m))
/// </summary>
/// <param name="m">最小値。</param>
/// <param name="M">最大値。</param>
/// <returns>素数のコレクション。</returns>
public static long[] GetPrimes(long m, long M)
{
var rM = (int)Math.Ceiling(Math.Sqrt(M));
var b1 = new bool[rM + 1];
var b2 = new bool[M - m + 1];
if (m == 1) b2[0] = true;
for (long p = 2; p <= rM; ++p)
if (!b1[p])
{
for (var x = p * p; x <= rM; x += p) b1[x] = true;
for (var x = Math.Max(p, (m + p - 1) / p) * p; x <= M; x += p) b2[x - m] = true;
}
var r = new List<long>();
for (int x = 0; x < b2.Length; ++x) if (!b2[x]) r.Add(m + x);
return r.ToArray();
}
}
}
view raw Primes.cs hosted with ❤ by GitHub

使用例を次に示します。

using System;
using AlgorithmLib.Numerics;
namespace AlgorithmConsole
{
class Program
{
static void Main(string[] args)
{
// 素因数分解
Console.WriteLine(string.Join(" * ", Primes.Factorize(2020)));
// 2 * 2 * 5 * 101
// 約数の列挙
Console.WriteLine(string.Join(" ", Primes.Divisors(2020)));
// 1 2 4 5 10 20 101 202 404 505 1010 2020
// 素数判定
Console.WriteLine(Primes.IsPrime(1000000007));
// True
// n 以下の素数の列挙
Console.WriteLine(string.Join(" ", Primes.GetPrimes(100)));
// 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
// m 以上 M 以下の素数の列挙
Console.WriteLine(string.Join(" ", Primes.GetPrimes(1000000000, 1000000100)));
// 1000000007 1000000009 1000000021 1000000033 1000000087 1000000093 1000000097
}
}
}
view raw Program.cs hosted with ❤ by GitHub

前回: 二分探索のライブラリ化

作成したサンプル

バージョン情報

  • C# 7.0
  • .NET Standard 2.0

参照

カテゴリー: アルゴリズム, 数学. タグ: . 2 Comments »

transform と deviceorientation における回転の表現 (HTML)

CSS の transform プロパティと JavaScript の deviceorientation イベントではともに 3D の回転状態 (姿勢、傾き) が登場しますが、
その扱い方に差があるため検証しました。
deviceorientation は、デバイスのジャイロ センサーが回転状態を通知することで発生するイベントです。

HTML の 3 次元座標系では、2 次元スクリーン座標系の x 軸および y 軸に加えて、スクリーンに垂直な z 軸が存在します。
デバイス (スマートフォンなど) を水平に持ち、北を向いた状態を基準に考えます。

このとき、CSS の transform プロパティと JavaScript の deviceorientation イベントにおける、
回転に関する性質の違いを下の表にまとめました。例えば z 軸を中心とする回転の角度は、現在は北を向いていたとしたら、
transform では東側 (時計回り) を向くと正、deviceorientation では西側を向くと正になります。

  transform deviceorientation
座標系 左手系 右手系
x 軸 右が正 右が正
y 軸 下 (手前) が正 上 (奥) が正
z 軸 (画面が水平のとき) 鉛直の上が正 鉛直の上が正
回転角度 回転軸の正方向に左ねじを進める場合が正 回転軸の正方向に右ねじを進める場合が正

 

検証のため、CSS の transform プロパティと JavaScript の deviceorientation イベントを利用して、
デバイスの回転状態を画面内の立方体オブジェクトに同期させるサンプルを作成しました。

DeviceOrientation

ジャイロ センサーを搭載した端末であれば、こちらのテストページで確認できます。
HTML のソースは以下の通りです。

<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width" />
<title>Device Orientation</title>
<style type="text/css">
.container {
width: 2px;
height: 2px;
margin: 100px auto;
perspective: 300px;
}
#cube {
transform-style: preserve-3d;
}
.face {
width: 100px;
height: 100px;
position: absolute;
left: -50px;
top: -50px;
background: rgba(255, 102, 0, 0.5);
border: 2px solid gray;
text-align: center;
font-size: 60px;
line-height: 100px;
}
.front {
transform: translateZ(50px);
}
.back {
transform: rotateY(180deg) translateZ(50px);
}
.right {
transform: rotateY(90deg) translateZ(50px);
}
.left {
transform: rotateY(-90deg) translateZ(50px);
}
.top {
transform: rotateX(90deg) translateZ(50px);
}
.bottom {
transform: rotateX(-90deg) translateZ(50px);
}
</style>
</head>
<body>
<div id="log"></div>
<div class="container">
<div id="cube">
<div class="face front">1</div>
<div class="face back">6</div>
<div class="face right">2</div>
<div class="face left">5</div>
<div class="face top">3</div>
<div class="face bottom">4</div>
</div>
</div>
<script type="text/javascript">
var logEl = document.querySelector("#log");
var cubeEl = document.querySelector("#cube");
if (window.DeviceOrientationEvent) {
window.addEventListener("deviceorientation", function (e) {
logEl.innerHTML = `alpha: ${e.alpha}<br />beta: ${e.beta}<br />gamma: ${e.gamma}`;
cubeEl.style.transform = `rotateZ(${-e.alpha}deg) rotateX(${-e.beta}deg) rotateY(${e.gamma}deg)`;
});
}
</script>
</body>
</html>
view raw cube-sync.html hosted with ❤ by GitHub

以下は、各技術についての説明です。

transform プロパティ

transform プロパティで rotateX などを利用して回転状態を指定する場合、次に示すように複数の回転を重ね合わせることができます。

CSS:
transform: rotateX(45deg) rotateY(30deg) rotateZ(60deg);

JavaScript:
element.style.transform = "rotateX(45deg) rotateY(30deg) rotateZ(60deg)";

ただし、座標系ごと回転させながら左から順に適用します (オイラー角)。
これは、以前に 3D における回転の表現と相互変換で書いた通り、元の座標系のまま右から順に適用する、と考えても同じです。

以下に rotateX(45deg)rotateY(45deg) を組み合わせた例を載せておきます。

初期状態

rotateX(45deg) (左)        rotateX(45deg) rotateY(45deg) (右)

rotateY(45deg) (左)        rotateY(45deg) rotateX(45deg) (右)

deviceorientation イベント

window.addEventListenerdeviceorientation に対するイベントリスナーを登録します。
デバイスの回転状態が変化すると、イベントリスナーが呼び出されます。

引数の alpha, beta, gamma はそれぞれ z 軸、x 軸、y 軸を中心とした回転の角度を表し、
座標系ごと回転させながらこの順に重ね合わせたものが回転状態を表します。
それぞれの値の範囲は次の通りです。

  • z 軸: 0 ≦ alpha < 360
    • 北を向いたとき、alpha = 0
  • x 軸: -180 ≦ beta < 180
  • y 軸: -90 ≦ gamma < 90

 

回転状態の同期

以上から、デバイスの回転状態を画面内のオブジェクトに同期させるには次のようにします。

cubeEl.style.transform = `rotateZ(${-e.alpha}deg) rotateX(${-e.beta}deg) rotateY(${e.gamma}deg)`;

正負の符号に注意します。
結果として、z 軸および x 軸における回転角度の正負は異なり、y 軸では同じになります。

 

次回は、回転状態をネットワーク経由で同期させます。

次回: ASP.NET SignalR でデバイスの回転状態を同期する

作成したサンプル
参照
transform
deviceorientation

SIR 感染症モデルのシミュレーター

// この投稿は ライフゲーム Advent Calendar 2017 の 19 日目の記事です。

2017年11月11日の稲葉寿先生還暦記念祝賀研究集会「数理人口学・数理疫学・構造化個体群モデル」で講演してきました。
そのときに発表したシミュレーション ツールを紹介します。

数理生物学で「SIR 感染症モデル」という数理モデルがあり、

S(t): 未感染者の人口
I(t): 感染者の人口
R(t): 回復者の人口

としたとき、

S’ = θR – βSI
I’ = βSI – γI
R’ = γI – θR

のような微分方程式系で表されるものを指します。ここで、各定数は

β: 感染率
γ: 回復率
θ: 免疫喪失率

を表します (このように免疫喪失を考慮する場合は SIRS ともいう)。

このモデルをもとに、感染症の伝播を視覚的に表現するセルオートマトンを WPF で作成しました。
各ファイルは GitHub にあります。

EpidemicSimulator.exe を実行し、右下のトグルスイッチを押せばシミュレーションを開始できます。
感染率などのいくつかのパラメーターは実行時にリアルタイムに変更できます。

S と I が隣り合っているとき、一定の割合で感染が発生するようになっています (したがって、数式で表した状況とは厳密には異なります)。
また、I と R は一定の割合でそれぞれ R と S に移動します。

実行前に設定するパラメーター:

  • 高さ
  • SIR の初期人口比

実行中も設定できるパラメーター:

  • 感染率
  • 回復率
  • 免疫喪失率
  • Looping Map: マップの端でループするかどうか
  • ターンの時間間隔

Epidemic Simulator

 

実装方法については技術的に目新しいところはありませんが、特徴としては以下が挙げられます。

(1) シミュレーションの演算は UI スレッドではなく、バックグラウンド スレッドで実行
(2) 各フレームで画像データを生成して、Image コントロールで表示

(1) については、重い処理を UI スレッドで実行するとアプリケーションがフリーズしてしまうため、
各フレームで非同期的にデータのスナップショットを取得しています。
技術的な説明は、以前に

で書いた通りです。

バージョン情報
.NET Framework 4.5
ReactiveProperty 3.6.0
ToggleSwitch 1.1.2

参照
稲葉寿先生還暦記念祝賀研究集会「数理人口学・数理疫学・構造化個体群モデル」

カテゴリー: ツール, 数学. タグ: . Leave a Comment »

3D における回転の表現と相互変換

以前に投稿した WPF で 3D オブジェクトを回転させるではオブジェクトの回転の状態を行列で表していましたが、
3 次元空間における回転を表現する方法は、次のように何通りか考えられます。
なお、原点を中心に回転させるものとし、回転角度は、回転軸を表すベクトルの方向に右ねじを押し込む場合を正とします。

  1. 行列
    • 3 次正方行列。とくに、回転を表すものは直交行列であり、P^{-1} = {}^t P が成り立つ
    • WPF では Matrix3D 構造体
  2. 四元数 (クォータニオン)
  3. オイラー角 (ロール、ピッチ、ヨー)
    • 任意の回転を 3 回の単純な回転の合成により表す。詳細の定義はロール・ピッチ・ヨー ― Kinectで学ぶ数学を参照
    • 座標系ごと回転させる
    • ロール、ピッチ、ヨーをそれぞれどの軸に対応させるか、どの順番で作用させるかで結果が変わってしまうため注意が必要
  4. 特定の 2 点の回転前後の座標
    • とくに、直交する 2 つのベクトルを用いるとよい

3D のプログラミングをしていると、API によってどの表現を利用するかが異なることがあります。
以下では、これらの表現を互いに変換する方法について考えます。

■ 行列 → 2 点の座標

任意のベクトルに対して行列を作用させれば回転後のベクトルが求められます。
WPF では演算子 * が用意されています。

■ 四元数 → 行列

四元数の成分から行列の各成分を計算できます (詳細は四元数による回転行列の表現を参照)。
WPF では Matrix3D.Rotate メソッドとして用意されています。

■ オイラー角 → 行列または四元数

ここでは、ヨーは y 軸、ピッチは x 軸、ロールは z 軸のまわりの回転を表し、この順に適用されるとします。
(元の座標系における) ヨー、ピッチ、ロールを表す回転行列をそれぞれ R_y, R_p, R_r とし、それぞれの回転角度を \theta_y, \theta_p, \theta_r とします。
すなわち、

R_y = \left( \begin{array}{ccc} \cos \theta_y & 0 & \sin \theta_y \\ 0 & 1 & 0 \\ -\sin \theta_y & 0 & \cos \theta_y \end{array} \right), R_p = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \theta_p & -\sin \theta_p \\ 0 & \sin \theta_p & \cos \theta_p \end{array} \right), R_r = \left( \begin{array}{ccc} \cos \theta_r & -\sin \theta_r & 0 \\ \sin \theta_r & \cos \theta_r & 0 \\ 0 & 0 & 1 \end{array} \right)

このとき、座標系ごとヨー、ピッチ、ロールの順に適用した回転は、元の座標系で R_y R_p R_r で表されます (適用の順序が逆になる)。
証明は次のようにできますが、実際の回転をイメージするとわかりやすいと思います。

証明

ベクトル {\bf x} にヨーを作用させると、
R_y {\bf x}
次に、これにピッチを作用させるには、いったん座標系を戻して R_p を作用させるから、
R_y R_p R_y^{-1} \cdot R_y {\bf x} = R_y R_p {\bf x}
同様に、これにロールを作用させると、
(R_y R_p) R_r (R_y R_p)^{-1} \cdot R_y R_p {\bf x} = R_y R_p R_r {\bf x}
(証明終)

WPF での実際の演算では四元数を使うとよいでしょう。

■ 2 点の座標 → オイラー角

ここでは、2 点として {\bf e}_z = (0, 0, 1), {\bf e}_y = (0, 1, 0) を選びます。
これらの回転後のベクトルがそれぞれ {\bf u}, {\bf t} で与えられたとすると、オイラー角は以下の手続きにより求められます。

まず、{\bf u} の x 要素と z 要素の比はピッチおよびロールの影響を受けないから、
\tan \theta_y = \dfrac{{\bf u}_x}{{\bf u}_z}
により、\theta_y, R_y が決まる。

ピッチおよびロールの決め方から、
\tan \theta_p = \dfrac{- (R_p {\bf e}_z)_y}{(R_p {\bf e}_z)_z}, \tan \theta_r = \dfrac{- (R_r {\bf e}_y)_x}{(R_r {\bf e}_y)_y}

また、前項と同様に考えて、{\bf u} = R_y R_p {\bf e}_z, {\bf t} = R_y R_p R_r {\bf e}_y であるから、
R_p {\bf e}_z = R_y^{-1} {\bf u}, R_r {\bf e}_y = R_p^{-1} R_y^{-1} {\bf t}

したがって、\theta_p, \theta_r も順に決まる。
なお、この決め方の場合、回転角度の範囲は、
- \pi < \theta_y \le \pi, - \dfrac{\pi}{2} \le \theta_p \le \dfrac{\pi}{2}, - \pi < \theta_r \le \pi
となる。
(終)

arctan を求めるには、Math.Atan2 メソッドを使うとよいでしょう。

以上により、行列 → 2 点の座標 → オイラー角 → 四元数 → 行列と変換する方法が与えられたので、
回転の表現を相互に変換できるようになります。

では、これらを実装してみます。
ベクトル、行列、四元数を扱うための System.Numerics.Vectors という、WPF のライブラリよりも高機能なライブラリもあるのですが、
今回は WPF のライブラリのみを利用して実装したいと思います。

using System;
using System.Diagnostics;
using System.Windows.Media.Media3D;
namespace RotationTest
{
public static class Rotation3DHelper
{
public static readonly Vector3D UnitX = new Vector3D(1, 0, 0);
public static readonly Vector3D UnitY = new Vector3D(0, 1, 0);
public static readonly Vector3D UnitZ = new Vector3D(0, 0, 1);
public static double ToDegrees(double radians) => radians * 180 / Math.PI;
public static double ToRadians(double degrees) => degrees * Math.PI / 180;
public static Quaternion CreateQuaternionInRadians(Vector3D axis, double angleInRadians) => new Quaternion(axis, ToDegrees(angleInRadians));
// 四元数 → 行列
public static Matrix3D ToMatrix3D(this Quaternion q)
{
var m = new Matrix3D();
m.Rotate(q);
return m;
}
// オイラー角 → 四元数
public static Quaternion ToQuaternion(this EulerAngles e) =>
CreateQuaternionInRadians(UnitY, e.Yaw) *
CreateQuaternionInRadians(UnitX, e.Pitch) *
CreateQuaternionInRadians(UnitZ, e.Roll);
// 2 点の座標 → オイラー角
public static EulerAngles ToEulerAngles(Vector3D rotatedUnitZ, Vector3D rotatedUnitY)
{
var y_yaw = Math.Atan2(rotatedUnitZ.X, rotatedUnitZ.Z);
var m_yaw_inv = CreateQuaternionInRadians(UnitY, -y_yaw).ToMatrix3D();
rotatedUnitZ = rotatedUnitZ * m_yaw_inv;
rotatedUnitY = rotatedUnitY * m_yaw_inv;
var x_pitch = Math.Atan2(-rotatedUnitZ.Y, rotatedUnitZ.Z);
var m_pitch_inv = CreateQuaternionInRadians(UnitX, -x_pitch).ToMatrix3D();
rotatedUnitY = rotatedUnitY * m_pitch_inv;
// 本来は -rotatedUnitY.X だけでよいはずです。しかし、X=0 のときに π でなく -π となってしまうため、場合分けします。
var z_roll = Math.Atan2(rotatedUnitY.X == 0 ? 0 : -rotatedUnitY.X, rotatedUnitY.Y);
return new EulerAngles { Yaw = y_yaw, Pitch = x_pitch, Roll = z_roll };
}
}
[DebuggerDisplay(@"\{Yaw={Yaw}, Pitch={Pitch}, Roll={Roll}\}")]
public struct EulerAngles
{
public double Yaw { get; set; }
public double Pitch { get; set; }
public double Roll { get; set; }
}
}
view raw Rotation3DHelper.cs hosted with ❤ by GitHub

全体のソースコードは RotationTest (GitHub) にあります。
これらのメソッドを呼び出すテストが付いています。

前回: WPF で 3D オブジェクトを回転させる
次回: Leap Motion で手の回転状態を取得する

バージョン情報
.NET Framework 4.5

参照

命題論理の複雑な問題を解く (C#)

// この投稿は 数学 Advent Calendar 2016 の 23 日目の記事です。

前回の命題論理を実装する (C#) では命題論理を扱うためのライブラリを作成しました。
今回は、このライブラリを利用して複雑な問題をいくつか解いてみます。

騎士と悪漢

まず、前回でも紹介した書籍 Raymond Smullyan「記号論理学: 一般化と記号化」から、問題を 2 つ引用します。


すべての住民は騎士または悪漢のいずれかであり、騎士はつねに本当のことを言い、悪漢はつねに噓をつく。

問題 1.3
住民 A, B がおり、A は「私たちは二人とも悪漢だ」と言った。
A, B はそれぞれ騎士か悪漢か?

問題 1.5
住民 A, B がおり、A は「私たちは同種 (二人とも騎士または二人とも悪漢のいずれか) だ」と言った。
A, B はそれぞれ騎士か悪漢か?


普通に考えて解くこともできますが、この状況を記号で表すことを考えます。
住民 X が騎士であるという命題を k_X とし、X が命題 p を主張したとすると、
k_Xp の真偽性は等しいということになり、k_X \equiv p が成り立ちます。
したがって、問題 1.3 では k_A \equiv (\lnot k_A \land \lnot k_B) が、問題 1.5 では k_A \equiv (k_A \equiv k_B) が成り立ちます。

前回作成したライブラリ Blaze を使って、次のように解くことができます。

using System;
using Blaze.Propositions;
using static System.Console;
using static Blaze.Propositions.Formula;
namespace PropositionsConsole
{
class Program
{
static void Main(string[] args)
{
Knights_1_3();
Knights_1_5();
}
static void Knights_1_3()
{
var ka = Variable("kA");
var kb = Variable("kB");
var knowledge = Equivalent(ka, !ka & !kb);
knowledge.Determine(ka);
knowledge.Determine(kb);
WriteLine("Q 1.3");
WriteVariable(ka);
WriteVariable(kb);
}
static void Knights_1_5()
{
var ka = Variable("kA");
var kb = Variable("kB");
var knowledge = Equivalent(ka, Equivalent(ka, kb));
knowledge.Determine(ka);
knowledge.Determine(kb);
WriteLine("Q 1.5");
WriteVariable(ka);
WriteVariable(kb);
}
static void WriteVariable(VariableFormula v)
{
WriteLine($"{v}: {(v.Value.HasValue ? v.Value.ToString() : "Null")}");
}
}
}
view raw Program.cs hosted with ❤ by GitHub

実行すると、問題 1.3 では A が悪漢、B が騎士と確定します。
問題 1.5 では B が騎士と確定しますが、A は確定できません。

PropositionsConsole

 

数当てゲーム

さて次に、David Gale「Tracking the Automatic ANT: And Other Mathematical Explorations」という書籍から、
数当てゲームを紹介します。


  • プレイヤーは 2 人
  • それぞれのプレイヤーに数が割り当てられており、自身の数を知っているが、相手の数は知らない
  • 2 人の持つ数は、連続する 2 つの自然数 (つまり、差は 1)
  • ターン制で、相手の数を当てれば終了 (わかったら宣言し、わからなければ発声しない)
  • 先攻・後攻の区別はなく、いずれが宣言してもよい

どちらのプレイヤーが何ターン目に相手の数を当てることができるか?

人型ロボットと数当てゲームをする人のイラスト


自分の数が 9 だとすると、与えられた条件から相手の数は 8 または 10 の 2 通りに絞られます。
ここまではすぐにわかるのですが、いずれかを確定させるには相手の将来の行動あるいは実際の行動の原因を推測する必要があります。
相手がどう考えているかを考える、という日常生活においても重要な能力を問われているようでもあります。

小さい数値でいくつか試してみましょう。
2 人のプレイヤーには対称性がありますが、ここでは大小関係を固定して、プレイヤー A が n、プレイヤー B が n+1 を持つとします。
以下では、各プレイヤーの心理状況を書くことにします。

n=1 の場合、

  • A「B=2 しかありえない。」
  • B「A=1 または A=3。A=1 ならば、A は 1 ターン目に当てるはず。A=3 ならば、A は 1 ターン目に静観するはず。」

したがって、A が 1 ターン目で当てて終了します。
n=2 の場合、

  • A「B=1 または B=3。B=1 ならば、B は 1 ターン目に当てるはず。B=3 ならば、B は 1 ターン目に静観するはず。」
  • B「A=2 または A=4。A=2 ならば、A は 2 ターン目に当てるはず (B が 1 ターン目に静観するため)。
    A=4 ならば、A は 2 ターン目まで静観するはず。」

したがって、1 ターン目は両方静観し、A が 2 ターン目で当てて終了します。
このように、相手がどのように行動してくるかをあらかじめシミュレートしておけば、実際の行動に応じて相手の数を確定できることになります。

この方針でプログラミングしたものが NumberGuessConsole (GitHub) です (コードが長いため、ここには記載しません)。
このプログラムを A=8, B=9 で実行した結果です。A が 8 ターン目に当てました。
(「X @ n」は X が n ターン目に当てることを表します。)

NumberGuessConsole

各プレイヤーの脳内であらかじめシミュレーターを実行することで、値を確定させるための条件を得ます。
A の Knowledge には ((B = 7)≢(B = 9))∧((B = 7)⇒(B @ 7))∧((B = 9)⇒~(B @ 7)) と入っています。
B が 7 ターン目に静観したため、A は B=9 だとわかりました。

NumberGuessConsole

結果として、小さいほうの数 n を持つプレイヤーが n ターン目に相手の数 n+1 を当てることになります。
(プログラミングをしなくても、同様の方針により少し複雑な数学的帰納法で証明できます。)

前回: 命題論理を実装する (C#)

作成したライブラリ
Blaze (GitHub)
Blaze (NuGet Gallery)

作成したサンプル
PropositionsConsole (GitHub)
NumberGuessConsole (GitHub)

素材
不気味の谷現象のイラスト (いらすとや)

バージョン情報
.NET Framework 4.5

参照
Raymond Smullyan「記号論理学: 一般化と記号化」
David Gale「Tracking the Automatic ANT: And Other Mathematical Explorations」

カテゴリー: 数学. タグ: . 1 Comment »