累乗和の公式をプログラムで導出してみる[C#]

　水曜日からずっとネットが繋がらなくなっていた。知りたいことが知れないというのは，昔は当たり前だったが，今は当然になりすぎていて苦痛だった。さて本題本題。

1.累乗和の公式とは

べき乗和の公式とも呼ばれ，高校生の数列で習う内容もので，整数の足し算に関する公式です。深く説明する前に次の式を見てください。

この式は数学好きな小・中学生は知っているかもしれないですね。有名な話だとガウスが1から100までの合計を一瞬で計算したという話がありますね。実際にn=100を入れてみるとS1=5050となります。

　同様にして次のような式もあります。二乗和の公式と呼ばれるものです。

さらにさらに三乗和の公式というものもあります。

このように，自然数のn乗の和を求める公式のことを総じて累乗和の公式と呼ぶらしいです。

2.累乗和の公式の求め方

上の式をみて，なんとかして一般化して，n乗和の公式を求める方法を求めたいと思うかと思います。その方法はいくつかありますが，高校数学などでは(n+1)^k-n^kを計算することで求める方法が紹介されています。

　またベルヌーイ数を用いて一般化した式も存在するらしいです。詳しくはこちらのWikipediaとか参照されたし。

　上の式の説明をなぜ省いたかというと，今日作ったプログラムでは別の方法で実装しているからです。しかも個人的には上の二つよりもわかりやすい方法だと思っています（大げさ）。

　まず一乗和の公式と二乗和の公式を見てみてください。

見てみるとわかると思いますが，一乗和の公式はnの二次式，二乗和の公式はnの三次式となっています。このことから，

k乗和の公式はnのk+1次式で表わされるということが予想できます。これの証明は(n+1)^k-n^kによってk乗和の公式が求められるということから証明は出来ると思います。

　で，上の式をどう料理するかという話ですが，未定係数aの連立方程式を作ればいいわけです。わかりづらいと思うので三乗和の公式，つまりk=3で説明します。
　
　今，上の式より，S_3はnの4次式で表すことができるから，

とあらわすことができます。また，S3は三乗和の公式であるから，

　となる。したがって，

という4元連立方程式を立てることができます。これをa1，a2，a3，a4について解けば3乗和の公式の導出ができるという算段です。

　ただ自分で考えたやり方なので正式な名称は知らない。でもこういう初歩的な数学って誰かが絶対にやってるというお約束があるから正式名称を知ってる人がいたら教えて欲しいところです。

3.プログラム

連立方程式を解けばいいと書きましたが，その方法としてはガウスの消去法と呼ばれる方法を用いています。ただ下で書いたプログラムでは累乗和の公式を求めるだけだから，と一般化して書いていない（具体的には解がない連立方程式を入れると挙動がおかしくなる）ので注意されたし。

　あと浮動小数点を用いているので分数にはならなりません。なのでN=100とかはオーバーフローして求められないし，N=3であっても1/6は0.166667と表示されます。そのへんはまた機会があったら修正しようと思います。

/**
 * 2015/10/10
 * 連立方程式で累乗和の公式を求めるプログラム
 */
using System;

class Program
{
   static void Main()
   {
      Console.WriteLine("何乗和の公式を求めますか？");
      var N = int.Parse(Console.ReadLine()) + 1;
      // 係数を求める連立方程式の拡大係数行列を作る
      var formula = MakeFormula(N);
      // 連立方程式をガウス消去法で計算
      Solve( formula , N);
      // 結果を表示
      Show(formula,N);
   }
   
   // 連立方程式を解くメソッド
   static void Solve(double[,] formula,int N)
   {
      double temp;
      
      // ガウス消去法
      
      // 前進消去
      for(int i = 0 ; i < N ; i++ )
      {
         // 対角対角成分をを1にする
         temp = formula[i,i];
         for(int j = 0 ; j < N + 1 ; j++ )
            formula[i,j] /= temp;
         
         // i+1行以降のi列目を0にする
         for(int j = i+1 ; j < N ; j++ )
         {
            temp = formula[j,i];
            for(int k = 0 ; k < N+1 ; k++ )
               formula[j,k] -= temp * formula[i,k];
         }
      }
      
      // 後進消去
      for(int i = N-1 ; i >= 0 ; i-- )
      {
         // 対角成分以外を0にする
         for(int j=i-1;j>=0;j--)
         {
            formula[j,N] -= formula[j,i] * formula[i,N];
            formula[j,i] = 0;
         }
      }
   }
   
   // 累乗和の公式を表示するメソッド
   static void Show(double[,] d,int N)
   {
      Console.WriteLine("{0}乗和の公式は……",N-1);
      for(int i = 0 ; i < N ; i++ )
      {
         Console.Write("({0:0.0000})*n^{1}",d[i,N],i+1);
         if( i + 1 < N )
            Console.Write("+");
      }
      Console.WriteLine();
   }
   
   // 累乗和の公式を求める連立方程式の拡大係数行列を作るメソッド
   static double[,] MakeFormula(int N)
   {
      var ret = new double[N,N+1];
      var sum = 0.0;
      
      for(int i = 0 ; i < N ; i++ )
      {
         // 係数部分を作る
         for(int j = 0 ; j < N ;j++ )
            ret[i,j] = Math.Pow(i+1,j+1);
         
         // 答えの部分を作る
         sum = sum + Math.Pow(i+1,N-1);
         ret[i,N] = sum;
      }
      
      return ret;
   }
}

実行結果。

　実行すると余計分数で書いて欲しいと思いたくなるね。

追記：
　分数で実装してみました。→累乗和の公式をプログラムで実装してみる（分数編）[C#]
　255乗和の公式のような，任意のべき乗和を求めるプログラムが完成しました→BigIntegerを使って累乗和の公式を計算[C#]