As usual, it's a trivial story, but it also serves as a memo for myself.
Floating-point numbers have a special number called NaN (Not a Number) that represents anomalous real numbers. This happens when you perform calculations that cannot be represented by real numbers, such as infinity-infinity, indeterminate forms such as 0.0 / 0.0, square roots of negative numbers, and logarithms of negative numbers.
Below are some examples.
nan1.cpp
// 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;
}
Nan1.java
// 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);
}
}
nan1.swift
// 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)
nan1.py
# Python
import numpy as np
inf = np.inf
a = inf-inf
# 0.0/0.0 is ZeroDivisionError
b = np.sqrt(-1.0)
c = np.log(-1.0)
print(a, ", ", b, ", ", c)
Depending on the language, "nan" may be "NaN" or there are subtle differences, but the results are roughly as follows.
nan, nan, nan, nan
In Python, 0.0 / 0.0 will cause ZeroDivisionError.
Unlike floating point numbers, ** integers have no way to store infinity or other anomalous values in the usual way. ** This means that, to be precise, the bit representation of such a number is not assigned in the bit representation of an integer. As a result, integer division by zero is undefined in ** C and C ++, with exceptions in most other languages. ** **
For example, consider the following C ++ code.
zero_division1.cpp
// C++
#include <iostream>
int main() {
int a = 1/0;
std::cout << "1/0 = " << a << "\n";
return 0;
}
It's out-of-the-box code, but if you compile it, you'll get the following warning, for example: For GCC.
zero_division1.cpp: In function 'int main()':
zero_division1.cpp:5:12: warning: division by zero [-Wdiv-by-zero]
int a = 1/0;
~^~
For Clang.
zero_division1.cpp:5:12: warning: division by zero is undefined [-Wdivision-by-zero]
int a = 1/0;
^~
1 warning generated.
In the case of GCC, if you execute the special warning without any notice,
$ ./a.out
Floating point exception: 8
And run-time exceptions, For 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
And every time it is executed, a random value is output. I think this depends on the environment, but it's undefined behavior anyway, so I can't complain no matter what happens.
Other slightly more secure (?) Languages will raise an exception, for example Java:
ZeroDivision.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)
Let's get back to the floating point NaN. NaN has the following properties.
-** The four arithmetic operations with NaN are always NaN **
a + b
a - b
a * b
a / b
Always returns NaN when a or b is NaN. For Java, it is as follows.
Nan2.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
So, when doing numerical calculations, it often happens that NaN generated at a certain point is transmitted in a blink of an eye and the output data is filled with NaN in the meantime.
-** Comparison operations with NaN other than! = Are always false,! = Are always true **
a == b
a > b
a < b
a >= b
a <= b
Always returns false when a or b is NaN, and also
a != b
Always returns true when a or b is NaN. For Java, it is as follows.
Nan3.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
In particular, ** NaN == NaN is false! ** ** So a common mistake when checking if x is NaN is
if (x == Double.NaN) {
// Do something...
} else {
// Do something...
}
Can't detect NaN with the code (the else is always executed). To be correct, most languages have their own methods such as isnan (x) and Double.isNaN (x).
if (Double.isNaN(x)) {
// Do something...
} else {
// Do something...
}
Let's write like this. Also, taking advantage of the above properties
if (x != x) {
// Do something...
} else {
// Do something...
}
Can also detect NaN. In fact, JavaScript before ECMAScript 6 used a similar method (About NaN determination, NaN is Not a Number but Number type. The story is). Also note that if a or b can be NaN, then a> b and! (A <= b) are not equivalent.
Based on the above, I would like to introduce one story that I actually experienced in the numerical calculation class at the university.
One of the algorithms for solving simultaneous equations is the Gauss-Seidel method (1 Solving simultaneous equations (iteration method). ),, [Simultaneous linear equations: Iterative method --PukiWiki for PBCG Lab](http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%CF%A2 % CE% A91% BC% A1% CA% FD% C4% F8% BC% B0% A1% A7% C8% BF% C9% FC% CB% A1)). Let's take a look at the following two implementations in C ++, leaving the details to the link above.
--Implementation 1
gauss_seidel1.cpp
#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;
}
--Implementation 2
gauss_seidel2.cpp
#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;
}
You don't have to worry about the details. What I would like you to see is the treatment of d. This d represents the "distance (to be exact, the L1 norm)" when the vector of the previous step is x1 and the vector of the current step is x2. And the former implementation continues iterating while this d is greater than epsilon, and the latter conversely stops iterating when d is less than or equal to epsilon, and in the end, which one is doing? Looks almost the same. In fact, both implementations
\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)
Solution of
\bf{x} =
\left(
\begin{array}{c}
-1 \\
1 \\
2
\end{array}
\right)
Is given correctly. However, this Gauss-Seidel method has its drawbacks. that is,
\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)
Even for simultaneous equations that have unique solutions, the solution may not be found unless the coefficient matrix is diagonally dominant. Specifically, the x1 [i] and x2 [i] values in the above code diverge, eventually becoming NaN. Then, since x1 and x2 are used to find x1 and x2, once NaN is generated from the above property 1., it remains NaN after that. Therefore, d obtained from x1 and x2 will always be NaN. Therefore, from property 2, both d> epsilon and d <= epsilon remain false.
......
Have you noticed it yet? If you can't find a solution well and it diverges, you can get out of the loop in implementation 1, but you end up in an infinite loop in implementation 2. As you can see, even if the algorithm is correct, it may not work depending on the value you enter, so if there is a risk of doing so, you need to take appropriate measures (in the above case, set the maximum number of iterations). , Throw an exception if it doesn't work).