FIXME. This has not been presented to the students and has several problems. I'll come back to this later.
The smallest number representable by a 64-bit real is 2.2250738585072020*10-308. Its value is mainly limited by the number of digits you can store in the exponent (multiplied by a normalised mantissa whose value has the range 1.0..9.999).
Another small number of interest to computer programmers is the value of the smallest change you can make to a number (let's say number==1.0). This is determined by the number of bits in the mantissa. If the 23 bit mantissa for 1.0 on a 32 bit machine (being normalised) is 000000000000000000000000, then the smallest change will be to 000000000000000000000001 (only the last bit changes). The number represented by this difference is called the machine's epsilon (http://en.wikipedia.org/wiki/Machine_epsilon). The number is called the "machine's epsilon" because before IEEE 754, every machine had a different way of representing reals, and each machine would have its own epsilon. Now with most machines (and except for specialised embedded machines) being IEEE 754 compliant, the machine's epsilon will be the IEEE 754 epsilon. Let's ask our machine what its epsilon is.
>>> epsilon = 1 #has a simple representation in binary >>> while (epsilon+1>1): ... epsilon = 0.5*epsilon #still has a simple representation in binary ... print epsilon ... 0.5 0.25 0.125 . . 7.1054273576e-15 3.5527136788e-15 1.7763568394e-15 8.881784197e-16 4.4408920985e-16 2.22044604925e-16 1.11022302463e-16 |
If to the computer, 1+epsilon==1, the loop exits. The machine's epsilon for 64-bit IEEE 754 reals is 1.11022302463e-16 (for comparison, the smallest representable real is 2.2250738585072020*10-308). (note: epsilon is a lot bigger than the smallest number representable by the machine's reals. Do you understand the difference between these two numbers?)
What this section is saying is that the computer can't differentiate reals that are more closely spaced than the spacing represented by the last bit in the mantissa.
After shifting x=1/5 one place to the right (in the previous section) to give the result x=1/10, note that the right hand "1" falls off the end and is lost. In the following sequence of calculations (right shift, losing the last 1, then left shift, which will insert a 0, in the place which previously had a 1)
x=0.2 y=x/2 #y=0.1 z=y*2 #z=0.2 |
would the resulting z have the same value as the original x? If you were using fixed point operations, you would get a different value. However a computer uses a mantissa and an exponent to represent a real, so dividing by a power of 2 just changes the exponent and no precision is lost. The mantissa will be unchanged following the above operations, and z would have the same value as x.
>>> x=0.2 >>> y=x/2 >>> z=y*2 >>> print z-x 0.0 |
Taking this example to the extreme...
>>> n_128=2**128 >>> print n_128 340282366920938463463374607431768211456 >>> x=0.1 >>> y=x/n_128 >>> z=y*n_128 >>> print "%30.30g, %30.30g, %30.30g" %(x,y,z) 0.100000000000000005551115123126, 2.93873587705571893305445304302e-40, 0.100000000000000005551115123126 |
What if you did an operation that changed the mantissa (ignoring whether or not the exponent was changed)? We need a number that doesn't have a finite representation in binary.
>>> x=0.1 #x doesn't have a finite represenation in binary >>> n_prime=149.0 #1/n_prime doesn't have finite representation in binary either >>> y=x/n_prime >>> z=y*n_prime >>> print "%30.30g, %30.30g, %30.30g, %30.30g " %(x,y,z,z-x) 0.100000000000000005551115123126, 0.000671140939597315508424735241988, 0.100000000000000005551115123126, 0 |
I would have thought that z and x would be different, but you get back the same number. I puzzled about this, till a co-worker Ed Anderson suggested the following as a method to find out what was going on: try a whole bunch of numbers. In particular, try some irrational numbers, like square roots (irrational numbers cannot be represented by a finite number of digits, in any base system). All square roots, except integers, are irrational (for the range 1..10, only 1,4 and 9 would have rational square roots). So use all the square roots from 1..some_large_number for both x and the divisor/multiplier. Ed's initial test was a 100x100 grid. Here's Ed's code for a 64*64 grid. There is a blank when the numbers are identical, and a letter/number if they were different.
#! /usr/bin/python
# fp_ed.py
# Joseph Mack (C) 2008, GPL v3.
# With suggestions from Ed Anderson
# takes a number, divides it by a 2nd, stores the result, then multiplies the result the 2nd number.
# in a perfect world, the result should be the original answer.
# due to limited precision of computer reals, sometimes a different answer will be returned.
from math import sqrt;
number_that_are_different=0;
end = 100;
for j in range(1,end):
output_string = "" # can't output strings without \n in python and can't supress a blank either
output_string += str('%3d ' %j); #put blank after string
for i in range(1,end):
x = sqrt(i);
y = x/j;
z = y*j;
if (x-z == 0.0):
output_string += " "
else:
#print x-z ;
output_string += "x";
number_that_are_different += 1;
print output_string
print number_that_are_different
#- fp_ed.py -----------------------
|
Here's the output
1 2 3 4 5 6
1
2
3 x2 33 3 w 3 3
4
5 1 2 2 333 w 3 3
6 x2 33 3 w 3 3
7 x w w3 w 3
8
9 w w
10 1 2 2 333 w 3 3
11 x 2 w 3w3w ww w 3w 3 3
12 x2 33 3 w 3 3
13 xx ww w w w w 3w
14 x w w3 w 3
15 x w3w
16
17 w
18 w w
19 x w w ww
20 1 2 2 333 w 3 3
21 w 3 w 3 3 3 3
22 x 2 w 3w3w ww w 3w 3 3
23 x 3 3 3 w www3
24 x2 33 3 w 3 3
25 x x w w w w
26 xx ww w w w w 3w
27 xx w ww w w3 w
28 x w w3 w 3
29 2x 33w w3w
30 x w3w
31 w
32
33
34 w
35 2 3 3 3w 3
36 w w
37 y 2 x 33 w w w
38 x w w ww
39 2 33 3 33 w
40 1 2 2 333 w 3 3
41 x 3 3 3 w w 3
42 w 3 w 3 3 3 3
43 1y 2 x2x 2 w 3 w w 3 ww w 3 w 33
44 x 2 w 3w3w ww w 3w 3 3
45 3 w w 3 3
46 x 3 3 3 w www3
47 2x2 2 w333 w 3 w33 w w
48 x2 33 3 w 3 3
49 0 1 2 w w 3 3 3w
50 x x w w w w
51 2 x 3ww33 w 3 w
52 xx ww w w w w 3w
53 22 3 3 3 3 3 3
54 xx w ww w w3 w
55 y xx w3 3w 3 3 w33
56 x w w3 w 3
57 2 3ww3 w
58 2x 33w w3w
59 2x 33w3 w
60 x w3w
61 w w
62 w
63
|
What if we do the multiplication first
1 2 3 4 5 6
1
2
3 y x x w w 3w 3 3 w w
4
5 2 x 3 3w w ww 3
6 y x x w w 3w 3 3 w w
7 x w w3 3
8
9 x x w 3w w3 3
10 2 x 3 3w w ww 3
11 3 w 3 w w3
12 y x x w w 3w 3 3 w w
13 y x 2 2 3 3w w 3 3 3 w 3
14 x w w3 3
15 x w 3 w
16
17 w3 w3
18 x x w 3w w3 3
19 22 w 3 3w 3w3 3 w3
20 2 x 3 3w w ww 3
21 y 2 x2 3w3 w 3 w 3 3w 3
22 3 w 3 w w3
23 2 2 33 3 3 3 33
24 y x x w w 3w 3 3 w w
25 ww ww w
26 y x 2 2 3 3w w 3 3 3 w 3
27 1 2 2 w w 3 w 3 3 w 3 3
28 x w w3 3
29 w w w3
30 x w 3 w
31 3 3
32
33 w3
34 w3 w3
35 3 3
36 x x w 3w w3 3
37 2 3w3 w 3
38 22 w 3 3w 3w3 3 w3
39 y xx 2 ww3 w w ww3 3 3
40 2 x 3 3w w ww 3
41 x x 3 w 333 w 3 www
42 y 2 x2 3w3 w 3 w 3 3w 3
43 2 3 ww 3 w 3
44 3 w 3 w w3
45 y 2xx w w 3 ww w w
46 2 2 33 3 3 3 33
47 y x 2x w 33 3 3www 3
48 y x x w w 3w 3 3 w w
49 x 3 3 w ww 3w 3
50 ww ww w
51 1 2 22 3 3 33 3 w
52 y x 2 2 3 3w w 3 3 3 w 3
53 2x 33 w 3 3 w
54 1 2 2 w w 3 w 3 3 w 3 3
55 x w w 3 3 w 3
56 x w w3 3
57 3 w 3 w
58 w w w3
59 x 3 w 3 w
60 x w 3 w
61 w
62 3 3
63 3
|
The figures are different. It seems that multiplication and division isn't always commutative.
Surprisingly (to me) there is a pattern, but then I didn't understand what was going on so any result would have been a surprise. "x"s appear in blocks on boundaries which are multiples of 8,16,32 and 64 (increasing the value of end shows that the block pattern continues).
With a slight modification of the code, the differences, if not 0.0, were all seen to be (1.11022302463e-16)*2n.
Is there anything in particular about 1.11022302463e-16? It turns out that it's the machine's epsilon (http://en.wikipedia.org/wiki/Machine_epsilon) (also see machine's epsilon).
What is 1.11022302463e-16 in binary?
dennis:~# echo "obase=2;1.11022302463*10^-16" | bc -l .0000000000000000000000000000000000000000000000000000011111111111111 |
If you round this up by one bit, it's 2^-53. A 64 bit double has a 53 bit mantissa, indicating that the epsilon calculation is using a 64 bit real.
From the wiki entry, a single precision (32bit) real has a 24 bit mantissa. 1.0 then has an exponent of 0 and mantissa of 1.000000000000000000000002. The next biggest number has an exponent of 0 and a mantissa of 1.000000000000000000000012. The difference is 0.000000000000000000000012 or 2-23.
Using IEEE 754 Converter (http://www.h-schmidt.net/FloatApplet/IEEE754.html), to do the conversion to a 32-bit real, I found
1.11022302463e-16=0x2500000=00100101000000000000000000000000. |
Breaking the binary representation into sign, exponent and mantissa (with the implied 1 for a normalised mantissa) gives
1.11022302463e-16 = 0 01001010 (1)00000000000000000000000
decimal s exp mantissa
|
Calculating this all gives: The sign bit: (bit 32, the leftmost bit) is 0 (positive).
The exponent: the IEEE-754 real math convention doesn't use the signed int representation. a 32 bit real, the exponent is a 9 bit integer (here 001001010). Floating point math is usually done on floating point hardware (or emulated on the the integer registers of the CPU, but the incremental cost of adding a floating point processor/unit (FPU) is small, so all x86 chips, starting with the 486, have a floating point processor built-in to the CPU chip). Since real math doesn't use the CPU's registers (each composed of 2,4,8 or 16 bytes), but instead uses the math co-processor's registers, then real math doesn't have to use the CPU's conventions (e.g. a signed int) to process a 9 bit int.
For 32 bit floating point, the convention is to subtract 127 from the value of the exponent. is is called "excess-127" (excess-1023 for 64 bit) and allows the exponent and the mantissa, when taken together, to vary monotonically (see IEEE-754 References http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html). (What this gets you, I don't know.) The 8-bit scheme 8-bit set of normalised reals I made up for this class is not monotonic in this respect (see the table and the figure, where there are jumps in value every 16 numbers).
In the current example, the exponent is -53 (from the next 9 bits). To get the actual value of the exponent, you subtract 127.
echo "obase=10; ibase=2; 001001010" | bc 74 then subtract 127 74-127=-53 |
The FPU int then is an unsigned int, from which is subtracted 127, so that -ve numbers are possible (it's called signed magnitude). Here are a few comparison numbers (assuming an 8 bit int):
bits value, signed int value, unsigned int value, signed magnitude 00000000 0 0 -127 00000001 1 1 -126 01111111 127 127 0 10000000 -127 128 1 11111111 -1 255 128 |
Since all of the evaluations can done with the same speed, with hardware designed to do these calculations, there is nothing to pick between them. Why IEEE-754 uses their particular representation of the exponent, rather than the signed int representation, I don't know.
The mantissa is 1.0 (after adding the implied 1 at the beginning of 23 zeroes).
Multiplying this all out gives the original number.
echo "1*2^-53" | bc -l .00000000000000011102 #just checking that I can count 0s echo "1*2^-53/10^-16" | bc -l 1.11020000000000000000 |
The following table shows the number of times an entry occured for a grid size. The 2nd entry says that in an 16*16 grid, that 206 numbers were not changed and that -(1.11022302463e-16)*23 was returned 9 times, +(1.11022302463e-16)*2*2 was returned 2 times and that +(1.11022302463e-16)*2*3 was returned 8 times. In the 16*16 grid 15*15=225 numbers are tested (=9+206+2+8).
grid size\n -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 8*8 0 0 0 0 0 0 0 0 0 0 0 0 48 0 1 0 0 0 0 0 0 0 0 0 0 16*16 0 0 0 0 0 0 0 0 0 9 0 0 206 0 2 8 0 0 0 0 0 0 0 0 0 32*32 0 0 0 0 0 0 0 0 6 22 0 0 911 0 3 13 6 0 0 0 0 0 0 0 0 64*64 0 0 0 0 0 0 0 0 163 44 3 0 3566 1 6 33 153 0 0 0 0 0 0 0 0 128*128 0 0 0 0 0 0 0 99 346 82 10 0 15062 4 14 80 349 83 0 0 0 0 0 0 0 256*256 0 0 0 0 0 0 0 3085 723 163 24 0 56990 14 39 174 750 3063 0 0 0 0 0 0 0 512*512 0 0 0 0 0 0 1679 6215 1505 355 56 0 241228 35 87 381 1570 6360 1650 0 0 0 0 0 0 1024*1024 0 0 0 0 0 0 50447 12720 3101 752 121 0 911185 91 196 799 3209 12807 51101 0 0 0 0 0 0 2048*2048 0 0 0 0 0 26941 101676 25687 6251 1500 250 0 3863963 230 452 1669 6449 25610 102498 27033 0 0 0 0 0 4096*4096 0 0 0 0 0 822872 204267 51363 12576 3001 525 0 14575322 504 969 3366 12919 51254 205106 824981 0 0 0 0 0 8192*8192 0 0 0 0 432519 1647991 409597 102641 25183 5958 1057 0 61830620 1105 2100 6885 25992 103109 411309 1652318 434097 0 0 0 0 16384*16384 0 0 0 0 13222697 3300406 821156 205169 50246 11827 2083 0 233145705 2345 4366 13903 52320 206826 824241 3306218 13233181 0 0 0 0 32768*32768 0 0 0 6931293 26447787 6604323 1645180 409899 100085 23405 4099 0 989292832 4872 9022 28195 105323 414810 1651250 6614158 26462945 6926811 0 0 0 65536*65536 0 0 0 211603215 52895451 13214901 3295139 820314 200000 46646 8163 0 3730616481 9841 18209 56688 211726 830962 3306011 13229067 52912230 211561181 0 0 0 131072*131072 262144*262144 524288*524288 1048576*1048576 8*8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.980 0.0 0.020 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 16*16 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.040 0.0 0.0 0.916 0.0 0.009 0.036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 32*32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.006 0.023 0.0 0.0 0.948 0.0 0.003 0.014 0.006 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 64*64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.041 0.011 0.001 0.0 0.898 0.0 0.002 0.008 0.039 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 128*128 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.006 0.021 0.005 0.001 0.0 0.934 0.0 0.001 0.005 0.022 0.005 0.0 0.0 0.0 0.0 0.0 0.0 0.0 512*512 0.0 0.0 0.0 0.0 0.0 0.0 0.006 0.024 0.006 0.001 0.0 0.0 0.924 0.0 0.0 0.001 0.006 0.024 0.006 0.0 0.0 0.0 0.0 0.0 0.0 1024*1024 0.0 0.0 0.0 0.0 0.0 0.0 0.048 0.012 0.003 0.001 0.0 0.0 0.871 0.0 0.0 0.001 0.003 0.012 0.049 0.0 0.0 0.0 0.0 0.0 0.0 2048*2048 0.0 0.0 0.0 0.0 0.0 0.006 0.024 0.006 0.001 0.0 0.0 0.0 0.922 0.0 0.0 0.0 0.002 0.006 0.024 0.006 0.0 0.0 0.0 0.0 0.0 4096*4096 0.0 0.0 0.0 0.0 0.0 0.049 0.012 0.003 0.001 0.0 0.0 0.0 0.869 0.0 0.0 0.0 0.001 0.003 0.012 0.049 0.0 0.0 0.0 0.0 0.0 8192*8192 0.0 0.0 0.0 0.0 0.006 0.025 0.006 0.002 0.0 0.0 0.0 0.0 0.922 0.0 0.0 0.0 0.0 0.002 0.006 0.025 0.006 0.0 0.0 0.0 0.0 16384*16384 0.0 0.0 0.0 0.0 0.049 0.012 0.003 0.001 0.0 0.0 0.0 0.0 0.869 0.0 0.0 0.0 0.0 0.001 0.003 0.012 0.049 0.0 0.0 0.0 0.0 32768*32768 0.0 0.0 0.0 0.006 0.025 0.006 0.002 0.0 0.0 0.0 0.0 0.0 0.921 0.0 0.0 0.0 0.0 0.0 0.002 0.006 0.025 0.006 0.0 0.0 0.0 65536*65536 0.0 0.0 0.0 0.049 0.012 0.003 0.001 0.0 0.0 0.0 0.0 0.0 0.869 0.0 0.0 0.0 0.0 0.0 0.001 0.003 0.012 0.049 0.0 0.0 0.0 131072*131072 262144*262144 524288*524288 1048576*1048576 |
After the multiplication/division, how many numbers are returned with a different value? For any particular j you have a 50% chance of losing a 0 off the right end, in which case you'll get back the original answer. So one possibility is 50% will be different. How about the divisor?; you only have a 50% chance of getting that right when you drop off extra digits. With two ways of getting it wrong, you only have 25% chance of getting it right, so 75% will be different. Here's the results.
grid_size number_mul/div number_changed %changed largest_difference 10*10 10^2 2 2 + 4 * 1.11e-16 100*100 10^4 713 7 +- 16 * 1.11e-16 1000*1000 10^6 127980 13 +- 32 * 1.11e-16 10000*10000 10^8 8060303 8 +- 128 * 1.11e-16 100000*100000 10^10 927798399 9 +- 512 * 1.11e-16 |
Let's look at this example with the IEEE-754 converter referred to above
>>> sqrt(14) 3.7416573867739413 >>> x=sqrt(14) >>> y=x/3 >>> z=y*3 >>> z-x -4.4408920985006262e-16 |
#from IEEE 754 converter
x = 3.7416573867739413 = 2 * 1.8708287 = 0 10000000 (1) 11011110111011101010001
3.0 = 2 * 1.5 = 0 10000000 (1) 10000000000000000000000
#from python
y = x/3.0 = 1.247219128924647 = 0 01111111 (1) 00111111010010011100000
#from bc
echo "obase=10;ibase=2; 1.11011110111011101010001/1.10000000000000000000000" | bc -l
1.24721916516621907552
echo "obase=2;ibase=2; 1.11011110111011101010001/1.10000000000000000000000" | bc -l
1.0011111101001001110000010101010101010101010101010101010101010101010
#from python
z = 3.7416573867739409
from bc
echo "obase=10;ibase=2; 1.0011111101001001110000010101010101010101010101010101010101010101010 *11" | bc -l
3.7416574954986572265489474728439311945749068399891257286071777343750
echo "obase=2;ibase=2; 1.0011111101001001110000010101010101010101010101010101010101010101010 *11" | bc -l
11.101111011101110101000011111111111111111111111111111111111111111111
The result after right shifting by one (to give an exponent of 1), and dropping the implied leading 1 from the mantissa give
a mantissa of
11011110111011101010000 #new mantissa (after dropping extra digits)
11011110111011101010001 #original mantissa
1 #difference
since the exponent is 1 (a factor of 2) then the difference is
|