Remember this image from the lecture:
These values can be played with directly inside a debugger, such as
gdb:
(gdb) list
3 main()
4 {
5 float a = 1.1, b = 1.2;
6 int i;
7
8 for ( i = 0 ; i < 100000000 ; i++ ) {
9 a = expf(1.0);
10 a*=b;
11 b*=a;
12 a*=b;
(gdb) ptype a
type = float
(gdb) whatis a
type = float
(gdb) x/wt &a
0xbf98f228: 01000000010010010000111111010000
(gdb) set var a=3.1415926
(gdb) x/wt &a
0xbf98f228: 01000000010010010000111111011010
(gdb) p a
$3 = 3.1415925
(gdb)
- 00000000000000000000000000000000
0.0
- 00111111100000000000000000000000
1.0
sign (符号部) (one bit): 0
significand (仮数部) (twenty-three bits): 0b00000000000000000000000 = 1.0....(binary) (ひとつのビットはかくれていると覚えて)(remember the hidden bit)
exponent (指数部) (eight bits): 0b01111111 (バイアスは127で覚えて)(remember the bias of 127)
- 10111111100000000000000000000000
-1.0
sign (符号部) (one bit): 1
significand (仮数部) (twenty-three bits): 0b00000000000000000000000 = 1.0....(binary) (ひとつのビットはかくれていると覚えて)(remember the hidden bit)
exponent (指数部) (eight bits): 0b01111111 (バイアスは127で覚えて)(remember the bias of 127)
- 00111111000000000000000000000000
0.5
sign (符号部) (one bit): 1
significand (仮数部) (twenty-three bits): 0b00000000000000000000000 = 1.0....(binary) (ひとつのビットはかくれていると覚えて)(remember the hidden bit)
exponent (指数部) (eight bits): 0b01111110 (バイアスは 127で覚えて)(remember the bias of 127)
- 01000000000000000000000000000000
2.0
sign (符号部) (one bit): 0
significand (仮数部) (twenty-three bits): 0b10000000000000000000000=1.0...(binary) (ひとつのビットはかくれていると覚えて)(remember the hidden bit)
exponent (指数部) (eight bits): 0b10000000 (バイアスは127で覚えて)(remember the bias of 127)
- 01000000010000000000000000000000
3.0
sign (符号部) (one bit): 0
significand (仮数部) (twenty-three bits): 0b00000000000000000000000=1.10...(binary) (ひとつのビットはかくれていると覚えて)(remember the hidden bit)
exponent (指数部) (eight bits): 0b10000000 (バイアスは127で覚えて)(remember the bias of 127)
- 01000000010000000000000000000001
3.0+2^(-23)
(gdb) set var b = 3.0000002 (gdb) x/wt &b 0xbf98f22c: 01000000010000000000000000000001
- 01000000010010010000111111011010
3.1415926
Here is an example for integer addition:
[rdv@dhcp-143-221 architecture]$ cat int.c
main()
{
int a = 7, b = 2, c = 5, i;
for ( i = 0 ; i < 100000000 ; i++ ) {
c += a;
a += b;
b += c;
a += a;
a += a;
b += b;
c += a;
a += b;
b += a;
c += a;
c += a;
a += b;
b += c;
a += a;
a += a;
b += b;
c += a;
a += b;
b += a;
c += a;
a = 7;
b = 2;
c = 5;
}
}
[rdv@dhcp-143-221 architecture]$ gcc -o int int.c
[rdv@dhcp-143-221 architecture]$ time ./int
real 0m3.839s
user 0m3.828s
sys 0m0.010s
[rdv@dhcp-143-221 architecture]$ time ./int
real 0m3.859s
user 0m3.846s
sys 0m0.011s
[rdv@dhcp-143-221 architecture]$ time ./int
real 0m3.815s
user 0m3.809s
sys 0m0.003s
The inner loop is ten instructions, and the loop is executed 100,000,000 times, so it runs 2 billion integer additions. 2E9 / 3.8 seconds = 525M integer additions/second. On a 2GHz processor, that is one result every 3.8 clock cycles.
Similar programs are easily written for the other operations.
[rdv@dhcp-143-221 architecture]$ cat fp.c
#include <math.h>
main()
{
float a = 1.1, b = 1.2;
int i;
for ( i = 0 ; i < 100000000 ; i++ ) {
a*=b;
b*=a;
a*=b;
b*=a;
a*=b;
b*=a;
a*=b;
b*=a;
a*=b;
b*=a;
a = 1.1; b = 1.2;
}
}
[rdv@dhcp-143-221 architecture]$ gcc -o fp fp.c
[rdv@dhcp-143-221 architecture]$ time ./fp
real 0m3.026s
user 0m3.016s
sys 0m0.009s
[rdv@dhcp-143-221 architecture]$ time ./fp
real 0m3.070s
user 0m3.059s
sys 0m0.011s
[rdv@dhcp-143-221 architecture]$ time ./fp
real 0m2.938s
user 0m2.936s
sys 0m0.002s
[rdv@dhcp-143-221 architecture]$ time ./fp
real 0m2.917s
user 0m2.907s
sys 0m0.010s
The inner loop is ten instructions, and the loop is executed 100,000,000 times, so it runs one billion floating point multiplies. One billion FP multiplies in 3 seconds is 3.3E8 FP multiplies/second, or one every 6.6 clock cycles with a 2.0GHz processor.
However, there are a bunch of data dependencies in that code: if we double the number of FP multiplies, the amount of time less than doubles.
[rdv@dhcp-143-221 architecture]$ cat fp2.c
#include <math.h>
main()
{
float a = 1.1, b = 1.2, c = 1.3, d = 1.4;
int i;
for ( i = 0 ; i < 100000000 ; i++ ) {
a*=b;
c*=d;
b*=a;
d*=c;
a*=b;
c*=d;
b*=a;
d*=c;
a*=b;
c*=d;
b*=a;
d*=c;
a*=b;
c*=d;
b*=a;
d*=c;
a*=b;
c*=d;
b*=a;
d*=c;
a = 1.1; b = 1.2; c = 1.3; d = 1.4;
}
}
[rdv@dhcp-143-221 architecture]$ gcc -o fp2 fp2.c
[rdv@dhcp-143-221 architecture]$ time ./fp2
real 0m5.073s
user 0m5.068s
sys 0m0.002s
[rdv@dhcp-143-221 architecture]$ time ./fp2
real 0m5.142s
user 0m5.139s
sys 0m0.001s
[rdv@dhcp-143-221 architecture]$ time ./fp2
real 0m5.085s
user 0m5.079s
sys 0m0.002s
In that example, 2 billion instructions are executed in 5 seconds, for 4E8 FP multiplies/second, one every 5 clock cycles.
Notice the line from the example above:
(gdb) set var a=3.1415926 (gdb) x/wt &a 0xbf98f228: 01000000010010010000111111011010 (gdb) p a $3 = 3.1415925
The last digit is different. That suggests that there is a problem, but at the point it is not identified. It might be in the routine that converts floating point binary to human-readable ASCII.
You can see from the code and output below that it is possible for floating point operations to become inaccurate after only a few dozen divisions, or additions of small numbers to large ones.
[rdv@dhcp-143-221 architecture]$ cat double-double-double.cpp
// first, the standard libraries
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <set>
#include <map>
#include <cmath>
#include <iterator>
#include <vector>
#define FALSE 0
#define TRUE 1
main()
{
float f = 1.0;
double a = 1.0;
long double b = 1.0;
long double c = 1.0;
int i;
int reported, reported2;
std::cout.precision(50);
std::cout << "float (" << sizeof(f)*8 << "bits):" << std::endl;
i = 0; reported = FALSE; reported2 = FALSE;
while (1) {
float fprime = f/2;
// n.b.: where this reports will depend on the architecture, the compiler,
// and the compiler optimization switches! Might optimize to 1.0-f^2,
// in which case the question is, when does f^2 == 0
// but an expression similar to this likely occurs when manipulating
// fidelities in a quantum simulation, since the fidelities sum to 1.0.
if (!reported && (1.0-f)*f == f) {
std::cout << "After " << i << " iterations, ((1.0-f)*f == f) at f = " << f << std::endl;
reported = TRUE;
}
// try f^2 directly.
// Note that on many machines this will *never* trip because
// the CPU internal floating format is double, and we'll run out of precision
// on the divide by two first
if (!reported2 && f*f == 0) {
std::cout << "f^2 becomes zero at " << i << " iterations" << std::endl;
reported2 = TRUE;
}
if (f != 2*fprime) {
std::cout << "After " << i << " iterations: " << f << std::endl;
break;
}
f = fprime;
i++;
}
std::cout << "double (" << sizeof(a)*8 << "bits):" << std::endl;
i = 0; reported = FALSE; reported2 = FALSE;
while (1) {
double aprime = a/2;
if (!reported && (1.0-a)*a == a) {
std::cout << "After " << i << " iterations, ((1.0-a)*a == a) at a = " << a << std::endl;
reported = TRUE;
}
// try a^2 directly
if (!reported2 && a*a == 0) {
std::cout << "a^2 becomes zero at " << i << " iterations" << std::endl;
reported2 = TRUE;
}
if (a != 2*aprime) {
std::cout << "After " << i << " iterations: " << a << std::endl;
break;
}
a = aprime;
i++;
}
std::cout << "long double (" << sizeof(b)*8 << "bits):" << std::endl;
i = 0; reported = FALSE; reported2 = FALSE;
while (1) {
long double bprime = b/2;
if (!reported && (1.0-b)*b == b) {
std::cout << "After " << i << " iterations, ((1.0-b)*b == b) at b = " << b << std::endl;
reported = TRUE;
}
// try b^2 directly
if (!reported2 && b*b == 0) {
std::cout << "b^2 becomes zero at " << i << " iterations" << std::endl;
reported2 = TRUE;
}
if (b != 2*bprime) {
std::cout << "After " << i << " iterations: " << b << std::endl;
std::cout << b << std::endl;
break;
}
b = bprime;
i++;
}
std::cout << "long double (" << sizeof(c)*8 << "bits):" << std::endl;
i = 0; reported = FALSE; reported2 = FALSE;
while (1) {
long double cprime = c/2;
if (!reported && (1.0-c)*c == c) {
std::cout << "After " << i << " iterations, ((1.0-c)*c == c) at c = " << c << std::endl;
reported = TRUE;
}
// try c^2 directly
if (!reported2 && c*c == 0) {
std::cout << "c^2 becomes zero at " << i << " iterations" << std::endl;
reported2 = TRUE;
}
if (c != 2*cprime) {
std::cout << "After " << i << " iterations: " << c << std::endl;
std::cout << c << std::endl;
break;
}
c = cprime;
i++;
}
}
[rdv@dhcp-143-221 architecture]$ g++ -o double-double-double double-double-double.cpp
[rdv@dhcp-143-221 architecture]$ ./double-double-double
float (32bits):
After 65 iterations, ((1.0-f)*f == f) at f = 2.710505431213761085018632002174854278564453125e-20
After 149 iterations: 1.4012984643248170709237295832899161312802619418765e-45
double (64bits):
After 65 iterations, ((1.0-a)*a == a) at a = 2.710505431213761085018632002174854278564453125e-20
After 1074 iterations: 4.9406564584124654417656879286822137236505980261432e-324
long double (96bits):
After 65 iterations, ((1.0-b)*b == b) at b = 2.710505431213761085018632002174854278564453125e-20
b^2 becomes zero at 8223 iterations
After 16445 iterations: 3.6451995318824746025284059336194198163990508156936e-4951
3.6451995318824746025284059336194198163990508156936e-4951
long double (96bits):
After 65 iterations, ((1.0-c)*c == c) at c = 2.710505431213761085018632002174854278564453125e-20
c^2 becomes zero at 8223 iterations
After 16445 iterations: 3.6451995318824746025284059336194198163990508156936e-4951
3.6451995318824746025284059336194198163990508156936e-4951
