2015年12月6日日曜日

三角関数のテイラー展開をmatplotlibで見てみる[Python]

久々の投稿。matplotlibの勉強にやってみた。マクローリン展開と言ったほうが良いかもしれないけれど。

1.描いたグラフ

グラフだけ見たいという人向けにこちらをどうぞ。-5 <= x <= 5の範囲でやってみました。

当然といえば当然だけれども,一度は自分で描いてみるとすうがくのちからってすげーってなる。近似する次数が大きくなれば大きくなるほど目的の関数に近づいてくるのは不思議な感覚。

2.プログラムの予備知識

とりあえず数式はこいつら。ちゃんとした奴はWikipediaさんとか見てください。
まぁこれをPythonで書くとこんな感じ。sin関数の近似のところはコメント多めにしといた。そこのコメントを見ればnumpyの関数の意味もある程度わかると思う。
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np

pi = np.pi

# サイン
# -5から5までの数字を100個の数字でプロットする
x = np.linspace(-5,5,num=100)
plt.plot(x,np.sin(x),label="sin(x)")
# 1次の項
y = x
plt.plot(x,y,label="n=1")
# 3次の項
y = y - x**3 / 6 
plt.plot(x,y,label="n=3")
# 5次の項
y = y + x**5 / 120
plt.plot(x,y,label="n=5")
# 7次の項
y = y - x**7 / 5040
plt.plot(x,y,label="n=7")
# スケールを指定
plt.xlim(-5,5)
plt.ylim(-5,5)
# タイトルを指定
plt.title("sine approximation formula")
# x=0,y=0の線を書く
plt.axhline(0, color='black')
plt.axvline(0, color='black')
# ラベルを指定
plt.xlabel("x")
plt.ylabel("y")
# 判例を表示
plt.legend()
plt.show()

# コサイン
x = np.linspace(-5,5,num=100)
plt.plot(x,np.cos(x),label="cos(x)")
y = x**0
plt.plot(x,y,label="n=0")
y = y - x**2/2
plt.plot(x,y,label="n=2")
y = y + x**4/24
plt.plot(x,y,label="n=4")
y = y - x**6/720
plt.plot(x,y,label="n=6")
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.title("cosine approximation formula")
plt.axhline(0, color='black')
plt.axvline(0, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

# タンジェント
x = np.linspace(-5,5,num=100)
plt.plot(x,np.tan(x),label="tan(x)")
y = x
plt.plot(x,y,label="n=0")
y = y + x**3 / 3
plt.plot(x,y,label="n=3")
y = y + x**5 * 2.0/15
plt.plot(x,y,label="n=5")
y = y + x**7 * 17.0/315
plt.plot(x,y,label="n=7")
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.title("tangent approximation formula")
plt.axhline(0, color='black')
plt.axvline(0, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
matplotlibありがたし。

3.おまけ

記事のタイトルをマクローリン展開とすべきかテイラー展開とすべきかGoogleトレンドで調べてみた。

7月は期末テストかな。大学生も大変だー。

2015年10月27日火曜日

画像のORB特徴量を抽出[Python]

1.実行結果

実行してみた結果を載せてみる。とりあえず定番Lenaさん。左右反転して横にしてみたが,まぁそこそこの精度で認識できている。
本当は詳しい説明などをしたいのだが,残念ながらそんな知識はないのでとりあえず公式サイト最強ということでごまかす。リンク先は英語です。

2.ソースコード

まぁコメントは大目に書いたつもりだけれど,「画像を読み込む」→「特徴量抽出」→「類似する特徴量を探す」→「類似する部分を描く」という感じ。類似する部分を描くという作業は実際に何かを作るときは必要ないのだけれども,なんかそれ専用のdrawMatches関数なるものが定義されているので使ってみた。
# -*- coding: utf-8 -*-
import numpy as np
import cv2

# ORB特徴量を抽出するためのもの
orb = cv2.ORB_create()

# 画像を読み込み
img1 = cv2.imread("img1.jpg")
img2 = cv2.imread("img2.jpg")

# 各画像のキーとなる場所(keypoint)とその詳細(description)を計算
kp1, des1 = orb.detectAndCompute(img1,None)
kp2, des2 = orb.detectAndCompute(img2,None)

# 特徴量がどれだけ似ているかを調べるためのもの
bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True)

# 一致している(類似している)ものを判定
matches = bf.match(des1,des2)

# 似ている点を線で引いてくれる関数
img3 = cv2.drawMatches(img1,kp1,img2,kp2,matches[:10],None, flags=2)

# 結果の保存
cv2.imwrite("match.jpg",img3)

特定の色の物体の座標を検出[Python]

1. 概要

PythonのOpenCVを使って特定の色の物体の座標を検出するプログラムのサンプルを書いた.青色の物体がどこにあるかを検出できるようにしてあります.
手法は以下のような手順です

  1. HSV色空間に変換して
  2. 青色を部分を抽出して
  3. 抽出した各部分をラベリングして
  4. ラベリングした塊の面積,座標を計算
  5. 最も面積の大きい座標を返す

2. ソースコード


3. 結果

test()でsamle.jpgに対して処理を行います.またmain()でウェブカメラの画像をリアルタイムで処理します.
以下はsample.jpgに処理を適用した結果です.



4.参考




2015年10月25日日曜日

ウラムの螺旋を描く[Processing]

 素数っていいですよね。孤高な感じがして。でもなんか覚えにくいイメージがありますね。私の場合は,数字を覚えるときに,数字が小さいとその数字を因数分解して覚えます。そうすると「ああ,2^2+3^3=108なのか。2,2,3,3で覚えやすいな」みたいなことがあるのですが,素数はそういった覚え方ができないのがまた逆に良いですね。

1.ウラムの螺旋とは


 ウラムの螺旋はウラムさんが見つけたものです。とりあえずこんな感じの図形のこと。なんとなくn斜めに線が浮かび上がっているのがわかるでしょうか。
 この図形は素数に当たる数字の場所を緑色に塗ったものです。素数というと,なんだか法則性があるようでないような感じがしますが,このように図にしてみると法則性がありそうだということがわかります。
 この図形の書き方は,数字をある規則によって並べていき,素数を塗ると幾つかの線が現れるというものです。以下にその数字の並べ方を書いておきます。自分の手で一回書いてみると面白いですよ。
素数を青色に塗ってみました。なんとなく線が浮かび上がってくるのがわかるでしょうか。なぜこのような性質があるのかに対して,証明は確かされていなかったはず。しかしこのウラムの螺旋から,素数を多く生産する数式がたくさん発見された。

 ちなみに二次式で素数を大量生産できる式といえばこんな式がある。
n=0なら41(素数),n=1なら43(素数),n=2なら47(素数),n=3なら53(素数)……と,n=39までは全部素数になる。この数式であらわされる素数をオイラー素数と言ったりする。昔この式であらわされる自然数は全部素数になるよ,というデマに騙されたことがある。n=41を入れれば合成数になることは明らかなのにね。

2.ソースコード


 余談はさておき,プログラムはこんな感じ。素数はエラトステネスの篩で計算した。色を塗る場所の計算がちょっと雑で,上の表とは違っているかもしれない。まぁでも見るだけなら問題ない。
 ホントよくこんなこと考えたよね。まず数字を並べたりしないし,素数を縫ったりもしないからなぁ……。

boolean[] Board;

void setup() {
  /* 適当な大きさでウィンドウを開く */
  size(512, 512);

  /* indexが素数ならfalseの配列を作る */
  Board = new boolean[width*height+1];
  Board[0] = Board[1] = true;

  /* 素数かどうかをあらかじめエラトステネスの篩で計算 */
  for (int i = 2; i < sqrt(Board.length) ; i++ ) {
    /* iが素数ならば */
    if ( Board[i] == false ) {
      /* その倍数は合成数 */
      for (int j = 2 * i; j < Board.length; j += i ) {
        Board[j] = true;
      }
    }
  }
}

void draw() {
  /* 背景を黒にする */
  background(0);

  /* 描画する色を緑にする */
  stroke(0, 255, 0);
  fill(0, 255, 0);
  
  /* 原点を中心に持ってくる */
  translate(width/2, height/2);

  /* 描画する座標 */
  int x = 0;
  int y = 0;
  /* 座標をどう動かすかを計算するための変数 */
  int dx = 1;
  int dy = 0;
  int n = 1;
  
  for (int i = 1; i < width * height + 1; i++ ) {
    /* 素数ならば描く */
    if ( Board[i] == false )
      rect(x, y, 1, 1);  

    /* iに応じて次の描画場所を決める */
    if ( n * n + 1 == i ) {
      dy = (n % 2)*2-1;
      dx = 0;
      n++;
    } else if ( n * n - n + 1 == i) {
      dx = (n % 2)*2-1;
      dy = 0;
    }
    x += dx;
    y += dy;
  }
}

遺伝的アルゴリズムでオセロのAIを学習させてみた[C#]

 最近何かとネットに上がっている遺伝的アルゴリズムを自分でも実装してみた。

 本当は先々週ぐらいにもう完成してたけれどネットが繋がらなくなったりパソコンが起動しなくなったりと,結構面倒なことが起きたので遅れてしまった。まるで何か大きな力が私を止めようとしていたのか……?何事もなければいいのだが。

1.動画

 自分が人工知能や人工生命みたいなワードに興味を持ったきっかけの一つが,むにむに教授さんがやってた遺伝的アルゴリズムシリーズ。今回はそれに便乗してオセロのAIを作ってみて,その評価関数を最適化するのに遺伝的アルゴリズムを使ってみた。そんでもって,せっかくなので動画にしてみた。この動画に刺激されてこういう分野に興味を持ってくれる人が出てくると嬉しい。

 動画編集は初めてで苦戦したし自分でもうまくできていないのはわかっているが,まぁ一度やってみたかったからいいんだ。楽しかったし。

2.解説

 まぁ詳しい内容は割りと動画内でしたつもりだが,もう少しプログラムの詳細をここでは書く。あ,一応述べておくけれど動画内では学習とかなんとか言ってたけどただの最適化だからね。学習って聞くと誤解する人がいそうだから本当は避けるべきだったのかもしれない。

 評価に際して対戦させたわけだが,実際はバックグラウンドで行った。あとAIはすべて1手先しか読んでいない。ミニマックス法で実装したら割りと学習には時間がかかりそうだったしめんどくさいし……。

 交配は64の個体のうち,勝った回数が最も多かった上位10個体から行った。その上位10個体の遺伝子はそのまま引き継ぐ。

 今回は150世代までやったわけだけれども,実際には条件をちょっと変えて行ってみたのだが,大体同じ感じになった。だが対戦相手が変わるため,遺伝子が収束しないというのがちょっと問題かな。あれ?なんか問題多いな。まぁいいか。

3.学習結果

 一応どうなるのかを予想してみた。イメージではこんな感じになると思っていた。

 カドの点数は絶対高いだろうというのは間違いないだろうと思っていた。カドの周りの点数が下がり,自分がそこに置かないようにするというのもまぁ想定の範囲内だった。だが端(一番外側)の部分は,ほとんど高くなると思っていたが実際は違った。

 端(一番外側)の部分でも,中央寄りの部分は点数が下がっている。また中央の部分が下がっている。その辺が驚きだった。いやだって取れるもんはとったほうがいいじゃん?と思っていましたので。

 これについてちょっと考えていたのだが,おそらく「すぐひっくり返されるところ」の点数が下がっていき,「ひっくり返されにくいところ」の得点が高くなっているのだろうと思う。彼らの対戦をすべて保存して状態を解析してみないとわからないのだが,間違いなく相関関係はあるだろうね。

 まぁ実を言うとオセロのAIを作るのは人生で2回目。1回目のプログラムはちゃんとミニマックス法を実装してた。その時の経験から言わせてもらうと,終盤は正直カドを取るということよりも量を取る方に評価関数を変えたほうが強くなった。確か40手ぐらい打ったらもうカドよりも量を取る方に評価関数を変えて,深読みさせたほうが強かった。終盤だとカドを取るべきと取らざるべきがあるからね。そう考えるとこの学習もミニマックス法を実装したほうが良かったかもしれない。

 ソースコードは公開しようと思ったけれど,いじっている間にスパゲティになったので公開はしないつもり。でもオセロはここで書いたやつを使いまわしてる。

2015年10月21日水曜日

黄金比でひまわりの種の配列を描く[Processing]

 現状すべての記事に「プログラミング」タグが付いていることに困惑を隠せない今日この頃です。

 今日は黄金比やフィボナッチ数列の紹介でよく見るこの図形を描きます。


 この図形よく見かけるけれど描き方が今まで全くわからんかったのです。某専門学校に行ってる友人が授業で出されたというソースコードを見て知りました。今日はこれを紹介します。

1.描き方


 描き方は簡単。点を222.5°ずつ回転させながら打っていくだけ。
極座標表示するならこんな感じ。kを整数とし,θ0 = 222.5°としたとき,
という感じになる。この辺りで「黄金比使ってないやんけ!」という突っ込みが来そうだ。だが実は222.5という数字の中に隠されている。ちょっとここで222.5に黄金比φ=1.618...を掛けてみると,

222.5 × 1.618 = 360.005

 360という数字は一周の角度360°のことです。つまり222.5という数字は360÷(黄金比)の値だったのです。 さらにここでラジアン表記してみるとθ0は次のように表されます。
円周率と黄金比のイリュージョンやぁぁぁ!と叫びたくなりますね。非常に美しいですね。それと同時になぜ円周率を2πと定義しなかったのかというあの話を思い出しますね。

2.ソースコード


 Processingではラジアンでsin,cosは定義されているのでラジアンで計算しています。
/* 黄金比 */
float golden_ratio = (sqrt(5)-1.0)/2.0;

void setup() {
  /* 適当な大きさでウィンドウを表示 */
  size(512, 512);
}

void draw() {
  /* 背景を黒くする */
  background(0);
  /* 原点を画面中央へ寄せる */
  translate(width/2, height/2);
  rotate(PI);

  /* 描画の色を緑にする */
  stroke(0, 255, 0);
  fill(0, 255, 0);

  /* 角度は2π/φ[rad]ごとにプロットする */
  /* つまり円一周をを黄金比で割った値 */
  float theta = 2*PI/golden_ratio;
  float x, y;

  for (int i = 0; i < width/2; i++ ) {
    /* iは中心からの距離であり,また点の番号でもある */
    /* 原点からiの距離,角度i*thetaにプロット */
    x = i * cos(i*theta);
    y = i * sin(i*theta);

    ellipse(x, y, 8, 8);
  }
}

2015年10月16日金曜日

オセロのコマンドプロンプト上ゲームを作った[C#]

 まぁGUIで作ったのだが,それを入れるとかなり長くなってしまうからコマンドプロンプト上のものにした。近々GUIで動くやつを公開する予定なのでその時は参考にしてください。

1.簡単な説明

プログラムを見てもらえればGUIに移行するのは容易だと思います。一手先しか読まないAIも実装してあります。
 このAIは評価ボードと呼ばれるものにもとづいて自分の手を決めています。どの場所がどのくらい価値が有るかを示したものです。
この評価ボードが一番高くなる手を打つというわけです。しかしそのため,2手先で全部ひっくり返されてしまうような手を打ってくるので,もしそれを直したいときはミニマックス法やアルファベータ法などを用いるといいと思います。

2.プログラム

まぁ長々と話してもプログラムの説明は難しいのでとりあえずコードです。コメントは多めにしてあるからそれを見てください。
using System;
using System.Collections.Generic;

class Program
{
    static void Main()
    {
        // ゲーム用クラスを用意
        var r = new Reversi();
        int x,y;
        
        while(r.CheckFinish() == false )
        {
            Show(r);
            
            if( r.Turn )
            {
                // 人間からの入力
                Console.WriteLine("左からの番号を入力");
                if ( int.TryParse(Console.ReadLine(),out x) == false )continue;
                Console.WriteLine("上からの番号を入力");
                if ( int.TryParse(Console.ReadLine(),out y) == false )continue;
                r.PutStone(x,y);
            }
            else
            {
                // AIに打たせる
                r.AIPut();
            }
        }
    }
    
    static void Show(Reversi r)
    {
        Console.Clear();
        
        Console.Write(" ");
        for(int i = 0 ; i < Reversi.N ; i++ )
            Console.Write(" {0}",i);
        Console.WriteLine();
        
        for(int i = 0 ; i < Reversi.N ; i++ )
        {
            Console.Write(i);
            for(int j = 0 ; j < Reversi.N ; j++ )
            {
                switch( r.Board[j,i] )
                {
                    case Reversi.BLACK:
                        Console.Write("●");
                        break;
                    case Reversi.WHITE:
                        Console.Write("○");
                        break;
                    default:
                        Console.Write(" ");
                        break;
                }
            }
            Console.WriteLine();
        }
    }
}

class Reversi
{
    /// <summary>
    /// 状態を保存するボード
    /// </summary>
    public int[,] Board;
    /// <summary>
    /// 一辺のマスの数
    /// </summary>
    public const int N = 8;
    /// <summary>
    /// 何もない状態
    /// </summary>
    public const int NONE = 0;
    /// <summary>
    /// 白の石
    /// </summary>
    public const int WHITE = 1;
    /// <summary>
    /// 黒の石
    /// </summary>
    public const int BLACK = -1;
    /// <summary>
    /// どちらの順番がを示す変数(trueなら黒)
    /// </summary>
    public bool Turn;
    /// <summary>
    /// ボードの状態を保存する変数
    /// </summary>
    private List<int[,]> BoardHistory;
    /// <summary>
    /// 手番を保存する変数
    /// </summary>
    private List<bool> TurnHistory;
    /// <summary>
    /// デフォルトの評価ボード
    /// </summary>
    private int[,] EvaluationBoard = new int[,]{
                {  60, -25, 15, 15, 15, 15,-25, 60},
                { -25, -50,-30,-30,-30,-30,-50,-25},
                {  15, -30, 15, 15, 15, 15,-30, 15},
                {  15, -30, 15, 25, 25, 15,-30, 15},
                {  15, -30, 15, 25, 25, 15,-30, 15},
                {  15, -30, 15, 15, 15, 15,-30, 15},
                { -25, -50,-30,-30,-30,-30,-50,-25},
                {  60, -25, 15, 15, 15, 15,-25, 60}
    };

    /// <summary>
    /// コンストラクタ
    /// </summary>
    public Reversi()
    {
        this.Init();
    }

    /// <summary>
    /// 初期化メソッド
    /// </summary>
    public void Init()
    {
        this.Board = new int[N, N];
        this.Board[3, 3] = WHITE;
        this.Board[4, 4] = WHITE;
        this.Board[3, 4] = BLACK;
        this.Board[4, 3] = BLACK;
        this.Turn = true;

        this.BoardHistory = new List<int[,]>();
        this.TurnHistory = new List<bool>();
    }

    /// <summary>
    /// 置ける場所かどうかを判定するメソッド
    /// </summary>
    /// <param name="x">判定するx座標</param>
    /// <param name="y">判定するy座標</param>
    /// <returns>置ける場合はtrue</returns>
    public bool CanPut(int x, int y)
    {
        //とりあえず実際に置いてみる
        var ret = Put(x, y);

        //おけなかったらおけない場所(当然)
        if (ret == false)
            return false;

        //勝手に置くのはダメなので元に戻す
        Undo();

        return true;
    }

    /// <summary>
    /// 元に戻すメソッド
    /// </summary>
    private void Undo()
    {
        // 一番最後の要素のindex
        int n = this.BoardHistory.Count - 1;

        if (n < 0)
            return;

        // 一個前の状態に戻す
        this.Board = this.BoardHistory[n];
        this.Turn = this.TurnHistory[n];

        // その時のボードの状態・手番は消す
        this.BoardHistory.RemoveAt(n);
        this.TurnHistory.RemoveAt(n);
    }

    /// <summary>
    /// 手番を変更するメソッド
    /// </summary>
    private void ChangeTurn()
    {
        // とりあえず順番変える
        this.Turn = !this.Turn;

        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < N; j++)
            {
                // おける場所が一か所でもあればOK
                if (CanPut(i, j) == true)
                    return;
            }
        }

        // おける場所がなかったので手番は元に戻る
        this.Turn = !this.Turn;
    }

    /// <summary>
    /// ゲームが終了したかどうかを判定するメソッド
    /// </summary>
    /// <returns>終了したらtrue</returns>
    public bool CheckFinish()
    {
        // ChangeTurnで!Turnの場合は置ける場所がないのはわかっている
        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < N; j++)
            {
                if (CanPut(i, j) == true)
                    return false;
            }
        }

        // 今の人も置く場所がないので終わり
        return true;
    }

    /// <summary>
    /// 石の数を数え上げるメソッド
    /// </summary>
    /// <param name="target">数える対象</param>
    /// <returns>石の数</returns>
    public int CountStone(int target)
    {
        int count = 0;
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                if (Board[i, j] == target)
                    count++;

        return count;
    }

    /// <summary>
    /// 石を置き,手番を変更するメソッド
    /// </summary>
    /// <param name="x">置く場所のx座標</param>
    /// <param name="y">置く場所のy座標</param>
    /// <returns>置くことが出来たらtrueを返す</returns>
    public bool PutStone(int x, int y)
    {
        // とりあえず置く
        var flag = Put(x, y);

        if (flag == false)
            return false;

        ChangeTurn();

        return true;
    }

    /// <summary>
    /// 石を置くメソッド
    /// </summary>
    /// <param name="x">置く場所のx座標</param>
    /// <param name="y">置く場所のy座標</param>
    /// <returns>置くことが出来たらtrueを返す</returns>
    private bool Put(int x, int y)
    {
        //範囲外なら何もせず返す
        if (InRange(x, y) == false)
            return false;
        // なにかあったら置けない
        if (Board[x, y] != NONE)
            return false;

        // ひっくり返したかどうかを格納するメソッド
        bool isChanged = false;
        // 現在の状態を一旦保存する
        var currentBoard = (int[,])(this.Board.Clone());
        // 現在の攻撃側はどちらかを一旦保存する
        var currentTurn = this.Turn;

        for (int i = 0; i < 9; i++)
        {
            //これでdx,dyは-1から1までの値が入る
            int dx = i / 3 - 1;
            int dy = i % 3 - 1;

            //両方共0じゃなければ
            //(dx,dy)方向へひっくり返せるかを調べる
            if (dx != 0 || dy != 0)
                isChanged = isChanged | Reverse(x, y, dx, dy);
        }

        // ひっくり返さなかった場合はfalse
        if (isChanged == false)
            return false;


        // ココに来てやっと置くことが出来る
        this.Board[x, y] = NowStone();

        // 手番とボードの状態を保存
        this.BoardHistory.Add(currentBoard);
        this.TurnHistory.Add(currentTurn);

        return true;
    }

    /// <summary>
    /// 石をひっくり返すメソッド
    /// </summary>
    /// <param name="x">石をおいたx座標</param>
    /// <param name="y">石をおいたy座標</param>
    /// <param name="dx">調べる方向のx</param>
    /// <param name="dy">調べる方向のy</param>
    /// <returns>ひっくり返せたらtrue</returns>
    private bool Reverse(int x, int y, int dx, int dy)
    {
        var attack = NowStone();
        var defense = -attack;

        // その方向が枠の外ならひっくり返せない
        if (InRange(x + dx, y + dy) == false)
            return false;
        // 一個先を見て敵の石じゃなかったらひっくり返せない
        if (Board[x + dx, y + dy] != defense)
            return false;


        // その先を見ていく
        for (int i = 2; i < N; i++)
        {
            int index_x = x + i * dx;
            int index_y = y + i * dy;

            if (InRange(index_x, index_y) == false)
            {
                //範囲外ならfalse
                return false;
            }
            else if (Board[index_x, index_y] == attack)
            {
                //探した先に攻撃側の駒があった場合はひっくり返す
                for (; i >= 1; i--)
                    Board[x + i * dx, y + i * dy] = attack;
                return true;
            }
            else if (Board[index_x, index_y] == NONE)
            {
                // その先に仲間の石がなかったのでfalse
                return false;
            }
        }

        // 仲間の石が見つからなかったのでfalse
        return false;
    }

    /// <summary>
    /// 現在攻撃側の石を置くメソッド
    /// </summary>
    /// <returns></returns>
    private int NowStone()
    {
        if (this.Turn)
            return BLACK;
        else
            return WHITE;
    }

    /// <summary>
    /// 特定の位置が配列の範囲内かどうかを判定するメソッド
    /// </summary>
    /// <param name="x">判定するx座標</param>
    /// <param name="y">判定するy座標</param>
    /// <returns>範囲内ならtrue</returns>
    private bool InRange(int x, int y)
    {
        if (x < 0 || x >= N)
            return false;
        if (y < 0 || y >= N)
            return false;

        return true;
    }

    /// <summary>
    /// 評価関数
    /// </summary>
    /// <param name="evaluationBoard">評価ボード</param>
    /// <returns>評価値(黒が有利なら正)</returns>
    private int Evaluate(int[,] evaluationBoard)
    {
        var point = 0;

        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < N; j++)
            {
                point += evaluationBoard[i, j] * Board[i, j];
            }
        }

        return point;
    }

    /// <summary>
    /// 人工知能による最適な手の選択し,置く
    /// </summary>
    /// <param name="evaluationBoard">評価ボード</param>
    public void AIPut(int[,] evaluationBoard)
    {
        int max = int.MinValue;
        int x = 0;
        int y = 0;
        // 自分の石を覚えとく
        int myStone = NowStone();

        // 終わっていたら関係ない
        if (CheckFinish() == true)
            return;

        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < N; j++)
            {
                // とりあえず石を置く
                if (PutStone(i, j) == true)
                {
                    // 置けたら盤面評価
                    // 自分の石の値をかけて常に正にする
                    var point = myStone * Evaluate(evaluationBoard);

                    // ポイントが高かったら
                    if (point > max)
                    {
                        // その手を保存する
                        x = i;
                        y = j;
                        max = point;
                    }

                    // 元に戻す
                    Undo();
                }
            }
        }

        // 置ける場所なら置く
        if (InRange(x, y) == true)
            PutStone(x, y);
    }
    
    /// <summary>
    /// 人工知能による最適な手の選択し,置く
    /// </summary>
    public void AIPut()
    {
        AIPut(this.EvaluationBoard);
    }
}

2015年10月15日木曜日

ニュートン法で黄金比を求める[Python]

 まだ慣れないPythonで書いたのでミスや書き方がPythonっぽくないかもしれないが許してつかあさい。

1.ニュートン法とは

平たく言えば漸化式にもとづいて計算していくと答えが出てくるというもので,細かい説明をはしょると次の式を用います。
X0にf(x) = 0の解に近い値を入れて計算をしていくと,f(x) = 0 の解に収束するという素晴らしい漸化式です(関数によって必ずしも収束するとは言い切れない)。

 例えばf(x) = x^2 - 2 を代入した次の式は有名。
大学入試の相加平均・相乗平均とかの問題に出てきそうな形してますね。X0 = 1.5 を代入すると,X1=1.41666……,X2=1.41421568……と,どんどんx^2 - 2 = 0の解である√2に近づいていきます。

2.黄金比とは

いろいろな書き方がありますが,黄金比とは,関数f(x)
の正の解とします。つまり,黄金比φは
と表せる事ができるわけです。黄金比については面白い性質が沢山あったりするわけですが,プログラムを書いて,さらにその話を書くにはちょっと余白が狭すぎます。

3.求める方法

1,2で説明したことをくっつければいいわけですから,1のf(x)に,2で示したf(x)を代入します。すると
という漸化式が得られます。φ=1.6180339……ですから,X0 = 1.5ぐらいで始めれば答えは出そうです。

 まぁちなみに手元の電卓を使ってX0=1.5を代入すると,X1 = 1.625,X2 = 1.618055……となっていますので正しそうです。

 ここで電卓でやればいいじゃん,という声が出てきそうですね。ですが今回は任意の桁数を計算したいのです。1000桁とか表示してくれる電卓があればいいのですがそうも行きません。

 なので一工夫を加えます。まずXnを次のように置きます。
ここでan,bnは整数です。なぜこんなことをしたのかというと,コンピュータは小数を扱うのが苦手です。なので小数を整数の組み合わせ,つまり分数で扱おうという作戦です。
 これを代入すると次のようになります。

4.ソースコード

前置きが長かったです。つまり言いたいことはこういうこと。本当にプログラミング言語は便利な言葉やなぁ。
# coding: utf-8
# 2015/10/15 黄金比を求めるプログラム

# python2.xは下の一行がいる
from __future__ import print_function
import math

if __name__ == "__main__":
   # 分母分子を表す変数を用意
   b0 = 3
   a0 = 2
   bn = 0
   an = 0

   # 桁数をN=1000とする
   N = 1000

   # 誤差は2次収束する
   # つまり正しい桁数が1回の計算で2倍になっていく
   # よってN桁正しい値を得るためにはlog[2]Nとなる
   # ここでは多めにlog[1.5]Nにした
   for i in range(0,int(math.log(N,1.5))):
      an = a0*(2*b0-a0)
      bn = a0*a0+b0*b0
      
      a0 = an
      b0 = bn
   
   # 整数部は消す
   bn = (bn - an)*10
   print("1.")
   
   # 一桁ずつ表示する
   for j in range(N):
      print(int(bn/an),end="")
      if( j % 10  == 9):
         print(" ",end="")
      if( j % 50 == 49):
         print()
      bn = (bn % an)*10
   

 実行結果

2015年10月13日火曜日

対数(log)を手計算で求める方法とプログラムで求める方法

 大学入試では頻出の対数。その数値を計算させることはないですが,知っておくと面白いと思うのでまとめておきます。
 余談ですが無理数の計算って何かワクワクしますよね。未知の数字を求めてる感じが凄く。

1.手計算で求める(その1)

 一番簡単だと思う方法です。ただしあまり精度は良くないです。

 まず2の10乗を計算してみましょう。
 2の10乗までは大学入試で覚えておいたほうがいいと思います。でも計算すればすぐ出てきます。ここでこんな近似をしてします。
 「≒」は大体等しいという意味です。1024は大体1000だろう,と大まかに決めてしまうのです。ここで両辺を底10の対数を取ります。

 するとlog[10]2≒0.3という結果が出てきます。電卓で計算すると0.3010…なので小数点以下2桁が正しいことがわかります。

 では他の場合はどうでしょう。3について考えてみます。さっきと同じように累乗していきます。
 3^4=81はおおよそ80,つまり80=2^3✕10と表せます(3^2=9≒10じゃないかと思う人もいるかもしれませんが累乗した数は大きい方が精度が上がります)。両辺の対数をとると,
 となって先ほど求めたlog[10]2 ≒ 0.3を代入すると
 となります。log[10]3 = 0.47712なので小数点以下2桁が正しいことがわかります。
 
 同様にして4,5,6,7,8,9,10,11,12…と求めていくことができますが,大きい数字になるに連れて誤差が大きくなっていくので気をつけてください。
 基本的には素数さえ求めれば合成数は素因数の対数の和なので簡単に求められます。以下に5,7,11,13の場合を示します。

 5の場合は10=5✕2ということを利用します(log[10]5=0.6989…)。
 7の場合は7^2 = 49≒50を利用します(log[10]7=0.8450…)。

 11の場合は11✕9 = 99 ≒ 100 を利用します。累乗してないので誤差は他に比べて大きいと思います(log[10]11=1.0413…)。

 13の場合は13^3 = 2197 ≒ 2200 を利用します。この辺になるとどこで近似するかは好みですかね(log[10]13=1.1139…)。
 極めると面白いものが見えてきそうですね。

2.手計算で求める(その2)

 この方法はあまりおすすめはしません……が,こういう方法もあるんだと思ってみてってください。

 対数は次の式によって表すことができます。ここで底はe(自然対数)です。
 この式はf(x) = a^xを定義によって微分したものと,両辺の対数をとって微分した結果を照らし合わせると見えてきますが,循環論法になりそうなのでこれが成り立つものとして見て下さい。

 ここで対数の重要な公式を使います。
 (2)にa=10,c=e(自然対数),b=aを代入すると(1)より次の式が求まります。
 a = 2について考えてみようと思うのですが,残念ながら2^h,10^hを求めるのがなかなか面倒なのです。h=0.5とすると,ルート2,ルート10の値を求めればよいのですが,
 となり,正直うまくいきません。平方根ならばニュートン法や開平法で求められることには求められるのですが,正直あまり合理的ではなさそうですね。ちなみにh=1/16の場合で0.2860...です。

3.積分を使った方法

 こっちは手計算で求めるというにはちょっと違うけれど,求められる方法。次の積分を用いる。
 手計算で求めるならばこれを区分求積法で求めればオッケーです。でもごめんなさい。数学好きの高校生の子が見ていたら申し訳ないけれどここからプログラミングの話にシフトチェンジしていくのです。

 上の式は積分ですが,積分を数値的に求める方法としてモンテカルロ法や台形公式なんかがまぁ楽でいいですね。それでちょっと実装してみました。ここでは自然対数で常用対数ではないのでそこんところよろしくお願いします。

/**
 * 2015/10/13
 * 対数を求めるプログラム
 */
using System;

class Program
{
      // 試行回数
      const int N = 10000;
      
      static void Main()
      {
            for(int i = 2 ; i < 100 ; i++ )
            {
                  Console.Write("log({0})\t",i);
                  Console.Write("{0:f5}\t",Math.Log(i));
                  Console.Write("{0:f5}\t",MyLog(i));
                  Console.Write("{0:f5}\t",MyLog2(i));
                  Console.WriteLine();
            }
      }
      
      // モンテカルロ法で対数を求めるプログラム
      static double MyLog(double d)
      {
            // 乱数生成用クラス
            var rand = new Random();
            var count = 0;
            
            // 1未満の場合は定義しない 
            if( d < 1 )
                  return -1;
            
            // 打点をする
            for(int i = 0 ; i < N ; i++ )
            {
                  // x = [1,d)の乱数生成
                  var x = rand.NextDouble() * ( d - 1 ) + 1;
                  // y = [0,1)の乱数生成
                  var y = rand.NextDouble();
                  
                  // グラフの下側ならtrue
                  if( y < 1.0/x )
                        count++;
            }
            
            return (d-1.0) * count / N;
      }
      
      // 台形公式で対数を求めるプログラム
      static double MyLog2(double d)
      {
            double integral = 0.0;
            double dx = ( d - 1.0 )/N;
      
            // 1未満の場合は定義しない 
            if( d < 1 )
                  return -1;
            
            // 積分していく
            for(int i = 1 ; i < N - 1 ; i++ )
                  integral += 1.0 / ( i * dx + 1);
            
            // y0/2,yn/2を加える
            integral += 1.0 / 2.0;
            integral += 1.0 / d;
            
            return integral * dx;
      }
}
 やっぱりモンテカルロ法には限界がありますね。

2015年10月12日月曜日

ナイトツアーの数え上げプログラム[C#]

 今日は作ろうと思ってたプログラムが完成しなかった。残念。せっかくの休日だったのになぁ……。まぁいいや。今日はナイトツアーの数え上げプログラムです。

1.ナイトツアーとは

ナイトツアーというのは,チェスのナイトの駒でチェスの盤をすべて一度だけ通るというパズルゲームです。レイトン教授のゲームの中にこのパズルゲームがあって,小学生の頃友達と競争していました。

 イメージがわかない人はWikipediaとか見たり,「ナイトツアー ゲーム」で検索すれば出てきますのでやってみると面白いと思います。
 
 チェスのナイトの動きはこんな感じ。「馬」がナイトです。丸のある場所に移動できます。

 注目すべきはこのマスの色。チェスの盤は色が付いているのですが,ナイトは2つの色を交互に移動します(なので盤の色が黒と白ならナイトは黒→白→黒→白と移動します)。両方の色に行けるということは,すべての場所に移動できる可能性があるということです。

 ちなみにビショップ(将棋で言う角)は同じ色の上にしかいけません。なのでビショップでこの問題を解くことができないのは容易に証明できます。
 
 今日作ったプログラムはナイトツアーを解くプログラムです。下は5✕5のあるパターンです。数字の順番にナイトが移動すると漏れ無く盤の上を移動できます。
 ちなみにコツは外側から攻めてくということです。上の答えも子供の頃にやったまんまの答えです。案外覚えているもんですねw。

 プログラムが正しければ5✕5の場合は304通りあるはずです。6✕6や7✕7は時間がかかりそうだったのでやりませんでしたが,もし実行して答えが出たという人は教えて下さい。

2.プログラム

 再帰関数を使っているからちょっと不安。あまり大きい数字を入れるとスタックオーバーフローとか出てくるかもしれないのでそこはご勘弁を。

/**
 * 2015/10/12
 * ナイトツアーの数え上げプログラム
 */
using System;
using System.Collections.Generic;

class Program
{
    static void Main()
    {
        //5x5のパターンを数え上げる
        var nt = new NightTour(5);

        Console.WriteLine("Search Start...");

        //すべて数え上げる
        nt.Start(true);

        //一個見つければいい人はこう
        //nt.Start(false);
        //nt.Show();

        Console.WriteLine("Search End...");
    }
}

class NightTour
{
    /// <summary>
    /// 枠の大きさ
    /// </summary>
    int N = 8;
    /// <summary>
    /// ボード
    /// </summary>
    int[,] Board;
    /// <summary>
    /// ナイトのX座標
    /// </summary>
    int X;
    /// <summary>
    /// ナイトのY座標
    /// </summary>
    int Y;
    /// <summary>
    /// 答えの個数
    /// </summary>
    int Count;

    /// <summary>
    /// コンストラクタ
    /// </summary>
    /// <param name="N">枠の大きさ</param>
    public NightTour(int N)
    {
        this.N = N;
        Board = new int[N, N];
        X = 0;
        Y = 0;

        // 左上を最初の場所にする
        Board[0, 0] = 1;
    }

    /// <summary>
    /// 旅を開始するメソッド
    /// </summary>
    /// <param name="forever">永遠に続けるかどうか</param>
    public void Start(bool forever)
    {
        Count = 0;
        Tour(2,forever);
    }

    /// <summary>
    /// 答えを表示するメソッド
    /// </summary>
    public void Show()
    {
        Console.WriteLine("Count = {0}",Count);
        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < N; j++)
            {
                Console.Write("{0,5}", Board[i, j]);
            }
            Console.WriteLine();
        }
        Console.WriteLine();
    }

    /// <summary>
    /// ナイトツアーを時効するメソッド
    /// </summary>
    /// <param name="depth">深さ</param>
    /// <param name="forever">永遠に続けるかどうか</param>
    /// <returns>旅が終わったかどうか</returns>
    private bool Tour(int depth,bool forever)
    {
        // すべて塗りつくした場合
        if (depth > N * N)
        {
            Count++;
            // 永遠に続けるならばここで表示する
            if (forever)
                Show();

            return true;
        }

        // 移動先を探す
        var nextPlace = SearchNextPlace();

        if (nextPlace.Length == 0)
            return false;

        foreach (Tuple<int, int> t in nextPlace)
        {
            // 現在の状態を保存しておく
            int temp_x = X;
            int temp_y = Y;

            // 移動する
            X = t.Item1;
            Y = t.Item2;
            Board[X, Y] = depth;

            //移動先で行動させる
            if (forever)
                Tour(depth + 1,forever);
            else
                if (Tour(depth + 1, forever) == true)
                    return true;

            // 元に戻す
            Board[X, Y] = 0;
            X = temp_x;
            Y = temp_y;
        }

        return false;
    }

    /// <summary>
    /// 移動できる座標を探すメソッド
    /// </summary>
    /// <returns>移動可能な位置座標</returns>
    private Tuple<int, int>[] SearchNextPlace()
    {
        var ret = new List<Tuple<int, int>>();

        // マスの大きさを無視してとりあえず移動先を考える

        ret.Add(Tuple.Create(X - 1, Y - 2));
        ret.Add(Tuple.Create(X - 1, Y + 2));
        ret.Add(Tuple.Create(X + 1, Y - 2));
        ret.Add(Tuple.Create(X + 1, Y + 2));
        ret.Add(Tuple.Create(X - 2, Y - 1));
        ret.Add(Tuple.Create(X - 2, Y + 1));
        ret.Add(Tuple.Create(X + 2, Y - 1));
        ret.Add(Tuple.Create(X + 2, Y + 1));

        // 移動できない場所を削除する
        for (int i = 0; i < ret.Count; i++)
        {
            if (CanMove(ret[i].Item1, ret[i].Item2) == false)
            {
                ret.RemoveAt(i);
                i--;
            }
        }

        // 配列にして返す
        return ret.ToArray();
    }

    /// <summary>
    /// 動ける場所かどうかを判定するメソッド
    /// </summary>
    /// <param name="x">x方向のindex</param>
    /// <param name="y">y方向のindex</param>
    /// <returns>動ける場所ならばtrue</returns>
    private bool CanMove(int x, int y)
    {
        if (InRange(x, y) == false)
            return false;

        if (Board[x, y] != 0)
            return false;

        return true;
    }

    /// <summary>
    /// インデックスが配列の範囲内かどうかを判定するメソッド
    /// </summary>
    /// <param name="x">x方向のindex</param>
    /// <param name="y">y方向のindex</param>
    /// <returns>範囲内ならtrue</returns>
    private bool InRange(int x, int y)
    {
        if (x < 0 || x >= N || y < 0 || y >= N)
            return false;
        return true;
    }
}
 6✕6とかやった人いるのかなぁ。魔法陣の数え上げは5✕5までだった気がする。こういう数え上げとかは競技プログラミングっぽくてもっといい方法がありそうな気がする。

2015年10月11日日曜日

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

 こっちのプログラムに,こっちで実装したFractionクラスを適用してみる。と言っても一昨日のやつのdoubleをFractionに変えてちょこっと調整しただけだけれども。

 ソースコードの前に実行結果を出しておきます。
 一乗和の公式:(1/2)*n^1+(1/2)*n^2
 二乗和の公式:(1/6)*n^1+(1/2)*n^2+(1/3)*n^3
 三乗和の公式:(0/1)*n^1+(1/4)*n^2+(1/2)*n^3+(1/4)*n^4
 四乗和の公式:(-1/30)*n^1+(0/1)*n^2+(1/3)*n^3+(1/2)*n^4+(1/5)*n^5
 五乗和の公式:(0/1)*n^1+(-1/12)*n^2+(0/1)*n^3+(5/12)*n^4+(1/2)*n^5+(1/6)*n^6
 六乗和の公式:(1/42)*n^1+(0/1)*n^2+(-1/6)*n^3+(0/1)*n^4+(1/2)*n^5+(1/2)*n^6+(1/7)*n^7
 七乗和の公式:(0/1)*n^1+(1/12)*n^2+(0/1)*n^3+(-7/24)*n^4+(0/1)*n^5+(7/12)*n^6+(1/2)*n^7+(1/8)*n^8
 八乗和の公式:出力はされるものの,間違っている。オーバーフローしてるっぽい。
 九乗和の公式:例外が発生する

 本当は100乗和の公式とか出してみたかったんだけれどなぁ。残念。int型は2^31-1までしか扱えないから,おそらくそのせい。BigInteger使うのもありかもしれないけれど,それをやるのはまた今度かな。

 ソースコード。昨日のFractionクラスは省略しています。必ず追加してくださいね。
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(Fraction[,] formula, int N)
    {
        Fraction 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(Fraction[,] d, int N)
    {
        Console.WriteLine("{0}乗和の公式は……", N - 1);
        for (int i = 0; i < N; i++)
        {
            Console.Write("({0})*n^{1}", d[i, N], i + 1);
            if (i + 1 < N)
                Console.Write("+");
        }
        Console.WriteLine();
    }

    // 累乗和の公式を求める連立方程式の拡大係数行列を作るメソッド
    static Fraction[,] MakeFormula(int N)
    {
        var ret = new Fraction[N, N + 1];
        var sum = 0;

        for (int i = 0; i < N; i++)
        {
            // 係数部分を作る
            for (int j = 0; j < N; j++)
                ret[i, j] = (int)(Math.Pow(i + 1, j + 1));

            // 答えの部分を作る
            sum = sum + (int)(Math.Pow(i + 1, N - 1));
            ret[i, N] = sum;
        }

        return ret;
    }
}

 ソースコードほぼ同じじゃねーか。ページ数稼ぎかよ,とは言わないで。

追記:
 255乗和の公式のような,任意のべき乗和を求めるプログラムが完成しました→BigIntegerを使って累乗和の公式を計算[C#]

分数を扱うFractionクラスの実装[C#]

 昨日作った累乗和の公式を求めるプログラムでは,分数を扱えなかったのでdouble型で実装した。今日は分数のクラスを定義したのでソースコードを貼っておく。

 基本的なオペレーターは定義したけれども足りないとか間違っているとかあるかもしれないので指摘or温かい目でご覧ください。

 だけどソースコード長い。GitHubとか始めたほうがいいかもしれんな……。
using System;

class Fraction
{
    /// <summary>
    /// 分子
    /// </summary>
    public int Upper
    {
        get { return _Upper; }
        set { this._Upper = value; CheckFraction(); }
    }
    /// <summary>
    /// 分母
    /// </summary>
    public int Lower
    {
        get { return _Lower; }
        set { this._Lower = value; CheckFraction(); }
    }
    /// <summary>
    /// 分子の値を格納する変数
    /// </summary>
    private int _Upper;
    /// <summary>
    /// 分母の値を格納する変数
    /// </summary>
    private int _Lower;

    public Fraction(int Upper, int Lower)
    {
        Assignment(Upper, Lower);
    }
    public Fraction(int n)
    {
        Assignment(n, 1);
    }
    public Fraction(Fraction f)
    {
        Assignment(f.Upper, f.Lower);
    }

    /// <summary>
    /// double型に変換するメソッド
    /// </summary>
    /// <returns>double型に変換した値</returns>
    public double ToDouble()
    {
        return 1.0 * Upper / Lower;
    }

    /// <summary>
    /// float型に変換するメソッド
    /// </summary>
    /// <returns>float型に変換した値</returns>
    public float ToFloat()
    {
        return 1.0f * Upper / Lower;
    }

    /// <summary>
    /// 約分や分母が0などのミスが無いか確認するメソッド
    /// </summary>
    private void CheckFraction()
    {
        int u = _Upper;
        int l = _Lower;
        int temp;
        bool sign;


        if (l == 0)
        {
            // 分母が0なのでゼロ除算の例外をぶん投げる
            throw new DivideByZeroException();
        }
        else if (u == 0)
        {
            // 分母が0ならばこうする
            _Upper = 0;
            _Lower = 1;
        }
        else
        {
            // 分数の値が負ならsign = true となる
            sign = l * u < 0;
            // とりあえずl,uを正にする(余りの計算をするため)
            l = Math.Abs(l);
            u = Math.Abs(u);

            // 常にu <= lとなるようにする
            if (u > l)
            {
                temp = u;
                u = l;
                l = temp;
            }

            // ユークリッドの互除法
            while (l % u != 0)
            {
                temp = l % u;
                l = u;
                u = temp;
            }

            // 約分
            _Upper = Math.Abs(_Upper)/u;
            _Lower = Math.Abs(_Lower)/u;

            // 正負を直す
            if (sign)
                _Upper = -_Upper;
        }
    }

    /// <summary>
    /// クラス内で代入を行うメソッド(ただしCheckFractonから呼んではならない)
    /// </summary>
    /// <param name="Upper">分母</param>
    /// <param name="Lower">分子</param>
    private void Assignment(int Upper, int Lower)
    {
        this._Upper = Upper;
        this._Lower = Lower;
        CheckFraction();
    }

    public static Fraction operator +(Fraction left, Fraction right)
    {
        int l = left.Lower * right.Lower;
        int u = left.Upper * right.Lower + right.Upper * left.Lower;

        return new Fraction(u, l);
    }

    public static Fraction operator -(Fraction left, Fraction right)
    {
        right.Upper = -right.Upper;

        return left + right;
    }

    public static Fraction operator *(Fraction left, Fraction right)
    {
        int l = left.Lower * right.Lower;
        int u = left.Upper * right.Upper;

        return new Fraction(u, l);
    }

    public static Fraction operator /(Fraction left, Fraction right)
    {
        // 逆数にする(0除算ならここで例外が発生する)
        right = new Fraction(right.Lower, right.Upper);
        return right * left;
    }
    public static Fraction operator /(Fraction f, int n)
    {
        return new Fraction(f.Upper, f.Lower * n);
    }

    public static Fraction operator *(Fraction f, int n)
    {
        return new Fraction(f.Upper * n, f.Lower);
    }
    
    public static implicit operator Fraction(int n)
    {
        return new Fraction(n);
    }

    public static bool operator <(Fraction left, Fraction right)
    {
        return left.ToDouble() < right.ToDouble();
    }

    public static bool operator >(Fraction left, Fraction right)
    {
        return left.ToDouble() > right.ToDouble();
    }

    public static bool operator <=(Fraction left, Fraction right)
    {
        return left.ToDouble() <= right.ToDouble();
    }

    public static bool operator >=(Fraction left, Fraction right)
    {
        return left.ToDouble() >= right.ToDouble();
    }

    public static implicit operator double (Fraction f)
    {
        return f.ToDouble();
    }

    public static implicit operator float (Fraction f)
    {
        return f.ToFloat();
    }

    public static bool operator ==(Fraction left, Fraction right)
    {
        if (left.Upper == null || right.Upper == null)
            return false;
        return left.Upper == right.Upper && left.Lower == right.Lower;
    }

    public static bool operator !=(Fraction left, Fraction right)
    {
        return !(left == right);
    }

    public override bool Equals(object obj)
    {
        if (obj == null || this.GetType() != obj.GetType())
            return false;

        return (this == (Fraction)obj);
    }

    public override int GetHashCode()
    {
        return this.Upper * this.Lower;
    }

    public override string ToString()
    {
        return String.Format("{0}/{1}", this.Upper, this.Lower);
    }
}

 サンプルコードはまた今度。

追記:
このプログラムを使って累乗和の公式を導出しました→累乗和の公式をプログラムで導出してみる(分数編)

2015年10月10日土曜日

累乗和の公式をプログラムで導出してみる[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#]