// この投稿は C# Advent Calendar 2014 の 7 日目の記事です。
ゲーム アプリケーションやセンサーを用いたアプリケーションでは、
力学系の方程式の計算や座標系の変換など、何らかの数値計算が必要になることがあります。
中でも一次方程式は頻度が高いと思いますが、開発者の多くは毎回、
あらかじめ手計算により解の形式 (x = ～) を求めてから実装しているのではないでしょうか？
例えば、ax + b = cx + d の形式のものを x = (d – b) / (a – c) に変形してからプログラムするとか。
処理速度に重点を置く場合はこれでよいのですが、
最初に脳内で導出された式とソースコード上の式とでは形式に乖離が生じることがあり、
確認のために再計算しなければならなくなるなど、保守性が高いとはいえないでしょう。
パターン化できるのであれば、面倒な作業は自動化したいところです。
そこで以下では、なるべく元の方程式のままの形から解を得る仕組みを考えます。
まず、数式を C# の標準的な算術演算子を用いて表現できるようにするため、
一変数の多項式を表す Polynomial 構造体を以下のように定義します。
それぞれの次数に対する係数を Dictionary で管理し、必要な演算子を用意します。
(全体のソースコードは GitHub の EquationConsole にあります。)
public struct Polynomial {
public static readonly Polynomial X
= new Polynomial(
new Dictionary<int,
double> { {
1,
1 } });
static readonly IDictionary<int, double> _coefficients_empty = new Dictionary<int, double>();
IDictionary<int, double> _coefficients;
IDictionary<int, double> Coefficients
{
get { return _coefficients == null ? _coefficients_empty : _coefficients; }
}
// The dictionary represents index/coefficient pairs.
// The dictionary does not contain entries whose coefficient is 0.
public Polynomial(IDictionary<int, double> coefficients)
{
_coefficients = coefficients;
}
public static implicit operator Polynomial(double value)
{
return value == 0 ? default(Polynomial) : new Polynomial(new Dictionary<int, double> { { 0, value } });
}
public static Polynomial operator +(Polynomial p1, Polynomial p2)
{
var coefficients = new Dictionary<int, double>(p1.Coefficients);
foreach (var item2 in p2.Coefficients)
{
AddMonomial(coefficients, item2.Key, item2.Value);
}
return new Polynomial(coefficients);
}
public static Polynomial operator –(Polynomial p1, Polynomial p2)
{
var coefficients = new Dictionary<int, double>(p1.Coefficients);
foreach (var item2 in p2.Coefficients)
{
AddMonomial(coefficients, item2.Key, –item2.Value);
}
return new Polynomial(coefficients);
}
public static Polynomial operator *(Polynomial p1, Polynomial p2)
{
var coefficients = new Dictionary<int, double>();
foreach (var item1 in p1.Coefficients)
{
foreach (var item2 in p2.Coefficients)
{
AddMonomial(coefficients, item1.Key + item2.Key, item1.Value * item2.Value);
}
}
return new Polynomial(coefficients);
}
public static Polynomial operator /(Polynomial p, double value)
{
var coefficients = new Dictionary<int, double>();
foreach (var item in p.Coefficients)
{
AddMonomial(coefficients, item.Key, item.Value / value);
}
return new Polynomial(coefficients);
}
// Power
public static Polynomial operator ^(Polynomial p, int power)
{
if (power < 0) throw new ArgumentOutOfRangeException("power", "The value must be non-negative.");
Polynomial result = 1;
for (var i = 0; i < power; i++)
{
result *= p;
}
return result;
}
public static Polynomial operator +(Polynomial p)
{
return p;
}
public static Polynomial operator –(Polynomial p)
{
return –1 * p;
}
static void AddMonomial(Dictionary<int, double> coefficients, int index, double coefficient)
{
if (coefficients.ContainsKey(index))
{
coefficient += coefficients[index];
}
if (coefficient != 0)
{
coefficients[index] = coefficient;
}
else
{
coefficients.Remove(index);
}
}
}
今回は等値演算子 (==) などは省きました。
暗黙的型変換演算子 (implicit operator Polynomial) を定義したため、
多項式 -x + 2 や 2x – 1 を次のように表すことができます。
var x = Polynomial.X;
var p1 = –x + 2;
var p2 = 2 * x – 1;
次に、この Polynomial 構造体に、方程式を解くためのメソッドを追加します。
方程式の左辺はこの多項式自身で、右辺は 0 であるとします。
一次方程式のほか、二次方程式についても実装しています。
public int Degree
{
get {
return Coefficients
.Count
== 0 ? 0 : Coefficients
.Max(c
=> c
.Key); }
}
// Substitution
public double this[double value]
{
get { return Coefficients.Sum(c => c.Value * Math.Pow(value, c.Key)); }
}
// Solve the equation whose right operand is 0.
public double SolveLinearEquation()
{
if (Degree != 1) throw new InvalidOperationException("The degree must be 1.");
// ax + b = 0
var a = GetCoefficient(1);
var b = GetCoefficient(0);
return –b / a;
}
// Solve the equation whose right operand is 0.
public double[] SolveQuadraticEquation()
{
if (Degree != 2) throw new InvalidOperationException("The degree must be 2.");
// ax^2 + bx + c = 0
var a = GetCoefficient(2);
var b = GetCoefficient(1);
var c = GetCoefficient(0);
var d = b * b – 4 * a * c;
return d > 0 ? new[] { (–b – Math.Sqrt(d)) / (2 * a), (–b + Math.Sqrt(d)) / (2 * a) }
: d == 0 ? new[] { –b / (2 * a) }
: new double[0];
}
double GetCoefficient(int index)
{
return Coefficients.ContainsKey(index) ? Coefficients[index] : 0;
}
インデクサーには代入の機能を割り当てました。
以上で Polynomial 構造体の最低限の実装はできたと思います。
これを使うと、例えば直線 y = x – 1 と直線 y = -2x + 5 の交点 P は次のように求められます。
ほぼ数式のままの形で記述できていると思います。
static class Program {
static readonly Polynomial x
= Polynomial.X;
static void Main(string[] args)
{
var l1 = x – 1;
var l2 = –2 * x + 5;
var p_x = (l1 – l2).SolveLinearEquation();
var p_y = l1[p_x];
Console.WriteLine("P({0}, {1})", p_x, p_y); // P(2, 1)
}
}
次の例では、特定の 2 点を通る直線上で、y 座標から x 座標を求めます。
static class Program {
static readonly Polynomial x
= Polynomial.X;
static void Main(string[] args)
{
var p1 = new Point2D(0, –300);
var p2 = new Point2D(1800, 300);
var y_to_x = GetFunc_y_to_x(p1, p2);
Console.WriteLine(y_to_x(0)); // 900
Console.WriteLine(y_to_x(–100)); // 600
}
// P1, P2 を通る直線上で、指定された y 座標に対応する x 座標を求めるための関数。
static Func<double, double> GetFunc_y_to_x(Point2D p1, Point2D p2)
{
// P1(x1, y1) および P2(x2, y2) を通る直線の方程式 (の公式):
// (x – x1) (y2 – y1) – (x2 – x1) (y – y1) = 0
return y => ((x – p1.X) * (p2.Y – p1.Y) – (p2.X – p1.X) * (y – p1.Y)).SolveLinearEquation();
}
}
struct Point2D
{
public double X { get; private set; }
public double Y { get; private set; }
public Point2D(double x, double y)
: this()
{
X = x;
Y = y;
}
}
座標系のマッピングも、簡単なものであればこれとほぼ同様にできます (一般的には行列を使いますが)。
ちなみに、二次方程式 x2 + x -1 = 0 の解もこの通り。
なお、今回は引数が (Polynomial, Polynomial) のみの +, -, * 演算子を定義しましたが、
(Polynomial, double) などの引数を受け付けるオーバーロードを用意すれば、演算量を若干削減できると思います。
作成したサンプル
EquationConsole (GitHub)
参照
オーバーロードされた演算子 (C# プログラミング ガイド)
