視覚的に見る浮動小数点数の誤差

数値計算分野の調べ物をしていたら、誤差を視覚化する話題があった。面白いなと思ったので書いておく。

題材

二次方程式 $ax^2+bx+c=0 (a \neq 0)$ の解の公式

$x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$

問題を単純化するために $a=1$ とする。まず真の解 $\alpha, \beta$ を定め、その解に合うような方程式の係数を、解と係数の関係から

$b = -(\alpha + \beta)$

$c = \alpha \beta$

で計算し、そこから解の公式で解を求める。真の解との差が、一連の計算の誤差になる。

実際に視覚化してみる

[-1, 1]の範囲でランダムに500組の真の解を生成し、そこから前述のように解を再計算する。真の解をx,y座標に、真の解と解から求めた相対誤差の和を高さとしてプロットする。

せっかくだから、またMathematicaに頑張ってもらう……と思ってたのだが、いまいち使い方が分からない*1。グラフ出力機としてだけ使おう。

仕方ないからCで書いたコードもなかなか動かない。何が悪いのかと悪戦苦闘するうちに、いつのまにかメルセンヌツイスタが入ったり色々と化けていた*2。突貫工事だが動くのでよしとしよう。

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#include "mt19937ar.c"

#define SWAP(type, a, b) do { type t = a; a = b; b = t; } while(0)

float frand(void) {
    return (float)genrand_real2();
}

int main() {
    unsigned long init[4] = {0x123, 0x234, 0x345, 0x456}, length = 4;
    init_by_array(init, length);

    for(int i = 0; i < 500; i++) {
        float alpha = (2*frand() - 1) + frand() * 10E-6;
        float beta = (2*frand() - 1) + frand() * 10E-6;
        float b, c, s, x, y, x1, y1, e;

        b = -(alpha + beta);
        c = alpha * beta;

        s = sqrt(b*b/4 - c);
        x = -b/2 + s; y = -b/2 - s
        if (alpha < beta) SWAP(double, x, y);

        e = fabs((x - alpha)/alpha) + fabs((y - beta)/beta);

        printf("{%f,%f,%f},\n", alpha, beta, e);
    }
    return 0;
}