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月18日

C#でコラッツ・角谷の問題のステップ数(2^20まで)

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

■問題 (出題者:ところてん さん)
任意の数nを与えたときに
・nが偶数ならば2で割る (n=n/2)
・nが奇数ならば3倍して1を足す (n = 3*n+1)
を繰り返すと、いづれは1になる。というものがあります。

数値計算の上ではかなりの数まで成り立つことが知られています。
(すべての数について成り立つかは不明) 参考リンク先参照

ある任意の数nがコラッツ・角谷の問題で1になるまでのステップ数をf(n)とします。
1〜2^20までの数でf(n)を求めて、f(n)が最大になるときのnとf(n)を表示してください。

たとえばn=9だと次のような数列をたどって、19ステップで1になります。
9->28->14->7->22->11->34->17->52->26->13->40->20->10->5->16->8->4->2->1
つまりf(9)=19です。

また、最大を求めた際の実行時間と環境を書いてください。

参考: コラッツの問題の成り立つ範囲
http://q.hatena.ne.jp/1115469752

C#のプログラミング的にもアルゴリズム的にも、普通のコードなので説明の必要はありませんね。
しいて言えば、コラッツ・角谷の予想が成立することが大前提で書かれたコードです。
成り立たないと無限ループに陥ります。まあ、その心配はなさそうですが...
また、オーバーフローは考慮していません。
もうすこし大きな値で確認するには、BigIntegerを使ったほうがよさげですね。


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

namespace Doukaku.Org {
    class Program {
        static void Main(string[] args) {
            Stopwatch sw = Stopwatch.StartNew();
            CollatzKakutani.PrintMaxStep((long)Math.Pow(2, 20));
            sw.Stop();
            Console.WriteLine(sw.ElapsedMilliseconds + "ミリ秒");
            Console.ReadLine();
        }
    }

    static class CollatzKakutani {
        public static void PrintMaxStep(long limit) {
            long maxn = long.MinValue;
            long maxfn = long.MinValue;           
            for (long i = 1; i <= limit; i++) {
                long r = Calc(i);
                if (r > maxfn) {
                    maxfn = r;
                    maxn = i;
                }
            }
            Console.WriteLine("f({0})={1}", maxn, maxfn);
        }

        public static long Calc(long n) {
            long result = 0;
            while (n != 1) {
                result++;
                if (n % 2 == 0)
                    n = n / 2;
                else
                    n = n * 3 + 1;
            }
            return result;
        }
    }
}

僕のPCでの実行結果 (Releaseモード)
f(837799)=524
401ミリ秒
  
Posted by gushwell at 22:30Comments(0)TrackBack(0)

2015年02月15日

C#で正整数のゲーデル数化

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

■問題 (出題者:nobsun さん)
正の整数 n を引数としてとり, 2^d1 * 3^d2 * 5^d3 ... * pk^dk を返す関数 goedel を定義してください.

ただし,n を10進表現で k 桁の数としたときの各桁の数が数列 [d1,d2,d3,...,dk] をなすとし,dk が 1 の位,d1 が 10^(k-1) の位です.また,pk は k番目の素数です.
goedel   9  ⇒ 2^9             ⇒  512
goedel  81  ⇒ 2^8 * 3^1       ⇒  768
goedel 230  ⇒ 2^2 * 3^3 * 5^0 ⇒  108


ゲーデル数とは何かについては、僕はほとんど何も分かってません(^^;;。
でも、定義が明確に書いてあるので、僕にもプログラムが書けますね。
コメントにも書いてありますが、素数を列挙するメソッドは効率無視です。
ここでは、整数 n は、int型としましたので、最大でも、intで表せる桁数分(つまり10個)しか素数を求めないから これで問題ないですね。

戻り値は、BigIntegerにしました。なので、.NET Framework4 以降が対象です。
BigIntegerは、論理上は大きさに上限、下限がないので、int.MaxValueを引数に与えても正しいゲーデル数を返します。

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

namespace Doukaku.Org {
    class Program {
        static void Main(string[] args) {
            Console.WriteLine(Goedel(9));
            Console.WriteLine(Goedel(81));
            Console.WriteLine(Goedel(230));
            Console.WriteLine(Goedel(4321));
            Console.WriteLine(Goedel(54321));
            Console.WriteLine(Goedel(654321));
            Console.WriteLine(Goedel(7654321));
            Console.WriteLine(Goedel(87654321));
            Console.ReadLine();
        }

        static BigInteger Goedel(int num) {
            BigInteger result = 1;
            string s = num.ToString();
            var ite = GetPrimes().GetEnumerator();
            foreach (var c in s) {
                ite.MoveNext();
                result *= BigInteger.Pow(ite.Current, c - '0');
            }
            return result;
        }
        
        // 素数を列挙する。効率は無視。
        static IEnumerable<int> GetPrimes() {
            for (int i = 2; i < int.MaxValue; i++) {
                bool isPrime = true;
                for (int j = 2; j < i; j++)
                    if (i % j == 0) {
                        isPrime = false;
                        break;
                    }
                if (isPrime)
                    yield return i;
            }
        }
    }
}

MoveNext, Current呼び出してるのは煩雑なので、LInq to Object の Zip メソッド使って Goedel メソッドを書き直してみました。
そうなると、当然ながら、Aggregate 使ったほうが楽ですね。
メソッドチェーンがいい感じです。LINQ最高!
ただ、Aggregateは、慣れるまでは、直感的ではないのが玉に瑕ですね。

        static BigInteger Goedel(int num) {
            var gedel = num.ToString()
                           .Select(c => c - '0')
                           .Zip(GetPrimes(), (d, p) => BigInteger.Pow(p, d))
                           .Aggregate((r, n) => r * n);
            return gedel;
        }
       

以下、結果です。

512
768
108
75600
174636000
5244319080000
2677277333530800000
25968760179275365452000000

  
Posted by gushwell at 21:30Comments(0)TrackBack(0)

2015年01月21日

C#で10000以下のダブル完全数をすべて求める

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

■問題 (出題者:にしお さん)
「完全数」とは、自分以外の約数の和が自分自身と等しいような整数のことです。
ここで「自分以外の約数の和が自分自身の2倍と等しいような整数」を「ダブル完全数」と呼ぶことにします。10000以下のダブル完全数をすべて求めるコードを書いてください。

最初に書いたコード。

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

namespace Doukaku.Org {
    class Program {
        static void Main(string[] args) {
           
            for (long i = 0; i <= 10000; i++)
                if (IsDoublePerfect(i)) 
                    Console.WriteLine(i);
                
            Console.WriteLine("end");
            Console.ReadLine();
        }

        static bool IsDoublePerfect(long n) {
            long r = 1;
            for (int i = 2; i <= n / 2; i++) {
                if (n % i == 0)
                    r += i;
            }
            return (n * 2 == r);
        }
    }
}


10000までのダブル完全数は、120, 672 の2つでした。
では、672の次のダブル完全数はいくつかな、と思って、 10000 を 1000000 までにしてみたら遅い。
いつまで待っても終わりません。 とんでもなく遅いです。いくらなんでも遅すぎます(T T)

ということで、ちょっと改良したコード。

static bool IsDoublePerfect(long n) {
    long r = 1;
    long limit = (long)Math.Sqrt(n);
    for (int i = 2; i <= limit; i++) {
        if (n % i == 0)
            r += i + (n / i);
    }
    return (n * 2 == r);
}

例えば、18 の約数を調べるとすると、2で割り切れたら、その商である 9も約数であることは明白なので、
r += i + (n / i);

で、2と9を加えています。当
然、約数かどうかは、(long)Math.Sqrt(n) まで調べれば 良いことになります。
たったこれだけで、随分と速くなりました。
3つ目のダブル完全数が知りたい方は、実際に動かしてみてください。
   
Posted by gushwell at 23:00Comments(0)TrackBack(0)

2010年12月06日

円周率を10000桁まで求める

   このエントリーをはてなブックマークに追加 Clip to Evernote
円周率を10000桁まで求めるプログラムをSilverlightで書いてみました。
BigNumber という多倍長整数演算クラスを作成し、マチンの公式をつかって円周率を求めています。

こちらで、プログラムを動作させることができます。
ソースコードも公開していますので、興味のある方はどうぞ。   
Posted by gushwell at 21:52Comments(0)TrackBack(0)