さんぽみち

なにか思いついたときとかに気まぐれに更新されます。

C#で Project Euler P11

問題

上の 20×20 の格子のうち, 斜めに並んだ4つの数字が赤くマークされている.

08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48
それらの数字の積は 26 × 63 × 78 × 14 = 1788696 となる.

上の 20×20 の格子のうち, 上下左右斜めのいずれかの方向で連続する4つの数字の積のうち最大のものはいくつか?
Problem 11 - PukiWiki

方針

C#で Project Euler P8 - さんぽみちが拡張された問題。
よって、素直にこれを2次元へ拡張する。

解答

using System;
using System.Collections.Generic;
using System.Linq;
using System.IO;

public static class P11
{
    public static void Main()
    {
        Console.WriteLine(Calc());
    }
    
    private const int MultiDepth = 4;
    private static Dictionary<Tuple<int, int, int>, Tuple<ulong, ulong, ulong, ulong>> multiMap = new Dictionary<Tuple<int, int, int>, Tuple<ulong, ulong, ulong, ulong>>();
    private static ulong Calc()
    {
        string rawSquare;
        using(var sr = new StreamReader("P11_source.txt"))
        {
            rawSquare = sr.ReadToEnd();
        }
        var square = rawSquare
            .Replace("\r\n", "\n")
            .Split('\n')
            .Select(line => line
                .Split(' ')
                .Select(elem => Convert.ToUInt64(elem))
                .ToArray()
            ).ToArray();
        MakeMap(square, MultiDepth);
        
        ulong max = 0;
        for(var y = 0; y < square.Length - MultiDepth; y++)
        {
            for(var x = MultiDepth; x < square[0].Length - MultiDepth; x++)
            {
                if(multiMap.ContainsKey(Tuple.Create(MultiDepth, x, y)))
                {
                    var val = multiMap[Tuple.Create(MultiDepth, x, y)];
                    max = val.Item1 > max ? val.Item1 : max;
                    max = val.Item2 > max ? val.Item2 : max;
                    max = val.Item3 > max ? val.Item3 : max;
                    max = val.Item4 > max ? val.Item4 : max;
                }
            }
        }
        return max;
    }
    
    private static void MakeMap(ulong[][] square, int depth)
    {
        MakeMapD1(square);
        for(int depthWork = depth >> 1, depthIndex = 2; depthWork > 0; depthWork >>= 1, depthIndex <<= 1)
        {
            var beforeDepth = depthIndex >> 1;
            for(var y = 0; y < square.Length - 1 - depthIndex; y++)
            {
                for(var x = 0; x < depthIndex; x++)
                {
                    var before = multiMap[Tuple.Create(beforeDepth, x, y)];
                    SetMap(
                        depthIndex, x, y,
                        before.Item1 * multiMap[Tuple.Create(beforeDepth, x + beforeDepth, y)].Item1,
                        before.Item2 * multiMap[Tuple.Create(beforeDepth, x, y + beforeDepth)].Item2,
                        before.Item3 * multiMap[Tuple.Create(beforeDepth, x + beforeDepth, y + beforeDepth)].Item3,
                        0
                    );
                }
                for(var x = depthIndex; x < square[0].Length - 1 - depthIndex; x++)
                {
                    var before = multiMap[Tuple.Create(beforeDepth, x, y)];
                    SetMap(
                        depthIndex, x, y, 
                        before.Item1 * multiMap[Tuple.Create(beforeDepth, x + beforeDepth, y)].Item1,
                        before.Item2 * multiMap[Tuple.Create(beforeDepth, x, y + beforeDepth)].Item2,
                        before.Item3 * multiMap[Tuple.Create(beforeDepth, x + beforeDepth, y + beforeDepth)].Item3,
                        before.Item4 * multiMap[Tuple.Create(beforeDepth, x - beforeDepth, y + beforeDepth)].Item4
                    );
                }
                for(var x = square[0].Length - 1 - depthIndex; x < square[0].Length; x++)
                {
                    var before = multiMap[Tuple.Create(beforeDepth, x, y)];
                    SetMap(
                        depthIndex, x, y,
                        0,
                        before.Item2 * multiMap[Tuple.Create(beforeDepth, x, y + beforeDepth)].Item2,
                        0,
                        before.Item4 * multiMap[Tuple.Create(beforeDepth, x - beforeDepth, y + beforeDepth)].Item4
                    );
                }
            }
            for(var y = square.Length - 1 - depthIndex; y < square.Length; y++)
            {
                for(var x = 0; x < square[0].Length - 1 - depthIndex; x++)
                {
                    var before = multiMap[Tuple.Create(beforeDepth, x, y)];
                    SetMap(
                        depthIndex, x, y, 
                        before.Item1 * multiMap[Tuple.Create(beforeDepth, x + beforeDepth, y)].Item1,
                        0,
                        0,
                        0
                    );
                }
            }
        }
    }
    private static void MakeMapD1(ulong[][] square)
    {
        for(var y = 0; y < square.Length - 1; y++)
        {
            SetMap(1, 0, y, square[y][0], square[y][0], square[y][0], 0);
            for(var x = 1; x < square[0].Length - 1; x++)
            {
                SetMap(1, x, y, square[y][x], square[y][x], square[y][x], square[y][x]);
            }
            SetMap(1, square[0].Length - 1, y, 0, square[y][square[0].Length - 1], 0, square[y][square[0].Length - 1]);
        }
        for(var x = 0; x < square[0].Length - 1; x++)
        {
            SetMap(1, x, square.Length - 1, square[square.Length - 1][x], 0, 0, 0);
        }
    }
    private static void SetMap(int depth, int x, int y, ulong h, ulong v, ulong s, ulong b)
    {
        multiMap.Add(Tuple.Create(depth, x, y), Tuple.Create(h, b, s, b));
    }
}

計算量
本プログラム：$O(n^2 \log m)$
総当たり：$O(n^2m)$
$ n $：格子の一辺あたりの長さ
$ m $：乗じる数字の長さ
一言だけ言いたいのが、オーダー記法では反映されなくて悲しいが、1つの層を処理するたびに計算する必要のある数字の一辺の長さは $2^{(階層)}$ ずつ減っていくんだよ！
$ \log n $ まで計算すると、斜面が対数関数の形をしたピラミッドの頂上となって数字は一つしか残っていないのを考えるとわかりやすいよ！

C#で Project Euler P10

C# Project Euler アルゴリズム

問題

Problem 10 「素数の和」
10以下の素数の和は 2 + 3 + 5 + 7 = 17 である.

200万以下の全ての素数の和を求めよ.

Problem 10 - PukiWiki

方針

C#で Project Euler P3 - さんぽみち で作成した素数生成器を使って足していくだけ。

解答

using System;
using System.Linq;

public static class P10
{
    public static void Main()
    {
        Console.WriteLine(Calc());
    }
    
    private static ulong Calc()
    {
        return PrimeGenerator
            .GetPrimeList((index, prime) => prime < 2000000)
            .Aggregate((sum, prime) => sum + prime);
    }
}

PrimeGeneratorのソースは割愛。
計算量もあんまり書く意味がない気がするので割愛。
素数関係の問題は生成するのがメインで一回素数生成器作っちゃうと楽しめない。

C#で Project Euler P9

C# Project Euler アルゴリズム

はじめに

　方針を書いてますが、あくまで方針であり証明などは省いてます。
　気になる方は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 : ピタゴラス数の和

C#で Project Euler P8

C# Project Euler アルゴリズム

はじめに

ソースの前に、簡単な方針を今回から書きます。

問題

次の1000桁の数字のうち, 隣接する4つの数字の総乗の中で, 最大となる値は, 9 × 9 × 8 × 9 = 5832である.
　
73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450
　
この1000桁の数字から13個の連続する数字を取り出して, それらの総乗を計算する. では、それら総乗のうち、最大となる値はいくらか.
　
EX 6桁の数123789から5個の連続する数字を取り出す場合, 12378と23789の二通りとなり, 後者の23789=3024が最大の総乗となる.
Problem 8 - PukiWiki

方針

構造を考慮して最適化する。
隣接した数字を乗じた数字は使いまわされる。
よって、各数字について右隣の数字と掛け合わせてこれを覚えておく。
上の手順で覚えた数字についても同じように右隣の数字と掛け合わせて覚えておく。
この手順を繰り返して作った数字列を層として考えると、
　1つの数字が掛け合わされた層(与えられた数字列の状態)
　2つの数字が掛け合された層
　4つの数字が掛け合された層
　8つの数字が掛け合された層
と、指数的に数字が掛け合わされ、木構像をとる。
その後、与えられた長さを2進数表記(2のべき乗の和に分解)し、ビットが立っている層のみを要素となった数字が重複しないようにずらして掛け合わせると、与えられた長さを乗じた値が得られる。
今回であれば長さは13であるため、
　$13_{(10)} = 1101_{(2)}$
となり、先頭の13文字で考えると
　$(7) × (3×1×6×7) × (1×7×6×5×3×1×3×3)$
と分割して考える。

解答

using System;
using System.Collections.Generic;

public static class P8
{
    public static void Main()
    {
        Console.WriteLine(Calc());
    }
    
    private const int MultiDepth = 13;
    private static Dictionary<Tuple<int, int>, long> multiMap = new Dictionary<Tuple<int, int>, long>();
    private static long Calc()
    {
        MakeMap(MultiDepth);
        int index = 0;
        long max = 0;
        while(multiMap.ContainsKey(Tuple.Create(MultiDepth, index)))
        {
            if(multiMap[Tuple.Create(MultiDepth, index)] > max)
            {
                max = multiMap[Tuple.Create(MultiDepth, index)];
            }
            index++;
        }
        return max;
    }
    
    
    private static void MakeMap(int depth)
    {
        var txt = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450";
        
        for(var lIndex = 0; lIndex < txt.Length; lIndex++)
        {
            multiMap.Add(Tuple.Create(1, lIndex), (long)((byte)(txt[(int)lIndex]) & 0x0f));
        }
        for(int depthWork = depth >> 1, depthIndex = 2; depthWork > 0; depthWork >>= 1, depthIndex <<= 1)
        {
            var beforeDepth = depthIndex >> 1;
            for(int lIndex = 0; lIndex < txt.Length - depthIndex; lIndex++)
            {
                multiMap.Add(
                    Tuple.Create(depthIndex, lIndex),
                    multiMap[Tuple.Create(beforeDepth, lIndex)] *
                    multiMap[Tuple.Create(beforeDepth, lIndex + beforeDepth)]
                );
            }
        }
        
        if(!multiMap.ContainsKey(Tuple.Create(depth, 0)))
        {
            for(var lIndex = 0; lIndex < txt.Length - depth; lIndex++)
            {
                multiMap.Add(Tuple.Create(depth, lIndex), 1);
            }
            var offset = 0;
            for(int depthWork = depth, depthIndex = 1; depthWork > 0; depthWork >>= 1, depthIndex <<= 1)
            {
                if((depthWork & 1) == 1)
                {
                    for(int lIndex = 0; lIndex < txt.Length - depth; lIndex++)
                    {
                        multiMap[Tuple.Create(depth, lIndex)] *= multiMap[Tuple.Create(depthIndex, lIndex + offset)];
                    }
                    offset += depthIndex;
                }
            }
        }
    }
}

計算量
本プログラム：$O(n\log m)$
総当たり：$O(nm)$
n : 数字列の長さ
m : 乗じる数字の長さ

C#で Project Euler P7

C# Project Euler アルゴリズム

問題

Problem 7 「10001番目の素数」
素数を小さい方から6つ並べると 2, 3, 5, 7, 11, 13 であり, 6番目の素数は 13 である.
10 001 番目の素数を求めよ.
Problem 7 - PukiWiki

解答

using System;

public static class P7
{
    public static void Main()
    {
        Console.WriteLine(Calc());
    }
    
    private static ulong Calc()
    {
        return PrimeGenerator.GetPrime(10000);  // インデックスは 0 始まり
    }
}

計算量　$O(nm(\log n)^2)$
C#で Project Euler P3 - さんぽみち で作った素数生成器を用いて10001番目の素数を得るだけ。
素数生成器のソースは割愛。

C#で Project Euler P6 別解

C# Project Euler アルゴリズム

はじめに

　C#で Project Euler P6 - さんぽみち　の効率のいい方法を見つけので実装しました。
　どう考えてももっと早くなりそうだったのでちゃんと考えました。
　手続的な実装というよりも一般項を求めることがメインです。

問題

Problem 6 「二乗和の差」
最初の10個の自然数について, その二乗の和は,
12 + 22 + ... + 102 = 385
最初の10個の自然数について, その和の二乗は,
(1 + 2 + ... + 10)2 = 3025
これらの数の差は 3025 - 385 = 2640 となる.
同様にして, 最初の100個の自然数について二乗の和と和の二乗の差を求めよ.
Problem 6 - PukiWiki

解答

public static class P6
{
    public static void Main()
    {
        Console.WriteLine(Calc());
    }
    
    private static long Calc()
    {
        return GetPowSub(100);
    }
    private static long GetPowSub(long n)
    {
        return (((3 * n + 2) * n - 3) * n - 2) * n / 12;
    }
}

計算量　$O(1)$

二乗の和を式１、和の二乗を式２とすると、 \begin{align*} 式２=&(N_1 + N_2 + N_3 + ... + N_n)^2 \\ = &N_1N_1+N_1N_2+N_1N_3+...+N_1N_n+N_2N_1+N_2N_2+N_2N_3+...+N_2N_n+ \\ &N_3N_1+N_3N_2+N_3N_3+...+N_3N_n+...+N_nN_1+N_nN_2+N_nN_3+...+N_nN_n \\ = &(N_1^2+N_2^2+N_3^2+...+N_n^2)+2\{N_1(N_2+N_3+...+N_n)+N_2(N_3+...+N_n)+...N_{n-1}(N_n)\}\\ =&Σ_{k=1}^n N_k^2 + 2Σ_{k=1}^{n-1}N_kΣ_{j=k+1}^{n}j \end{align*} ここで、 \begin{align*} 式１ = &Σ_{k=1}^n N_k^2 \end{align*} であるため、 \begin{align*} 式２-式１ = &Σ_{k=1}^n N_k^2 + 2Σ_{k=1}^{n-1}N_kΣ_{j=k+1}^{n}N_j - Σ_{k=1}^n N_k^2\\ =&2Σ_{k=1}^{n-1}N_kΣ_{j=k+1}^{n}N_j\\ =&2Σ_{k=1}^{n-1}k\frac{(k+1+n)}{2}(n-k)\\ =&Σ_{k=1}^{n-1}k(k+n+1)(-k+n)\\ =&Σ_{k=1}^{n-1}k(-k^2+n^2-k+n)\\ =&Σ_{k=1}^{n-1}(-k^3-k^2+kn^2+kn)\\ =&-Σ_{k=1}^{n-1}k^3-Σ_{k=1}^{n-1}k^2+Σ_{k=1}^{n-1}k(n^2+n)\ \ \ \ \ \ \ \ \ \ \ \ ...① \end{align*} と表せる。

まず、$-Σ_{k=1}^{n-1}k^3$について考える。 \begin{align*} -Σ_{k=1}^{n-1}k^3 =&-(Σ_{k=1}^{n}k^3 - n^3)\\ =&-Σ_{k=1}^{n}k^3 + n^3\\ =&-\{\frac{n}{2}(n+1)\}^2 + n^3\\ =&-\frac{n^2}{4}(n^2+2n+1) + n^3\\ =&-\frac{n^2}{4}(n^2+2n+1-4n)\\ =&-\frac{n^2}{4}(n^2-2n+1)\ \ \ \ \ \ \ \ \ \ \ \ ...② \end{align*} 次に、$-Σ_{k=1}^{n-1}k^2$について考える。 \begin{align*} -Σ_{k=1}^{n-1}k^2 =&-(Σ_{k=1}^{n}k^2 - n^2)\\ =&-Σ_{k=1}^{n}k^2 + n^2\\ =&-\frac{n}{6}(n+1)(2n+1) + n^2\\ =&-\frac{n}{6}(2n^2 + 3n + 1) + n^2\\ =&-\frac{n}{6}(2n^2 + 3n + 1 - 6n)\\ =&-\frac{n}{6}(2n^2 - 3n + 1)\ \ \ \ \ \ \ \ \ \ \ \ ...③ \end{align*} 最後に、$Σ_{k=1}^{n-1}k(n^2+n)$について考える。 \begin{align*} Σ_{k=1}^{n-1}k(n^2+n) =&(n^2+n)Σ_{k=1}^{n-1}k\\ =&(n^2+n)\frac{(1+n-1)}{2}(n-1) \\\ =&\frac{n}{2}(n^2+n)(n-1) \\\ =&\frac{n^2}{2}(n+1)(n-1) \\\ =&\frac{n^2}{2}(n^2-1) \ \ \ \ \ \ \ \ \ \ \ \ ...④ \end{align*} ①=②+③+④ であるため、 \begin{align*} &-Σ_{k=1}^{n-1}k^3-Σ_{k=1}^{n-1}k^2+Σ_{k=1}^{n-1}k(n^2+n)\\ =&-\frac{n^2}{4}(n^2-2n+1) -\frac{n}{6}(2n^2 - 3n + 1)+\frac{n^2}{2}(n^2-1)\\ =&\frac{n}{12}(-3n^3+6n^2-3n) +\frac{n}{12}(-4n^2 + 6n - 2)+\frac{n}{12}(6n^3-6n)\\ =&\frac{n}{12}(-3n^3+6n^2-3n-4n^2 + 6n - 2+6n^3-6n)\\ =&\frac{n}{12}(3n^3+2n^2-3n - 2)\\ =&(((3n+2)n-3)n - 2)n\frac{1}{12} \\ \end{align*} となり、一般項が求まる。

2016/05/09追記

　上記の方法でも導出できるが冷静に考えると $(Σ_{k=1}^nk)^2-Σ_{k=1}^nk^2$ を、等差数列の和及び平方数の和の公式から簡単にするだけでよかった気がした、両方の公式を上でも使ってるし。
　その方針でやってみると、 \begin{align*} (Σ_{k=1}^nk)^2=&(\frac{n + 1}{2}n)^2 \\ =&\frac{n^2 + 2n + 1}{4}n^2 \\ =&\frac{n^2}{4}(n^2 + 2n + 1) \\\\ Σ_{k=1}^nk^2=&\frac{n}{6}(n+1)(2n+1) \\ =&\frac{n}{6}(2n^2 + 3n + 1) \\\\ (Σ_{k=1}^nk)^2-Σ_{k=1}^nk^2 =& \frac{n^2}{4}(n^2 + 2n + 1) - \frac{n}{6}(2n^2 + 3n + 1) \\ =& \frac{n}{12}(3n^3 + 6n^2 + 3n) + \frac{n}{12}(-4n^2 - 6n - 2) \\ =& \frac{n}{12}(3n^3 + 6n^2 + 3n - 4n^2 - 6n - 2)\\ =& \frac{n}{12}(3n^3 + 2n^2 - 3n - 2)\\ =& (((3n + 2)n - 3)n - 2)n\frac{1}{12} \end{align*} 　前回のから引き継いで書いてたら出発点を間違えてえらい遠回りに。。。

C#で Project Euler P6

C# Project Euler アルゴリズム

問題

Problem 6 「二乗和の差」
最初の10個の自然数について, その二乗の和は,
12 + 22 + ... + 102 = 385
最初の10個の自然数について, その和の二乗は,
(1 + 2 + ... + 10)2 = 3025
これらの数の差は 3025 - 385 = 2640 となる.
同様にして, 最初の100個の自然数について二乗の和と和の二乗の差を求めよ.
Problem 6 - PukiWiki

解答

using System;

public static class P6
{
    public static void Main()
    {
        Console.WriteLine(Calc());
    }
    
    private static long Calc()
    {
        long res = 0;
        long sum = 0;
        for(int a = 100 - 1; a >= 1; a--)
        {
            sum += a + 1;
            res += sum * a;
        }
        return res * 2;
    }
}

計算量　$O(n)$

二乗の和を式１、和の二乗を式２とすると、 \begin{align*} 式２=&(N_1 + N_2 + N_3 + ... + N_n)^2 \\ = &N_1N_1+N_1N_2+N_1N_3+...+N_1N_n+N_2N_1+N_2N_2+N_2N_3+...+N_2N_n+ \\ &N_3N_1+N_3N_2+N_3N_3+...+N_3N_n+...+N_nN_1+N_nN_2+N_nN_3+...+N_nN_n \\ = &(N_1^2+N_2^2+N_3^2+...+N_n^2)+2\{N_1(N_2+N_3+...+N_n)+N_2(N_3+...+N_n)+...N_{n-1}(N_n)\}\\ =&Σ_{k=1}^n N_k^2 + 2Σ_{k=1}^{n-1}N_kΣ_{j=k+1}^{n}j \end{align*} ここで、 \begin{align*} 式１ = &Σ_{k=1}^n N_k^2 \end{align*} であるため、 \begin{align*} 式２-式１ = &Σ_{k=1}^n N_k^2 + 2Σ_{k=1}^{n-1}N_kΣ_{j=k+1}^{n}N_j - Σ_{k=1}^n N_k^2\\ =&2Σ_{k=1}^{n-1}N_kΣ_{j=k+1}^{n}N_j \end{align*} と表せる。
この形で実装しやすいように、$2Σ_{k=1}^{n-1}N_kΣ_{j=k+1}^{n}N_j$を降順にした。