Floating-Point Numbers: Format, Roundoff

• a 1 in the 2's place (2=2¹)
• a 0 in the 1's place (1=2⁰)
• a 0 in the (beyond the "binary point") 1/2's place (1/2=2^-1),
• a 1 in the 1/4's place (1/4=2^-2), and
• a 1 in the 1/8's place (1/8=2^-3)

Normalized Numbers

• You only need one bit to represent the sign--plus or minus.
• The exponent's just an integer, so you can store it as an integer.
• The 3.785746 part, called the "mantissa" or just "fraction" part, can be stored as the integer 3785746 (at least as long as you can figure out where the decimal point goes!)

1. There's exactly one digit to the left of the decimal point.
2. And that digit ain't zero.

Roundoff in Arithmetic

float f=0.73;
while (1) {
	volatile float g=f+1;
	if (g==f) {
		printf("f+1 == f  at f=%.3f, or 2^%.3f\n",
			f,log(f)/log(2.0));
		return 0;
	}
	else f=g;
}

double lil=1.0;
double big=pow(2.0,64);
printf(" big+(lil+lil) -big = %.0f\n", big+(lil+lil) -big);
printf("(big+lil)+lil  -big = %.0f\n",(big+lil)+lil  -big);

float gnats=1.0;
volatile float windshield=1<<24;
float orig=windshield;
for (int i=0;i<1000;i++)
	windshield += gnats;

if (windshield==orig) std::cout<<"You puny bugs can't harm me!\n";
else std::cout<<"Gnats added "<<windshield-orig<<" to the windshield\n";

float gnats=1.0;
volatile float windshield=1<<24;
float orig=windshield;
volatile float gnatcup=0.0;
for (int i=0;i<1000;i++)
	gnatcup += gnats;
windshield+=gnatcup; /* add all gnats to the windshield at once */

if (windshield==orig) std::cout<<"You puny bugs can't harm me!\n";
else std::cout<<"Gnats added "<<windshield-orig<<" to the windshield\n";

Roundoff in Representation

for (int i=1;i<1000000000;i*=10) {
	double mul01=i*0.1;
	double div10=i/10.0;
	double diff=mul01-div10;
	std::cout<<"i="<<i<<"  diff="<<diff<<"\n";
}

i=1  diff=5.54976e-18
i=10  diff=5.55112e-17
i=100  diff=5.55112e-16
i=1000  diff=5.55112e-15
i=10000  diff=5.55112e-14
i=100000  diff=5.55112e-13
i=1000000  diff=5.54934e-12
i=10000000  diff=5.5536e-11
i=100000000  diff=5.54792e-10
Program complete.  Return 0 (0x0)

Roundoff Taking Over Control

for (double k=0.0;k<1.0;k+=1.0/6.0) {
	printf("k=%a (about %.15f)\n",k,k);
}

(executable NetRun link)

• Don't use floating-point as your loop variable.  Loop over an integer i (without roundoff), and divide by six to get k.  This is the recommended approach if you care about the exact number of times around the loop.
  
• Or you could adjust the loop termination condition so it's "k<1.0-0.00001", where the "0.00001" provides some safety margin for roundoff.  This sort of "epsilon" value is common along floating-point boundaries, although too small and you can still get roundoff, and too big and you've screwed up the computation.
  
• Or you could use a lower-precision comparison, like "(float)k<1.0f".  This also provides roundoff margin, because the comparison is taking place at the lower "float" precision.
  

Bits in a Floating-Point Number

void print_bits(float f) {
	int i=*reinterpret_cast<int *>(&f); /* read bits with "pointer shuffle" */
	std::cout<<" float "<<std::setw(10)<<f<<" = ";
	for (int bit=31;bit>=0;bit--) {
		if (i&(1<<bit)) std::cout<<"1"; else std::cout<<"0";
		if (bit==31) std::cout<<" ";
		if (bit==23) std::cout<<" (implicit 1).";
	}
	std::cout<<std::endl;
}

int foo(void) {
	print_bits(0.0);
	print_bits(-1.0);
	print_bits(1.0);
	print_bits(2.0);
	print_bits(4.0);
	print_bits(8.0);
	print_bits(1.125);
	print_bits(1.25);
	print_bits(1.5);
	print_bits(1+1.0/10);
	return sizeof(float);
}

(Try this in NetRun now!)

/* IEEE floating-point number's bits:  sign  exponent   mantissa */
struct float_bits {
	unsigned int fraction:23; /**< Value is binary 1.fraction ("mantissa") */
	unsigned int exp:8; /**< Value is 2^(exp-127) */
	unsigned int sign:1; /**< 0 for positive, 1 for negative */
};

/* A union is a struct where all the fields *overlap* each other */
union float_dissector {
	float f;
	float_bits b;
};

float_dissector s;
s.f=8.0;
std::cout<<s.f<<"= sign "<<s.b.sign<<" exp "<<s.b.exp<<"  fract "<<s.b.fraction<<"\n";
return 0;

• Floats can store fractional numbers.
  
• Floats never overflow; they hit "infinity" as explored below.
  
• "double" has more bits than "int" (but less than "long").

Normal (non-Weird) Floats

Weird: Zeros and Denormals

Weird: Infinity

float z=0.0;
float f=1.0/z;
std::cout<<f<<"\n";
return (int)f;

(Try this in NetRun now!)

Weird: NaN

• 0.0/0.0 (neither zero nor infinity, because we'd want (x/x)==1.0; but not 1.0 either, because we'd want (2*x)/x==2.0...)
• infinity-infinity (might cancel out to anything)
• infinity*0

float f=sqrt(-1.0);
std::cout<<f<<"\n";
return (int)f;

(Try this in NetRun now!)

Bonus: Performance impact of special values

enum {n_vals=1000};
double vals[n_vals];

int average_vals(void) {
	for (int i=0;i<n_vals-1;i++) 
		vals[i]=0.5*(vals[i]+vals[i+1]);
	return 0;
}

int foo(void) {
	int i;
	for (i=0;i<n_vals;i++) vals[i]=0.0; 
	printf(" Zeros: %.3f ns/float\n",time_function(average_vals)/n_vals*1.0e9);
	for (i=0;i<n_vals;i++) vals[i]=1.0; 
	printf(" Ones: %.3f ns/float\n",time_function(average_vals)/n_vals*1.0e9);
	for (i=0;i<n_vals;i++) vals[i]=1.0e-310; 
	printf(" Denorm: %.3f ns/float\n",time_function(average_vals)/n_vals*1.0e9);
	float x=0.0;
	for (i=0;i<n_vals;i++) vals[i]=1.0/x; 
	printf(" Inf: %.3f ns/float\n",time_function(average_vals)/n_vals*1.0e9);
	for (i=0;i<n_vals;i++) vals[i]=x/x; 
	printf(" NaN: %.3f ns/float\n",time_function(average_vals)/n_vals*1.0e9);
	return 0;
}

Bonus: Arithmetic Tables for Special Floating-Point Numbers

Addition

Subtraction

Multiplication

Division

Equality

Less-Than