C#で Project Euler P9
はじめに
方針を書いてますが、あくまで方針であり証明などは省いてます。
気になる方はgoogle先生に聞くと答えてくれるので聞いてみてください。
問題
Problem 9 「特別なピタゴラス数」
ピタゴラス数(ピタゴラスの定理を満たす自然数)とは a < b < c で以下の式を満たす数の組である.a2 + b2 = c2
例えば, 32 + 42 = 9 + 16 = 25 = 52 である.a + b + c = 1000 となるピタゴラスの三つ組が一つだけ存在する.
これらの積 abc を計算しなさい.
Problem 9 - PukiWiki
方針
原始ピタゴラス数を生成する公式を経由して因数分解の問題に置き換え、互いに素な自然数の性質から探査範囲を絞ることで最適化する。
以下の式ではピタゴラス数の合計値である 1000 を一般化する為に定数 $s$ とする。
はじめに、原始ピタゴラス数 $a_0, b_0, c_0$ は自然数 $m, n$ を用いて次の式で表せる。
\begin{align*}
a_0 =& \ m^2 - n^2 \\
b_0 =& \ 2mn \\
c_0 =& \ m^2+n^2 \\
\end{align*}
ただし、$m, n$ は
\begin{align*}
・&m と n は互いに素 \\
・&m > n \\
・&m, nの偶奇が異なる \\
\end{align*}
を満たす。
また、ピタゴラス数 $a, b, c$ は原始ピタゴラス数の自然数倍であるため、自然数 $k$ を用いて次のように表せる。
\begin{cases}
a =& (ka_0)^2 \\
b =& (kb_0)^2 \ \ \ \ \ \ \ \ ①\\
c =& (kc_0)^2 \\
\end{cases}
これを与式へ代入すると、
\begin{align*}
s =& \ a + b + c \\
=& \ ka_0 + kb_0 + kb_0 \\
=& \ k(a_0 + b_0 + c_0) \\
=& \ k(m^2 - n^2 + 2mn + m^2 + n^2) \\
=& \ k(2m^2 + 2mn) \\
=& \ 2km(m + n) \\
\frac{s}{2} =& \ km(m + n) \\
\end{align*}
となる。
定義より、右辺は自然数となる為 $s$ は偶数でなければ解が無い。
$s$ が偶数の場合を考える。
$\gcd(x, y) = x と y の最大公約数$ として
\begin{align*}
s' =& \frac{s}{2} \\
n' =& \ m + n \\
k_m =& \gcd \bigl( \frac{s'}{m}, m \bigr) \\
k_{n'} =& \gcd(\frac{s'}{n'}, n') \\
\end{align*}
と置くと、$m > n$ より $n + m < m + m $ , $m\ と\ n は互いに素$ より$n'\ と\ mは互いに素$ から $n'$ は
\begin{align*}
2m >& \ n' > m \\
k =& \ k_mk_{n'} \\
s' =& \ k_mk_{n'}mn' \\
n' =& \ \frac{s'}{k_mk_{n'}m}
\end{align*}
を満たす解が存在する場合に解をもち、その時の $m, n, k$ を①に代入すれば、与式を満たすピタゴラス数が得られる。
また、この時の未知数 $m, k_{n'}$ は、${\rm divisor}(x) = \{ d \ | \ x \ {\rm mod} \ d = 0 \}$ を使って
\begin{align*}
m \in & \ {\rm divisor}(s') \\
k_{n'} \in & \ {\rm divisor}\bigl(\frac{s'}{k_m m}\bigr)
\end{align*}
を満たしている。
$|\ {\rm divisor}(x)\ | < \sqrt{x}$ である為、$n'$ を求めるには高々 $O(s)$ 程度の代入を行えばいいことがわかる。
解答
P9.cs
using System; public static class P9 { public static void Main() { Console.WriteLine(Calc()); } private static string Calc() { var res = Get(1000); if(res.Item1 != 0) { var r = GetPythagoreanNumber(res.Item1, res.Item2, res.Item3); return r.ToString() + " : " + (r.Item1 * r.Item2 * r.Item3).ToString(); } return "none"; } private static Tuple<ulong, ulong, ulong> GetPythagoreanNumber(ulong m, ulong n, ulong k) { return Tuple.Create(k * (m * m - n * n), k * 2 * m * n, k * (m * m + n * n)); } private static Tuple<ulong, ulong, ulong> Get(ulong s) { var none = Tuple.Create(0UL, 0UL, 0UL); if((s & 1) == 1) { return none; } var s_ = s / 2; foreach(var m in Multiple.Divisor(s_)) { var k_m = Multiple.Gcd(s_ / m, m); foreach(var k_n in Multiple.Divisor(s_ / (m * k_m))) { var n_ = s_ / (k_m * k_n * m); if(n_ < 2 * m && n_ > m) { return Tuple.Create(m, n_ - m, k_m * k_n); } } } return none; } }
Multiple.cs
using System; using System.Collections.Generic; using System.Linq; public static class Multiple { // 最大公約数を求める public static ulong Gcd(ulong a, ulong b) { return a * b == 0 ? 0 : MyMath.EuclideanAlgorithm(a, b); } // 最小公倍数を求める public static ulong Lcm(ulong a, ulong b) { return a / Gcd(a, b) * b; } // 素因数分解をする public static List<Tuple<ulong, ulong>> Factor(ulong n) { var ret = new List<Tuple<ulong, ulong>>(); for(ulong primeIndex = 0; n > 1; primeIndex++) { var prime = Prime.GetPrime(primeIndex); // ここで除数を指数的に増やせばO(log n)に if(n % prime == 0) { ulong primeCounter = 0; while(n % prime == 0) { n /= prime; primeCounter++; } ret.Add(Tuple.Create(prime, primeCounter)); } } return ret; } public static List<ulong> LinerFactor(ulong n) { return FactorToLiner(Factor(n)); } public static List<ulong> FactorToLiner(List<Tuple<ulong, ulong>> factor) { var ret = new List<ulong>(); foreach(var t in factor) { for(ulong i = t.Item2; i > 0; i--) { ret.Add(t.Item1); } } return ret; } // 素因数分解された集合から自然数を得る public static ulong ReverceFactor(List<Tuple<ulong, ulong>> factor) { return factor.Aggregate(1UL, (multi, t) => multi * MyMath.Pow(t.Item1, t.Item2)); } // 素因数分解の一意性を用いた掛け算 public static List<Tuple<ulong, ulong>> FactorMulti(List<Tuple<ulong, ulong>> a, List<Tuple<ulong, ulong>> b) { var ret = new List<Tuple<ulong, ulong>>(); var a_ = a.ToArray(); // O(n) var b_ = b.ToArray(); if(a.Count == 0) { return b; } else if(b.Count == 0) { return a; } if(a_[a_.Length - 1].Item1 < b_[b_.Length - 1].Item1) { var work = b_; b_ = a_; a_ = work; } int aIndex = 0; int bIndex = 0; // O(n) while(bIndex < b_.Length) { if(a_[aIndex].Item1 == b_[bIndex].Item1) { var p = a_[aIndex].Item1; ret.Add(Tuple.Create(p, a_[aIndex++].Item2 + b_[bIndex++].Item2)); } else if(a_[aIndex].Item1 < b_[bIndex].Item1) { ret.Add(a_[aIndex++]); } else { ret.Add(b_[bIndex++]); } } while(aIndex < a_.Length) { ret.Add(a_[aIndex++]); } return ret; } // デバッグ用、素因数分解の結果を表示する public static void ShowFactor(List<Tuple<ulong, ulong>> factor) { foreach(var t in factor) { Console.Write(t.Item1 + "^" + t.Item2 + ", "); } Console.Write("\r\n"); } // 約数を取得 O(sqrt(n)) public static IEnumerable<ulong> Divisor(ulong n) { var desk = new Stack<ulong>(); ulong i; for(i = 1; i * i < n; i++) { if(n % i == 0) { yield return i; desk.Push(n / i); } } if(i * i == n) { yield return i; } foreach(var d in desk) { yield return d; } } }
計算量
本プログラム：$O(n)$
総当たり：$O(n^2)$
n : ピタゴラス数の和
