例によってありきたりな話ですが......自分用のメモも兼ねて。
#NaNとは？
浮動小数点数にはNaN (Not a Number, 非数) と呼ばれる、実数の異常な値を表す特殊な数があります。
これは無限大-無限大、0.0/0.0といった不定形や、負数の平方根、負数の対数といった実数で表せない計算を行った場合に発生します。
以下いくつかの例。
// C++
#include <iostream>
#include <limits>
#include <cmath>
int main() {
double inf = std::numeric_limits<double>::infinity();
double a = inf-inf, b = 0.0/0.0, c = sqrt(-1.0), d = log(-1.0);
std::cout << a << ", " << b << ", " << c << ", " << d << "\n";
return 0;
}
// Java
class Nan1 {
public static void main(String[] args){
double inf = Double.POSITIVE_INFINITY;
double a = inf-inf, b = 0.0/0.0, c = Math.sqrt(-1.0), d = Math.log(-1.0);
System.out.println(a + ", " + b + ", " + c + ", " + d);
}
}
// Swift
import Foundation
let inf = Double.infinity
let a = inf-inf, b = 0.0/0.0, c = sqrt(-1.0), d = log(-1.0)
print(a, ", ", b, ", ", c, ", ", d)
# Python
import numpy as np
inf = np.inf
a = inf-inf
# 0.0/0.0 は ZeroDivisionError
b = np.sqrt(-1.0)
c = np.log(-1.0)
print(a, ", ", b, ", ", c)
言語によって"nan"が"NaN"だったり微妙な違いはありますが、結果は概ね以下のようになります。
nan, nan, nan, nan
なお、Pythonでは0.0/0.0はZeroDivisionErrorが発生します。
#ちなみに、整数の場合
浮動小数点数と違って、整数には通常の方法で無限大やその他異常な値を格納する方法がありません。これは、正確に言うと整数を表すbit表現の中にそういった数を表すbitが割り当てられていない、ということです。
そのため整数のゼロ除算はCやC++では未定義、その他大抵の言語では例外が発生します。
例えば、以下のC++のコードを考えてみます。
// C++
#include <iostream>
int main() {
int a = 1/0;
std::cout << "1/0 = " << a << "\n";
return 0;
}
見るからにアウトなコードですが、これをコンパイルすると、例えば以下のような警告を受けると思います。
GCCの場合。
zero_division1.cpp: In function 'int main()':
zero_division1.cpp:5:12: warning: division by zero [-Wdiv-by-zero]
int a = 1/0;
~^~
Clangの場合。
zero_division1.cpp:5:12: warning: division by zero is undefined [-Wdivision-by-zero]
int a = 1/0;
^~
1 warning generated.
せっかくの警告も無下に実行すると、GCCの場合、
$ ./a.out
Floating point exception: 8
と実行時例外が、
Clangの場合、
$ ./a.out
1/0 = 215109936
$ ./a.out
1/0 = 59375920
$ ./a.out
1/0 = 210141488
$ ./a.out
1/0 = 234668336
$ ./a.out
1/0 = 8167728
と実行の度にデタラメな値が出力されるようになってしまいました。
これは環境によっても違うと思いますが、ともかく未定義動作なので何が起きても文句は言えません。
その他のもう少し安全な(?)言語なら、例えば以下のJavaのように例外が発生します。
// Java
class ZeroDivision {
public static void main(String[] args) {
int a = 1/0;
System.out.println("1/0 = " + a);
}
}
$ java ZeroDivision
Exception in thread "main" java.lang.ArithmeticException: / by zero
at ZeroDivision.main(ZeroDivision.java:4)
#NaNの性質
話を浮動小数点数のNaNに戻しましょう。
NaNには以下のような性質があります。
- NaNとの四則演算は常にNaN
a + b
a - b
a * b
a / b
はaまたはbがNaNだったとき、常にNaNを返します。
Javaなら以下の通りです。
// Java
class Nan2 {
public static void main(String[] args) {
double a = 1.0, b = Double.NaN;
System.out.println((a+b) + ", " + (a-b) + ", " + (a*b) + ", " + (a/b));
}
}
$ java Nan2
NaN, NaN, NaN, NaN
ので、数値計算をやっていると、ある点で発生したNaNが瞬く間に伝染していっていつの間にか出力データがNaNで埋め尽くされていた......なんてことがよくあります。
- NaNとの!=以外の比較演算は常にfalse、!=は常にtrue
a == b
a > b
a < b
a >= b
a <= b
はaまたはbがNaNだったとき、常にfalseを返し、また
a != b
はaまたはbがNaNだったとき、常にtrueを返します。
Javaなら以下の通りです。
// Java
class Nan3 {
public static void main(String[] args) {
double a = 1.0, b = Double.NaN;
System.out.println((a == b) + ", " + (a > b) + ", " + (a < b) + ", "
+ (a >= b) + ", " + (a <= b) + ", " + (a != b));
}
}
$ java Nan3
false, false, false, false, false, true
特に、NaN == NaNはfalseです！
なので、xがNaNであるかどうか調べるときのありがちな間違いとして、
if (x == Double.NaN) {
// Do something...
} else {
// Do something...
}
というコードではNaNを検出できません(常にelseの方が実行されてしまいます)。
正しくは、大抵の言語にはisnan(x)とかDouble.isNaN(x)とかいった専用のメソッドが用意されていますので、
if (Double.isNaN(x)) {
// Do something...
} else {
// Do something...
}
というように書きましょう。
また、上の性質を逆手にとって、
if (x != x) {
// Do something...
} else {
// Do something...
}
としてもNaNを検出できます。
実際、ECMAScript6以前のJavaScriptではこれと似た方法が使われていました(NaNの判定について、NaNはNot a NumberだけどNumber型である話)。
また、aやbがNaNになりうる場合、a > b と !(a <= b) が等価でないことにも注意が必要です。
#数値計算での例
以上のことを踏まえて、大学の数値計算の授業で実際に経験した話を一つ紹介します。
連立方程式を解くアルゴリズムの一つに、ガウス=ザイデル法(1 連立方程式の解法 (反復法),
連立1次方程式：反復法 - PukiWiki for PBCG Lab)というものがあります。
詳細は上のリンクに譲るとして、C++による以下の二つの実装を見てみましょう。
- 実装1
#include <array>
#include <cmath>
template<size_t N>
using vec = std::array<double, N>;
template<size_t N>
using mat = std::array<std::array<double, N>, N>;
template<size_t N>
vec<N> gauss_seidel1(const mat<N>& A, const vec<N>& b) {
constexpr double epsilon = 1.0e-8;
vec<N> x1, x2;
for (size_t i = 0; i < N; i++) {
x1[i] = x2[i] = 0.0;
}
double d;
do {
d = 0.0;
for (size_t i = 0; i < N; i++) {
double sum = b[i];
for (size_t j = 0; j < N; j++) {
if (i != j) {
sum -= A[i][j]*x2[j];
}
}
x2[i] = sum/A[i][i];
d += abs(x1[i]-x2[i]);
}
x1 = x2;
} while (d > epsilon);
return x2;
}
- 実装2
#include <array>
#include <cmath>
template<size_t N>
using vec = std::array<double, N>;
template<size_t N>
using mat = std::array<std::array<double, N>, N>;
template<size_t N>
vec<N> gauss_seidel2(const mat<N>& A, const vec<N>& b) {
constexpr double epsilon = 1.0e-8;
vec<N> x1, x2;
for (size_t i = 0; i < N; i++) {
x1[i] = x2[i] = 0.0;
}
while (true) {
double d = 0.0;
for (size_t i = 0; i < N; i++) {
double sum = b[i];
for (size_t j = 0; j < N; j++) {
if (i != j) {
sum -= A[i][j]*x2[j];
}
}
x2[i] = sum/A[i][i];
d += abs(x1[i]-x2[i]);
}
if (d <= epsilon) break;
x1 = x2;
}
return x2;
}
あまり細かいところは気にしなくて構いません。見ていただきたいのはdの扱いです。このdは、一ステップ前のベクトルをx1、現在のステップのベクトルをx2としたときの「距離(正確にはL1ノルム)」を表しています。
そして、前者の実装はこのdがepsilonより大きい間は反復を継続し、後者は逆にdがepsilon以下になったら反復をやめるという形になっていて、結局のところ、やっていることはどちらもほとんど同じに見えます。
実際、どちらの実装も
\left(
\begin{array}{ccc}
3 & 1 & 1 \\
1 & 3 & 1 \\
1 & 1 & 3
\end{array}
\right)
\bf{x}
=
\left(
\begin{array}{c}
0 \\
4 \\
6
\end{array}
\right)
の解
\bf{x} =
\left(
\begin{array}{c}
-1 \\
1 \\
2
\end{array}
\right)
を正しく与えます。
しかしこのガウス=ザイデル法には欠点があります。それは、
\left(
\begin{array}{ccc}
1 & 2 & 2 \\
2 & 1 & 2 \\
2 & 2 & 1
\end{array}
\right)
\bf{x}
=
\left(
\begin{array}{c}
1 \\
0 \\
-1
\end{array}
\right) \\
\bf{x} =
\left(
\begin{array}{c}
-1 \\
0 \\
1
\end{array}
\right)
のような、解を一意にもつ連立方程式であっても、係数行列が対角優位でないと解が求められないことがあるということです。
具体的に言うと、上のコードのx1[i]やx2[i]の値が発散して、しまいにはNaNになってしまいます。
すると、x1、x2を求めるのにx1、x2を使っているので、上の性質1.より一度NaNが発生したらそれ以降はずっとNaNのままです。なので、x1とx2から求めたdも以降ずっとNaNになります。したがって、性質2.より d > epsilon と d <= epsilon はともにfalseのまま変わらなくなります。
......
もうお気づきでしょうか。
解をうまく求められず発散してしまった場合、実装1ではループから抜け出せますが、実装2では無限ループになってしまいます。
このように、アルゴリズムが正しくても入力する値によってはうまくいかないことがあるので、そうなる危険がある場合は適切に対策を取っておく必要があります(上の場合なら、最大反復回数を定めておく、うまく求められない場合は例外を投げるなど)。