ナムアナ

• 10 進整数を 2 で割って商と余剰を求める
• 商が 0 になるまで商に対して除算を繰り返す
• 手順1から手順2で求めた余剰を逆に並べる

実装

11進数以上の基数に変換する場合にはアルファベットが必要となる． そこで，配列に0からZまでの値を順に格納しておくことで，11以上の値をアルファベットに変換可能にする．

10 進整数を任意の基数(2進数から36進数まで)に変換するプログラムを以下に示す．

/* header files */
#include <stdio.h>
#include <stdlib.h>

/**
 * 基数変換を行う
 * @param[IN] x 変換する値
 * @param[IN] radix 変換する基数
 * @param[OUT] result 変換した値を逆順に格納
 * @return 変換後の桁数
 */
int ConvertRadix(unsigned int x, int radix, char *result) {
    char digits[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
    int n = 0;

    /* 変換する値が0ならそのまま0を格納する */
    if ( x == 0 ) {
        result[n] = 0;
        n++;

    } else {
        while ( x > 0 ) {
            result[n] = digits[x % radix]; /* 余剰を格納 */
            x /= radix;
            n++;
        }
    }
    return n;

}

/* main */
int main(void) {
    unsigned int x = 256; /* 変換する値 */
    int radix = 2; /* 基数 */
    char result[256] = {'\0'}; /* 変換後の値 */
    int n; /* 変換後の値の桁数 */
    int i;

    /* 基数変換 */
    n = ConvertRadix(x, radix, result);
    for ( i = n - 1; i >= 0; i-- ) {
        printf("%c", result[i]);
    }
    printf("\n");

    return EXIT_SUCCESS;
}

ConvertRadixは余剰を求めた順に配列に格納しているため，呼び出し側 (main) で逆順に表示していることに注意．

LU分解とは

行列Aを下三角行列(lower triangular matrix)Lと上三角行列(upper triangular matrix)Uの積に分解する。

LU分解してみる

例えば以下のような行列Aがあったとする。

これをLU分解すると以下のようになる。

LU分解を使って連立方程式を解く

簡単のため、以下のようなAx = bを解くことにする。

とりあえずAをLU分解する。

ここで以下のようにUx = yとおくと・・・

このように表せる

これを解くと、yが求まる。(y1 = 5, y2 = 3, y3 = 1)

Ux = yとおいたので、これに上で求めたyの値を代入すると以下のようになる。

これを解くと解が求まる。x1 = 3, x2 = 1, x3 = 1

ガウスの消去法とは

連立方程式の解を機械的に求めたい。こんなときに使うのがガウスの消去法

考え方

例えば以下のような連立方程式があったとする。

ガウスの消去法を使ってこの連立方程式を解いてみる。

ガウスの消去法前半部

(1)式と(2)式のx1の係数に注目する。

以下の計算を行う。

次に(1)式と(3)式のx1の係数に注目する。

以下の計算を行う。

このように計算を繰り返すと、以下のようになる。

(第1段階終了)

今度は(2)式と(3)式のx2の係数に注目する。

以下の計算を行う。

・・・

これらの計算を繰り返すと、以下のような式に変形できる。

ガウスの消去法後半部

(4)式からx4の値を求めることが出来る。(x4 = 2)
(3)式のx4に値を代入するとx3が求まる。(x3 = 1)
(2)式のx3に値を代入するとx2が求まる。(x2 = 1)
(2)式のx2に値を代入するとx1が求まる。(x1 = 2)

以上のように解を求めることができる。

実装

連立方程式は行列で表せる。行列は2次元配列で実装できる。 以下のように、連立方程式を2次元配列にセットし、前半部と後半部を実装すればOK

n = 4;

/* 前半部 */
for( k = 1; k <= n - 1; ++k ){
    for( i = k + 1; i <= n; ++i ){
        m = a[i][k] / a[k][k];
        /* 方程式から方程式を引く */
        for( j = k + 1; j  <= n + 1; ++j ){
            a[i][j] = a[i][j] - m * a[k][j];
        }
    }
}

/* 後半部 */
for( i = n; i >= 1; --i ){
    s = 0.0;
    for( j = i + 1; j <= n; ++j ){
        s = s + a[i][j] * x[j];
    }
    x[i] = (a[i][n + 1] - s) / a[i][i];
}

ニュートン法(Newton's method)とは

f(x) = 0の解を求める方法の1つ

考え方

例えば以下のような多項式方程式があったとする。

上記の式の解を求めたい。 これをニュートン法でやると、以下のような感じになる。

• 適当な1点x0を決める
• 関数f(x)上の点(x0,f(x0))での接線を引く
• x軸と上記で引いた接線の交点を新しくx1とする
• 関数f(x)上の点(x1,f(x1))での接線を引く
• x軸と上記で引いた接線の交点を新しくx2とする
• 関数f(x)上の点(x2,f(x2))での接線を引く
• x軸と上記で引いた接線の交点を・・・(繰り返し)

イメージ

ニュートン法の漸化式

x0からx1を求めるには、以下のような考え方をすればOK

なので、変形するとx1が求まる。

一般化すると以下のように書ける。

実装

ニュートン漸化式をプログラムに落とし込むだけ。

do {
    e = (exp(x) - x * x - 2.0) / (exp(x) - 2.0 * x);
    x = x - e;
} while ( fabs(e) > 1.0e - 8 );
printf("%16.8e\n", x);

2分法(bisection method)とは

f(x) = 0の解を求める方法の1つ

考え方

例えば以下のような多項式方程式があったとする。

上記の式の解を求めたい。こんな場合に使えるのが2分法。(ただし、速度は遅い)
2分法では次のような手順で求める。

• 適当な2点a,bを決める
• 線分a,bの中点cをとる
• f(c) > 0ならc → bに、f(c) < 0なら c → aに置き換える
• 手順2,3を繰り返す

イメージ

実装

手順1〜4をC言語に落とし込めばOK

/* f(x) */
double f( double x ) {
    return exp(x) - x * x - 2.0;
}

・・・

    do {
        c = (a + b) * 0.5;
        if ( f(c) * f(a) < 0.0 ) {
            b = c;
        } else {
            a = c;
        }
    } while( fabs(b-a) > 1.0e-12 );

ホーナー法とは

n次多項式の演算回数を減らす手法。

考え方

例えば以下のような多項式があったとする。

この多項式をC言語で実装すると、以下のように書ける。

P3 = a[3] * x * x * x + a[2] * x * x + a[1] * x + a[0];

でもこれだと、乗算が多すぎて実行速度が遅い。(乗算の回数: 6回)
そこで、もうちょっと計算量を減らせないものかとHornerさんは考えた。 そして以下のような式に変形したらいいじゃん。という結論に達したわけである。

P3 = ((a[3] * x + a[2]) * x + a[1]) * x + a[0];

上記の式に変形することで乗算の回数が3回に減少する。
なお、C言語のmath.hには、べき乗を計算するpow関数が用意されているが、一般にホーナー法を使ったほうが速いとされている。

実装

次数が決まっている場合は、上記のようにホーナー法によって変形した式をそのままプログラムに書いてしてしまえばよい。
次数が以下のように決まっていない場合(実数が変数になっている場合)はループで回してやる必要がある。

次数が未定の場合の実装例

/* 次数の入力 */
scanf("%d", &n);
p = a[n];
for (i = n - 1; i >= 0; --i) {
    p = p * x + a[i];
}