… and when does 35 – 35 not equal zero?
In the first case the answer is when one “35” is on an Excel spreadsheet and the other 35 is in VBA. The question arose from the ongoing continuous beam discussion with metroxx, who had, quite reasonably, generated a series of output points by calculating the length of one step (35/30) then adding this value to the previous length in 30 stages. This mysteriously caused an error in the function for no obvious reason.
Stepping through the code I found that the error was occurring when the last output point was compared to the beam length and found to be greater. The mysterious part was that both on the worksheet (formatted to 15 decimal places), and in the VBA Locals window, both values displayed as being exactly 35. Also entering the formula = A62 – 35 returned exactly zero.
To investigate what was going on I wrote the short VBA function shown below:
Function Diff(A As Double, B As Double) As Double Diff = A - B End Function
The function simply subtracts B from A, but using this function with the incrementally generated “35” as A and the “exact” 35 as B returned a value of 7.11E-15, even though stepping through the routine both variables displayed as exactly 35 in the locals window.
The answer to the second question is the same as the first (as would be expected), but also in some cases this behaviour can be seen on the spreadsheet with no VBA involved. As stated above, the formula = A62 – 35 returns exactly zero, but the formula = (A62 – 35) returns 7.11E-15, the same as the UDF! The screenshots below show more detail of this behaviour.
In the spreadsheet in column B I have generated 30 increments with the formula =(B7+$C$4), where C4 is =C2/C3, i.e. 35/30.
The value in column C are generated with =$C$2*A8/$C$3, where C2 is the beam length, A8 is the increment number, and C3 is the number of increments (30). Column D contains the the formula =B8 – C8, Column E: =Diff(B8, CB), and column F: =(B8 – C8).
Floating point value check - click to view full size
Error in Column D result
It can be seen that Column D has returned a difference of zero in all cases, even though the VBA function in Column E and the formula in Column F show a difference of up to 1.42E-14.
These differences are of course caused in part by the fact that all values are stored as binary floating point values, that cannot represent all decimal or fractional values exactly. This is described in several Microsoft documents (e.g. Understanding Floating Point Precision), which all claim that Excel follows the IEEE Standard for Binary Floating-Point Arithmetic. Clearly this is not the whole story though, since the example given generates errors up to an order of magnitude greater than the maximum difference between any exact value and the nearest floating point value, and also generates different results depending on how worksheet formulas are entered.
I would be interested if anyone has any more background on exactly how these things are handled in Excel, but for practical purposes I think the lessons are:
- When comparing non-integer values be aware that values that display as exactly equal may be stored as different values.
- When comparing doubles in VBA either round the values to a suitable precision, or check that the difference is less than some small value, rather than exactly zero.
- Calculations that involve the difference between two nearly equal values may give incorrect results. Consider carrying out this type of calculation entirely in VBA (or if necessary in a different language offering higher precision calculations).
Example added 27th Dec 2011:
Another example of the effect of brackets:
- Cell B4, =B2 + B3, displays as exactly 1, because the 1.6E-15 would be the 16th and 17th significant figure, but only 15 significant figures are displayed.
- Cell B5, =B4 – 1, displays as exactly 0, the contents of B4 being treated as the displayed value.
- Cell B6, =(B2 + B3), displays as exactly 1, the same as B4, as would be expected.
- But Cell B7, =B6-1, displays as 0.000000000000001554, using the decimal equivalent of the underlying value stored in B6, rather than the displayed value.
- Cell B8, =(B6-B4), displays exactly 0, using the displayed values in both cells
- Cell B9, = B7-B5, displays 1.554E-15 (i.e. the same as B7), so
((1+x))-(1+x) <> ((1+x)-1)-(1+x-1) !
Newton Excel Bach, not (just) an Excel Blog 
Hi all, here comes a short list of links that helped me to understand what IEEE 754 DP really means for computational arithmetics:
http://support.microsoft.com/kb/78113
here MS states where VBA (not XL!!) adheres to IEEE 754 and where not
http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
imho, good summary of all effects
http://www.cs.berkeley.edu/~wkahan/
homepage of Prof. Kahan, many papers available for download; among them this nice “fatwah” on Java’s IEEE implementation; imho, even VBA programmers can learn much from it
Click to access JAVAhurt.pdf
http://www.exploringbinary.com/
imho, a nice website on the subject
“staircase to hell” demonstrating one of XL’s dirty secrets, i.e., manipulating a wide range of numbers around the exact doubles
—
According to my experience, XL cannot be used to pre-test VBA code in detail because firstly, the results of even simple arithmetic operations may differ and secondly, even the standard mathematical functions sometimes yield different results in XL and VBA, resp. This holds true for XL10/VBA6/x32 as well as for XL14/VBA7/x64. Moreover, unlike C or Fortran VBA does not seem to take full advantage of the cache registers providing extended precision for intermediate results that all modern processors are equipped with. Therefore, I recommend to use VBA only to interface to routines written in a more suitable language like C or Fortran if accuracy really matters, for example when solving ODEs or PDEs, resp.
LikeLike
Agree it’s strange that parentheses around a formula can change the result?!
The differences in column F of the screenshot are exact binary powers: 2^-50, 2^-49, etc. as one might expect from the binary representations. This ties in with the fact that Excel truncates numbers beyond the fifteenth significant figure in calculations (http://en.wikipedia.org/wiki/Numeric_precision_in_Microsoft_Excel), whereas seventeen significant figures are needed to store numbers exactly in the xml file formats. In VBA results are displayed rounded to fifteen significant figures but are stored internally to greater precision.
LikeLike
After further investigation, it’s clear that formulas in column D return zeroes due to truncation but the differences in column F are generic IEEE double round-off errors. For instance this c code returns identical results:
double sum=0;
int i;
for(i=1; i<=30; i++)
{
sum+=35.0/30.0;
printf_s("%E\n",sum-i*35.0/30.0);
}
(The chart could also be converted to a step chart if desired.)
LikeLike
Ihm, the key needed to understand all those “strange” effects occuring with the DP numbers is to realise that the gap between two adjacent DP numbers doubles at each exact non-zero power of two. That is why in VBA if statements, e.g., you should only compare DP numbers to 0# and not to any other DP number, because in the latter case the interval yielding true/false is no more symmetric around your DP reference number.
Imho, the behaviour of the DP numbers can only be fully understood if we look at their binary representation and not at their (mostly) rounded decimal representation. If you know that 0.1 is internally approximated by the IEEE 754 DP binary number 3F B9 99 99 99 99 99 9A (little endian byte order), and that 0.5-0.4 is approximated by 3F B9 99 99 99 99 99 98, you will not be surprised that 2^55*(0.1-(0.5-0.4)) yields a result of 1, whereas 2^55*((0.5-0.4)-0.1) yields a result of -1, though mathematically, both should yield exactly 0. Afaik, forthcoming versions of xNumbers will provide functions that allow for conversions between the decimal and binary representation of a variable of type “Double” in VBA. That VBA lacks the support of subnormal numbers is a big difference between VBA and F, C,… Start with 2^(-1070) and successively divide by 2 in XL, VBA, C, F… and print out the results…
LikeLike
Sorry, I forgot to add the codes of intermediate results:
in XL: =0.1-(0.5-0.4) -> 00 00 00 00 00 00 00 00 -> 0
exactly 0, but this is one of XL’s dirty manipulations
in XL: =(0.1-(0.5-0.4)) -> 3C 80 00 00 00 00 00 00 -> 2.7756E-17
using VBA and Doug’s vdiff:
in XL: = vdiff(0.1;vdiff(0.5;0.4)) -> 3C 80 00 00 00 00 00 00 -> 2.7756E-17
Btw, the current version of xNumbers provides the functions vAdd, vSub, vMult, and vDiv. In order to test whether standard functions like the exponential yield the same result in XL as well as in VBA (using an obviously defined UDF vExp) one could iterate over random numbers in the cells of col. A and then use something like =vSub(exp(Annn);vExp(Annn)) to check for differences… =exp(Annn)-vExp(Annn) will always yield exactly 0 whereas =2^55*(exp(Annn)-vExp(Annn)) will not…
I wish everybody a Merry Christmas!
Georg
LikeLike
Hi Lori and Georg
Thanks for the discussion and interesting links. Some thoughts coming out of that:
– Firstly I was wrong about the errors being larger than would be given by the IEEE standard, I had forgotten to allow for the increasing magnitude of the base number.
– On the other hand the different results when enclosing the simple subtraction formula in brackets does seem surprising, and doesn’t seem to be explained by any of the documentation.
– Also I think that Microsoft are being a bit disingenuous in claiming that the loss of precision is all a result of following the standard, and that all other spreadsheets are the same. The standard does allow for 80 bit precision, and at least one spreadsheet (later DOS versions of Lotus 123) does/did provide that. I’m not sure why Lotus went back to 64 bit in their Windows versions. Excel compatability perhaps?
– That said, 15 significant figures is more than enough for almost all practical applications, even those involving inversion of matrices and the like. For instance the beam analysis program I’ve been working on recently agrees to 14 or 15 SF with the results of a commercial FEA program written in Fortran, using a completely different analysis method, which is at least 10 more SF that you really need!.
– Finally Georg – do you have any information on new versions of Xnumbers? I understood that development had stopped (at least by the original team).
LikeLike
Hi Doug,
this strange bracketing effect always occurs if a non-zero DP difference is tiny but still so small that XL would set it to exactly zero in order not to confuse people. Imho, the outer brackets make XL’s parser not manipulate the final result… Obtaining all the results by applying the v…-functions allows for correctly predicting what would happen in a VBA module.
Btw, some miracle has happened since September in my PCs: VBA6/VBA7 can now handle subnormal numbers. This can be demonstrated after defining a vPow as well as vStr UDF like this, for examle:
Function vPow#(ByVal myX#, ByVal myIPow%)
vPow = myX ^ myIPow
End Function
Function vStr$(ByVal myX#)
vStr = CStr(myX)
End Function
XL cell formula -> result of DBL2HEX UDF -> XL output
=vPow(2;-1022) -> 00 10 00 00 00 00 00 00 -> 2.2251E-308 (smallest normalised number)
=vPow(2;-1023) -> 00 08 00 00 00 00 00 00 -> 0.00000E-01
…
=vPow(2;-1026) -> 00 01 00 00 00 00 00 00 -> 0.00000E-01
=vPow(2;-1027) -> 00 00 80 00 00 00 00 00 -> 0.00000E+00
…
=vPow(2;-1074) -> 00 00 00 00 00 00 00 01 -> 0.00000E+00 (smallest subnormal number)
=vPow(2;-1075) -> 00 00 00 00 00 00 00 00 -> 0.00000E+00
=vStr(vPow(2;-1022)) -> 00 0F FF FF FF FF FF FD -> 2.2251E-308
=vStr(vPow(2;-1023)) -> 00 07 FF FF FF FF FF FF -> 1.1125E-308
…
=vStr(vPow(2;-1026)) -> 00 01 00 00 00 00 00 00 -> 1.3907E-309
=vStr(vPow(2;-1027)) -> 00 00 80 00 00 00 00 00 -> 6.9534E-310
…
=vStr(vPow(2;-1074)) -> 00 00 00 00 00 00 00 01 -> 4.9407E-324
=vStr(vPow(2;-1075)) -> 00 00 00 00 00 00 00 00 -> 0.00000E+00
You can see that XL somtimes manipulates the results of VBA as well as that the at maximum 15 digits of the XL/VBA output may not be sufficient do identify a DP number (in fact, it is well known that 17 digits are needed).
Doug, if different methods coded in two different languages and perhaps running on two different machines yield practically identical results then indeed there is strong evidence that they are precise as well as accurate. Be happy! Imho, when dealing with measured quantities obtained from a high-end 32bit-ADC, for example, all algorithms that are proven to be “harmless” will yield meaningfull results in DP. Contrary, some GPS applications use input quantities with 11 SF and thus require carefully designed algorithms. I’ve seen many standard libraries fail when being used to tackle this kind of problems.
Perhaps I’ve unintentionally created a wrong impression concerning Fortran. Even nowadays, there are many bad Fortran compilers. Imho, a modern compiler should handle subnormal numbers as well as exceptions. Among the free compilers I’ve tested so far, only the GNU gfortran compiler provides these functionalities but both the OpenWatcom and Salford FTN95 compilers do not. There is a nice test suite called ELEFUNT that can be used to test F (C and Java as well) compilers. Its latest version can be downloaded from here: http://www.math.utah.edu/~beebe/software/ieee/elefunt.tar.gz . Unfortunately, I have not found a transcription into VBA yet.
The current versions of xNumbers can be obtained from:
http://www.thetropicalevents.com/Xnumbers60.htm
According to my experience, John Beyers is a dedicated developer of xNumbers. Afaik, a new version will appear around Easter 2012. As VBA is out of the focus of my work at the moment, I’m only concerned with it during holidays.
I wish all of you a Happy New Year!
Georg
LikeLike
Sorry, using C&P in a hurry made me write meaningless sentences… I meant:
Imho, when dealing with measured quantities obtained from commonly used 20bit-ADCs (6 SF), for example, all algorithms that are proven to be “harmless” will yield meaningfull results in DP. Contrary, some GPS applications use input quantities with 11 SF (this even exceeds the accuracy of 9 SF of high-end 32bit-ADCs) and thus require carefully designed algorithms.
LikeLike
Extended precision data types don’t appear to be supported in any other MS language either, I don’t know why since FPU calculations are performed on the 80-bit wide stack anyway – at least on standard PCs. A VBA example where the extra bits make a difference is ?int(0.15*20) which returns 2, however ?int(cdbl(0.15*20)) gives the expected result of 3.
As Doug says, calculations are also affected by order of operations, for example using (35.0/30.0) with brackets either in column C of the spreadsheet or in the c code above gives consistent but slightly different results. Right-clicking the saved xlsm file, opening with winzip and dragging xml components into IE, allows one to see the precise cell values and calculation order.
For a thorough understanding of these matters one needs to look at the binary representations, as Georg already described, and also the disassembled code which is a lot less scary than one might think! This blog has some good insights, vba support for subnormal numbers is confirmed here.
Best wishes to all for 2012!
LikeLike
Hi Lori and Doug, to my great surprise, the example concerned with calculating 0.15*20 imho can be explained if it is assumed that VBA uses at least a guard digit in multiplication or even extended precision registers. If the result is turned into a DP “local” variable before further processing it, all is “fine”:
Function show_tests()
Dim a#
a = 0.15 * 20
MsgBox ("0.15*20 = " & CStr(0.15 * 20) & vbNewLine & _
"a = " & CStr(a) & vbNewLine & _
"int(0.15*20) = " & CStr(Int(0.15 * 20)) & vbNewLine & _
"int(a) = " & CStr(Int(a)) & vbNewLine & _
"int(cdbl((0.15*20)) = " & CStr(Int(CDbl(0.15 * 20)))) & vbNewLine & _
"int(cdbl((a)) = " & CStr(Int(CDbl(a)))
End Function
In extended precision, 0.15*20 is indeed smaller than 3 if it is calculated from two DP numbers:
1. 0.15 -> 3F C3 33 33 33 33 33 33.
2. split 20 into 16 + 4; both are powers of 2 => multiplication does not affect mantisssa bits
3. 0.15*4 = 0.6 -> 3F E3 33 33 33 33 33 33
4. 0.15*16 = 2.4 -> 40 03 33 33 33 33 33 33, which in binary format reads
O//IOOOOOO OOOO//OOII OOIIOOII OOIIOOII OOIIOOII OOIIOOII OOIIOOII OOIIOOII
5. multiply the exponent 3FE of 0.6 by 4 as well as shift the mantissa to the right by two places in order to obtain a number that can be added bitwise to 2.4. Take into account that the hidden 53rd bit now enters the mantissa of 0.6! The result is 0.6 =
O//IOOOOOO OOOO//OIOO IIOOIIOO IIOOIIOO IIOOIIOO IIOOIIOO IIOOIIOO IIOOIIOO.II
6. Add up and obtain 0.15*20 =
O//IOOOOOO OOOO//OIII IIIIIIII IIIIIIII IIIIIIII IIIIIIII IIIIIIII IIIIIIII.II
which is smaller than 3 -> 40 08 00 00 00 00 00 00 that in binary form reads:
O//IOOOOOO OOOO//IOOO OOOOOOOO OOOOOOOO OOOOOOOO OOOOOOOO OOOOOOOO OOOOOOOO
Of course, if the intermediate result is rounded to DP before it is used, 3 will be the result. But if the binary fraction is chopped by applying the int function, which in this case means to keep only the first displayed bit of the mantissa because it corresponds to 2^0, one obtains:
O//IOOOOOO OOOO//OOOO OOOOOOOO OOOOOOOO OOOOOOOO OOOOOOOO OOOOOOOO OOOOOOOO
which is 40 00 00 00 00 00 00 00 in hex notation and denotes a DP value of 2.
I’ve used “//” in order to indicate both the borders between sign bit and exponent bits as well as between exponent bits and mantissa bits.
phew….
LikeLike
Pingback: Comparing floating point numbers | Newton Excel Bach, not (just) an Excel Blog
Pingback: Floating Point Precision Problems | Newton Excel Bach, not (just) an Excel Blog