2015年03月01日

C#で立方根を求める(ニュートン法による近似解)

   このエントリーをはてなブックマークに追加 Clip to Evernote
どう書く?orgに感謝を込めて」シリーズ その73

■問題 (出題者:にしお さん)
xは0以上1000未満の実数です。 y * y * y = xになるような実数y(立方根)を小数点以下12桁以上の正確さで 求める関数cube_rootを作って下さい。
ただし、このお題の趣旨は実数区間での探索なので、 立方根関数があっても使ってはいけません。 指数関数と対数関数も禁止します。
Pythonで表現した入出力の例:
>>> cube_root(10.0)
2.1544346900318834
>>> _ ** 3
9.9999999999999947
>>> cube_root(100.0)
4.6415888336127793
>>> _ ** 3
100.00000000000003


まずは、立方根の前にイメージしやすい平方根について考ます。

※ ここでは、 x^n は、x の n乗を意味するものとします。

√aを求めることを考えて見ます。(√は、平方根の記号のつもり)
√aは、x^2 = a であるxの値ですので、
x^2 - a = 0 の方程式を解けば答えを求めることでができます。
y = f(x) = x^2 - a のグラフで考えてみると、y座標が0 である xの座標を求めればよいことになります。

では、どのようにして求めるかですが、コンピュータを使う場合は、ニュートン法を利用する方法が知られています。 こちら

考え方はこうです。
まず、解に近い値を x0 としたときに、座標 (x0,f(x0)) との接線を引き、その接線とx軸との交点をx1 とする。
同様に、この座標 (x1, f(x1))の接線を引き、x2を求める、続いて、(x2, f(x2))の接線を引き、x3を求める。
これを繰り返すことで, xnは、解に近づいてゆく、というものです。
xn と そのひとつ前の x の値の差が十分小さくなれば、解の近似値とみなせます。

図が無いので分かりにくいかもしれませんが、ここにある図3を見ていただければ、考え方は理解できると思います。


n 番目の x を xn 、n+1番目の x を xm とすれば、(x(n+1) だと紛らわしいので)

xm = xn - w

と表せます。このときの wは、xn の接線の傾きを g(xn) とすると、

w = f(xn) / g(xn)

です。

※ 傾き = 高さ/ 横幅 ですから、横幅 = 高さ / 傾き となります。

つまり、

xm = xn - (f(xn) / g(xn))

が導き出せます。

次に、グラフ f(x) の接線の傾き、g(xn) を求めます。

まず、y = f(x) のグラフの以下の2つの点を通る直線の傾きを考えてみます。
1: (x , f(x))
2: ((x+b) , f(x+b))

傾きは、

(f(x+b)-f(x)) / b

ですから、これを、平方根の式である f(x) = x^2 - a に当てはめてみると、

( f(x+b)-f(x) ) / b
= ( (x+b)^2 - a - (x^2 - a) ) / b
= ( x^2 + 2bx + b^2 - a - x^2 + a ) / b
= ( 2bx + b^2 ) / b
= ( 2x + b )


となります。ここで、b をどんどん小さくし、0に近づけていったときの値が、 x 地点の接線の傾きとなります。
つまり、

2x が傾きです。

よって、

xm = xn - (f(xn) / g(xn))
= xn - (xn^2 - a) / 2xn

※ xn でひとつの変数です。(x * nではありませんので注意してください)

となり、これを使って次々と新しい xm を求めていけばよいことになります。変化する値が十分に小さくなれば、 この xm が√a の 近似値となります。

これをコード化にしたものを以下に示します。

static double SquareRoot(double a, double xn = 1.0) {
    double e = (xn * xn - a) / (2 * xn);
    if (Abs(e) > 1.0e-14)
        return SquareRoot(a, xn - e);
    return xn - e;
}

static double Abs(double n) {
    return n < 0 ? -n : n;
}


√20 を求める場合は、以下のようにします。

Console.WriteLine(SquareRoot(20.0));

同様の考えで、立方根についても

g(x)
= ( f(x+b)-f(x) ) / b
= ( (x+b)^3 - a - (x^3 - a) ) / b
= ( x^3 + 3x^2b + 3xb^2 + b^3 - x^3 ) / b
= ( 3x^2b + 3xb^2 + b^3 ) / b
= ( 3x^2 + 3xb + b^2 )
となり、b を 0 にして、g(x) = 3x^2 となるので、

xm = xn - (xn^3 - a) / 3xn^2

が成り立ちます。
※ xn でひとつの変数です。(x * nではありませんので注意してください)

ここまで理解すれば。C#のコードにするのは難しくないですね。

static double CubeRoot(double a, double xn = 1.0) {
    double e = (xn * xn * xn - a) / (3 * xn * xn);
    if (Abs(e) > 1.0e-13)
        return CubeRoot(a, xn - e);
    return xn - e;
}

static double Abs(double n) {
    return n < 0 ? -n : n;
}

再帰使ってますが、ループ処理にして以下のように書いても答えは求まります。

static double CubeRoot(double a) {
    double x = 1.0;
    double e;
    do {
        e = (x * x * x - a) / (3 * x * x);
        x = x - e;
    } while (Abs(e) > 1.0e-13);
    return x;
}

  

Posted by gushwell at 22:00Comments(0)TrackBack(0)

2015年02月08日

C#で与えられた文字列でピラミッド

   このエントリーをはてなブックマークに追加 Clip to Evernote
どう書く?orgに感謝を込めて」シリーズ その68

「ピラミッドを作る」の続編です。 与えられた文字列を使って下の例のようなピラミッドを書いてください。 頂点は与えられた文字列の最後の一文字、 底辺は与えられた文字列の各文字の間に空白が入ったものになります。
>>> pyramid("hoge")
   e  
  g e 
 o g e
h o g e



>>> pyramid("abracadabra")
          a         
         r a        
        b r a       
       a b r a      
      d a b r a     
     a d a b r a    
    c a d a b r a   
   a c a d a b r a  
  r a c a d a b r a 
 b r a c a d a b r a
a b r a c a d a b r a


※ ごめんなさい、出題者がわかりません。

いろいろな書き方が考えられますね。 最初に書いたバージョンはこれ。
static void Pyramid1(string s) {
    int length = s.Length;
    for (int n = 0; n < length; n++) {
        // 左側の空白を表示
        int spCount = length - n - 1;
        Console.Write(new string(' ', spCount));
        // 右側の文字部分を表示
        var ls = s.Skip(spCount);
        foreach (var c in ls)
            Console.Write(c + " ");
        Console.WriteLine();
    }
}
メソッドの中で出力も行っているのが気に入らなかったので、 文字列の組み立てだけを行うメソッドにして見ました。 ついでに、上記メソッドの最後の foreach文の箇所をLINQ to Ojectを使って 書き直してみました。
static IEnumerable<string> Pyramid2(string s) {
    int length = s.Length;
    for (int n = 0; n < length; n++) {
        // 左側の空白の数
        int spCount = length - n - 1;
        // 右側の文字部分を表示
        yield return s.Skip(spCount)
                      .Aggregate(new string(' ', spCount), (x, c) => x + c + " ");
    }
}
出力までやるには、以下のようなコードを書きます。
        static void Main(string[] args) {
            Pyramid2("abracadabra").ToList().ForEach(Console.WriteLine);
            Console.ReadLine();
        }
ここまで書いたのだから、全部 LINQ使っちゃえ、って書いたのが次に示すコード。 ここまでやると、黒魔術的なにおいがしてきます。
static IEnumerable<string> Pyramid3(string s) {
    return Enumerable.Range(0, s.Length).Reverse()
            .Select(
                n => s.Skip(n).Aggregate(new string(' ', n), (x, c) => x + c + " ")               
            );
}
再帰処理だと、どんな感じになるのかなと想い、書いてみたのが次のコード。 これが一番すっきりしている感じがします。
static string Pyramid(string s) {
    return PyramidRec(s.Length - 1, s);
}

// indexは、abracadabra の場合は、10,9,8,7,... と変化していく。
// 左側に詰める空白の数
static string PyramidRec(int index, string s) {
    if (index < 0)
        return "";
    string line = new string(' ', index) +
                  s.Skip(index).Aggregate("", (x, c) => x + c + " ");
    return line + Environment.NewLine + PyramidRec(index - 1, s);
}
再帰処理でのyield return って煩雑なので、戻り値を string にして、 改行込みの結果を返すようにしています。 もちろん、コンソールに出力するには、以下のように書けばOK.
Console.WriteLine(Pyramid("abracadabra"));
  
Posted by gushwell at 21:30Comments(0)TrackBack(0)

2014年12月21日

C#で自然数をm個の和で表す全パターンを求める

   このエントリーをはてなブックマークに追加 Clip to Evernote
どう書く?orgに感謝を込めて」シリーズ その58

この問題のオリジナルタイトルは、「自然数の分割」です。

■問題 (出題者:herumiさん)
自然数nとm(n>=m>0)が与えられたとき,nをm個の非負の整数の和で表すやり方を全て出力してください. その際,和の組(x_1, ..., x_m)は大きい順に出力してください. ここでm = 3の時の「(a, b, c)が(A, B, C)より大きい」とは
(a > A)
(a == A) かつ (b > B)
(a == A) かつ (b == B) かつ (c > C)のいずれかが成り立つとき(つまりは辞書的順序)とします.
    
例:n = 5, m = 3が与えられたときは
5, 0, 0,
4, 1, 0,
4, 0, 1,
3, 2, 0,
3, 1, 1,
...
0, 1, 4,
0, 0, 5,
を出力する.
以下、C#のコード。
コメントにある通り、再帰処理で解を求めています。
Solveの第3引数は、答えが格納されるリストです。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Doukaku.Org {
    class Program {
        static void Main(string[] args) {
            Solver sol = new Solver();
            sol.Solve(5,2);
            Console.ReadLine();
        }
    }

    // 再帰呼び出しを使ってる。
    class Solver {
        public void  Solve(int n, int m) {
            Solve(n,m, new List<int>());
        }

        public void Solve(int n, int m, List<int> ans) {
            if (m == 0) {
                if (n == 0) {
                    // めんどくさいので、答えが見つかった時点で印刷。
                    // コンソールアプリ専用。
                    ans.ForEach(x => Console.Write("{0}, ", x));
                    Console.WriteLine();
                }
            } else {
                // 大きい順に出力するために、n, n-1, n-2... と処理していく
                for (int k = n; k >= 0; k--) {
                    ans.Add(k);
                    Solve(n-k, m-1, ans);
                    ans.RemoveAt(ans.Count - 1);
                }
            }
        }
    }
}

  
Posted by gushwell at 22:00Comments(0)TrackBack(0)

2014年08月17日

C#で隣り合う二項の差を求める

   このエントリーをはてなブックマークに追加 Clip to Evernote
どう書く?orgに感謝を込めて」シリーズ その33

■問題 (出題者:にしお さん)
整数のリストがxsが与えられたときに、隣り合う2要素の差のリストを作る関数diffを作ってください。
サンプル入力
[3, 1, 4, 1, 5, 9, 2, 6, 5]
サンプル出力
[-2, 3, -3, 4, 4, -7, 4, -1]  
どの問題もそうですが、いろんな書き方が考えられます。
この問題では4つのバージョンを書いてみました。

最初は、IList<int> を受け取って、添え字を使って各要素にアクセスするもの。 配列など他のコレクションを受け取れないのが欠点です。

■C#のコード(その1)
using System;
using System.Collections.Generic;
using System.Linq;

namespace Doukaku.Org {
    class Program {
        static void Main(string[] args) {
            var nums = new List<int> { 3,1,4,1,5,9,2,6,5 };
            var diff = Diff(nums);
            PrintList(diff);
            Console.ReadLine();
        }

        // IList<int>を受け取り、IList<int>を返す例
        static IList<int> Diff(IList<int> xs) {
            List<int> result = new List<int>();
            for (int i = 0; i < xs.Count() - 1; i++) {
                result.Add(xs[i + 1] - xs[i]);
            }
            return result;
        }

        static void PrintList(IEnumerable<int> nums) {
            string s = string.Join(", ", nums.Select(x => x.ToString()).ToArray());
            Console.WriteLine("[{0}]", s);
        }
    }
}


2つめは IEnumerable<int>を受け取り、 yield return で結果を列挙する例 (MoveNext、Currentを使用) やっぱり、MoveNext, Currentは面倒くさいですね。

■C#のコード(その2)
        static IEnumerable<int> Diff(IEnumerable<int> xs) {
            var ite = xs.GetEnumerator();
            if (ite.MoveNext())
                for (int prev = ite.Current; ite.MoveNext(); prev = ite.Current) {
                    yield return ite.Current - prev;
                }
        }

次は、IEnumerable<int>を受け取り、 yield return で結果を列挙する例 (foreachを使用) ループの中で、最初か最初でないかを判断しているところが気に入らないです。

■C#のコード(その3)   
        static IEnumerable<int> Diff(IEnumerable<int> xs) {
            bool first = true;
            int prev = int.MinValue;
            foreach (var n in xs) {
                if (!first)
                    yield return n - prev;
                prev = n;
                first = false;
            }
        }

そして最後が、再帰を使ったもの。 メソッドが2つに分かれてしまうのが、玉に瑕ですね。

■C#のコード(その4)
static IEnumerable<int> Diff(IEnumerable<int> xs) {
    if (xs.Any())
        return Diff5(xs.First(), xs.Skip(1));
    return new int[] { };
}

static IEnumerable<int> Diff(int prev,IEnumerable<int> xs) {
    if (xs.Any()) {
        int b = xs.First();
        yield return b - prev;
        foreach (var x in Diff5(b, xs.Skip(1)))
            yield return x;
    }
}

どれも一長一短ですね。あなたならどう書く?
  
Posted by gushwell at 23:00Comments(0)TrackBack(0)

2011年11月13日

パズル Magic Star をC#で解く

   このエントリーをはてなブックマークに追加 Clip to Evernote
WebSite に、
パズル Magic Star を解く
の記事を新たに追加しました。

MagicStar とは、下の図の12 個ある○に 1 から 12 までの数字をひとつずつ入れていき、
直線上の4個の数字の合計が、すべて 26 になるように、数字を配置するというパズルです。

MagicStar01


これを C# で解いています。
Silverlight で作成しているので、実際に動かして答えを見ることができます。
  
Posted by gushwell at 22:26Comments(0)TrackBack(0)

2011年09月20日

F#:どう書く?org 「文字列の均等分割」

   このエントリーをはてなブックマークに追加 Clip to Evernote
どう書く?orgの「文字列の均等分割」をF#で書いてみました。
一行の文字列を指定した数の行にできるだけ文字数が均等になるように分割してください.
ただし,除算や剰余算を使わないで書いてみてください.

sample = "ゆめよりもはかなき世のなかをなげきわびつゝあかしくらすほどに四月十よひにもなりぬれば木のしたくらがりもてゆく"

divid 4 sample =>
"ゆめよりもはかなき世のなかを"
"なげきわびつゝあかしくらすほ"
"どに四月十よひにもなりぬれ"
"ば木のしたくらがりもてゆく"


問題では、「除算や剰余算を使わないで」とありましたが、それは無視して書いています。
なので、1行の幅の文字数は、簡単に求められるので、
あとは、String.Substringメソッドで、文字列を分割して行くだけです。
気をつける点といえば、「できるだけ均等に」ということなので、1行の文字数は、
最初に求めるのではなく、再帰関数のなかでその都度求めている点でしょうか。


問題での例では、divide関数内で、プリントもしていましたが、
ここでは、分割する関数ではプリントはせずに、分割した結果を返すようにしています。
このdivide関数は、string listを返す関数で、再帰呼び出しにより実現しています。
:: 演算子で、 左側の要素と右側のリストを連結しています。


C#のネーミング規則に慣れてしまっているので、ついつい、関数名を大文字で始めてしまいます。
規模の大きなプログラムでは、すべてが小文字で始まると、どれが関数でどれが変数なのかの
判断が付きにくいと思うのですが、他の人はそんなふうには感じないのかな...
まあ、慣れの問題なのでしょうが、そもそも2,3週間に一回程度、F#をいじるだけですから、
慣れようがないです。

  
Posted by gushwell at 22:00Comments(0)TrackBack(0)

2011年09月13日

コッホ曲線をSilverlightで作成

   このエントリーをはてなブックマークに追加 Clip to Evernote
コッホ曲線を描くプログラムをSilverlight で作成しました。
コッホ曲線とは、フラクタル図形の一種で、線分を3等分し、分割した真ん中の線分を底辺とする
正三角形を描く(ただし底辺は消す)ことを無限に繰り返すことによって得られる図形です。


24


実際に動くプログラムと、そのソースコードは、Gushwell's C# Programming Page - コッホ曲線 をご覧ください。


ある直線(線分)から、次の点(正三角形の頂点)を求めるロジックが意外と面倒でした。
角度の変換とか、2点間の長さを求めるとか、Sin, Cos関数を使うとかは 普段のプログラムではやっていないので...

  
Posted by gushwell at 22:49Comments(0)TrackBack(0)

2011年07月03日

C#プログラマーのための再帰処理・超入門

   このエントリーをはてなブックマークに追加 Clip to Evernote
いげ太さんの記事に触発されて僕も書いてみた。

岩波国語辞典には、再帰→回帰 とあり、回帰の意味を見てみると、

(2) 〔名〕処理手続きや規則の定義に、それ自身を繰り返し使うような仕方。

とあります。
つまり、C#での再帰処理とは、ある処理をするのに、自分自身のメソッドを繰り返し呼び出す処理のことです。

なお、再帰は、処理だけではなく、構造においても使われます。コンピュータのフォルダ構造が
その身近な例ですね。

さて、それでは、1 から n までの整数の合計(これを sum(n)と表す)を求める処理
について考えてみます。


と表せます。

1 から (n-1) までの合計は、というと、sum(n-1) であらわすことができます。

つまり、


と定義できるわけです。
これが、分かれば、再帰コードを書くことが出来ます。

次のようなコードを書いてみました。


何をやっているかと言えば、
Sum(n) を求めるには、Sum(n-1) で自分自身のメソッドを呼び出し、1から n-1 の
合計を求め、その結果に n を加えた値を結果として、返しています。

しかし、これには、大きな欠陥があります。
実際に、動かしてもらえば、分かると思いますが、スタックオーバーフローの例外が発生してしまいます。

なぜならば、Sumメソッドが呼ばれたら、その中で、Sumメソッドを呼び出し、さらに、
その中で、Sumメソッドを呼び出し.... を繰り返すわけですから、無限に、Sum を呼び
続けてしまうわけです。

さて、ここで、While 文を使った、同様のコードを見てみましょう。


このコードでは、無限ループにならないように、i < n でループの終了判定を行っていますね。

再帰のコードも同様に、再帰の終了判定を行わないと、無限にメソッドを呼び続ける
ことになるのです。

ためしに、


とし、実行してみましょう。Sum メソッドに渡ってきた 引数 n をConsole.WriteLineで
表示するコードと、Sleepメソッドを追加しています。

※ 途中で、Ctrl+Cを押して、プログラムを終了させてください。


と表示され、Sumの引数に 0 以下の値が渡ってきていることが確認できます。
これを見れば、1 がきたところで、再帰呼び出しをストップすれば良いことがわかります。

それでは、再帰版のSumメソッドに終了判定を入れたいと思います。


※ 一時変数の resultは取り除いています。

問題は、コメントにも書いたように、n が 1の時に何を書くかです。

具体的に、Sum(1) の結果が何かを考えてみると、1 が返ればよいわけですから、
最終的なコードは、次のようになります。

なお、最初に、Sum(n)は、


と定義できると書きましたが、これはより正確には、


となります。
この定義をそのまま、C#のコードに書き下したものが、上記Sumメソッドとなります。
これは、メソッドの定義が、手続き的な「処理の記述」から、問題の構造そのものを表す
「問題の定義」に近くなっていることを示しています。(言い換えると、HowからWhatですね)
コード自体も、とてもすっきりしていますね。whileをつかったコードのごちゃごちゃ感が
こちらにはありません。

ただ、手続きを書くことに慣れているプログラマーには、これでは、どう動くのかが
良くわからないという方も多いと思います。実際に、僕が若いころに再帰を学んだ時はそうでした。

ということで、どう動いているのかを確認してみましょう。
以下のようにデバッグ行を追加してみます。


この結果は、
となります。
これを言葉で言い換えると、次のようになるかと思います。


1. Sum(5)を求めるには、そのひとつ前のSum(4)を知らなければならない、だからSum(4)を呼ぶ。
2. Sum(4)を求めるには、そのひとつ前のSum(3)を知らなければならない、だからSum(3)を呼ぶ。
3. Sum(3)を求めるには、そのひとつ前のSum(2)を知らなければならない、だからSum(2)を呼ぶ。
4. Sum(2)を求めるには、そのひとつ前のSum(1)を知らなければならない、だからSum(1)を呼ぶ。
5. Sum(1) は、1だから、1 を返す。
6. Sum(1)が求まったから、Sum(2) は、Sum(1) + 2 で 1 + 2 になるので、3を返す。
7. Sum(2)が求まったから、Sum(3) は、Sum(2) + 3 で 3 + 3 になるので、6を返す。
8. Sum(3)が求まったから、Sum(4) は、Sum(3) + 4 で 6 + 4 になるので、10を返す。
9. Sum(4)が求まったから、Sum(5) は、Sum(4) + 5 で 10 + 5 になるので、15を返す。


これからわかるように、C#における再帰呼び出しは、メソッドをどんどん深く呼び出し続けるため、
スタックを消費することになります。
大抵の場合は、これでも問題はありませんが、解くべき問題によっては、スタックオーバーが
起こりえます。この点は注意が必要です。

しかし、その欠点を補ってあまりあるパワーが再帰にはあります。
それは、再帰的な構造のデータを処理する際に、発揮されます。
先ほどあげた、フォルダ構造を扱ったり、家系図のような階層構造(Tree構造)のデータを扱う際には、
再帰はとても有効な手段となります。

ボードゲームのような先読みを処理を行うにも、再帰は有効(というか必須)です。
また、リンクリストを扱ったり、今回のような通常の繰り返し処理や、リトライ処理などでも
再帰が使えますね。

なお、いげ太さんの記事では、末尾再帰によるコードが示されていますが、
C#のコンパイラには、末尾再帰の最適化が組み込まれていません。 そのために、ここでは、末尾再帰のコードは示していません。
末尾再帰について知りたければ、こことかここを見てください。



ちなみに、僕のWebサイト「Gushwell"s C# Programming Page」の「C# プログラム小品集」
に掲載したプログラムでも、再帰を使っているプログラムがたくさんあります。
例えば、

階乗の計算
http://gushwell.ifdef.jp/etude/Factorial.html

すべての順列を求める
http://gushwell.ifdef.jp/etude/Permutation.html

ヒルベルト曲線
http://gushwell.ifdef.jp/etude/HirbertCurve.html

8-Queenパズル
http://gushwell.ifdef.jp/etude/nQueen.html

最大面積の領域を求める
http://gushwell.ifdef.jp/etude/MaxArea.html

などです。興味がありましたら、読んでみてください。

  
Posted by gushwell at 23:00Comments(0)TrackBack(1)