In an earlier post I discussed the issue of floating point accuracy in calculations. I have another example to share today that has been lurking in the Mapping Toolbox product. The issue is with the function wrapTo180.
Contents
Work-around
There is a work-around for those of you using Mapping Toolbox with releases R2007b and R2008a.
Numerical Issue
In R2007b and R2008a, the implementation for wrapTo180 causes results from wrapTo180(lon) to differ very slightly from lon even for for certain values of lon in the interval [-180 180]. The differences are on the order of 2*eps(lon).
Test Case
The following should evaluate to true.
lon = 115.8323;
but doesn't.
isequal(lon,wrapTo180(lon))
ans =
0
These quantities aren't equal because of the arithmetic being done: 115.8323 + 180 - 180 differs from 115.8323 by
(lon + 180 - 180) - lon
ans = 2.8422e-014
which happens to be
2*eps(lon)
ans = 2.8422e-014
The Code
Here's the code for wrapTo180
type wrapTo180function lon = wrapTo180(lon) %wrapTo180 Wrap angle in degrees to [-180 180] % % lonWrapped = wrapTo180(LON) wraps angles in LON, in degrees, to the % interval [-180 180] such that 180 maps to 180 and -180 maps to -180. % (In general, odd, positive multiples of 180 map to 180 and odd, % negative multiples of 180 map to -180.) % % See also wrapTo360, wrapTo2Pi, wrapToPi. % Copyright 2007 The MathWorks, Inc. % $Revision: 1.1.6.2 $ $Date: 2007/08/20 16:35:59 $ lon = wrapTo360(lon + 180) - 180;
Simple Fix
We can just skip the arithmetic operation for values,like 115.8323, that already are in the interval [-180 180]. A nice way to do this is with logical indexing and we can replace the line of code with this.
q = (lon < -180) | (180 < lon); lon(q) = wrapTo360(lon(q) + 180) - 180;
Comments?
Have you ever had a similar issue, where a seemingly innocent code expression causes erroneous results? How have you found this out and solved it? I'd love to hear your thoughts here.
Get
the MATLAB code
Published with MATLAB® 7.6
Interesting post. In grad school I would often search short lists I had made, but couldn’t find certain values. For example, S=[0.24:0.02:0.46]; find(S==0.34) returns an empty set. S(6)-0.34 yields -5.5e-17. Since then I’ve started searching using find(abs(S-0.34)<eps) which yields the correct answer.
Exactly - checking for equality is not appropriate for floating point calculations.
–Loren
Checking for equality is OK in certain circumstances: With “flints”, for example, (p = -10:10 for example), or for values that you know can be computed without error (such as 1./4.). But you have to know everything you’ve done to a number … which includes built-in functions for which you can’t see the code. So it’s not always possible.
The other place where checking for equality makes sense is when an entry is exactly zero. That happens all the time in sparse matrix computations … those entries get purged by MATLAB. Roundoff error can differ on on different machines, so nnz(A) can vary from machine to machine even though the code and the input data are the same.
Another example of checking for floating point equality with zero is gallery(’house’,…). If norm(x) is already zero, H*x does not need to do nothing so H=I is chosen. Try:
edit private/house.m
It checks for exact zero, not eps, because H is computed safely if norm(x) is O(eps).
Dave: the check for abs(S-c)<eps only works if the value c you’re checking against is abs(c)<.5 or so. If S is larger then that won’t work. You should use something like abs(S-c)<4*eps(c), for example, where “4″ is something small. Then if you compare S with c=34 instead of c=0.34, it will still work. Otherwise you might as well be testing “if(S==c)”. Type
help eps
for more details.
A collegue encountered this recently:
>> clear all
>> qsingle = single(0.4)
qsingle =
0.4
>> qdouble = double(qsingle)
qdouble =
0.400000005960464
>> whos
Name Size Bytes Class Attributes
qdouble 1×1 8 double
qsingle 1×1 4 single
I presume the “extra” digits are truncated for the single variable (as fewer bits for the mantissa).
Ljubomir
The extra digits are the representation error of the single precision number as compared to the double. Suppose we represent 199.2 in 3-digit decimal floating point (and I’ll normalize my mantissas however I want). We get 199×10^0. The loss of 0.2 here is representation error. Something like this also happens when you represent 4/10 as a double–double(4/10) is not exactly 4/10. Back to the example, now convert this 3-digit result to 2-digit decimal floating point. This time we round up: 20×10^1. Converting that back to 3-digit we get 200×10^0, not the 199×10^0 we had before. So here we picked up +1 by converting from 3-digits to 2-digits and then back to 3-digits.
Hm, in your example both 2-digit and 3-digit numbers are 200. In mine the single and the double differ. Analogy would be if you got something else then 200×10^0 when converting 20×10^1 to 3 digits. Don’t understand how the single -> double manissa extension (the extra bits were set to 0s?) yielded the difference? Maybe is just a printing issue in the command window (singles have less significant digits printed)?
Ljubomir-
It’s not just printing. Converting to a representation that can’t exactly hold the number causes some loss of real information. That’s what happens when you go from the mathematical number 0.4 to holding it in double precision, with 53 bits for the mantissa. Single has even fewer bits to hold the information so even more information gets “lost”. Then, converting back to double, the computer takes the now appproximate single value and doesn’t get the double value exactly as it was (which may not have been exact in the first place).
With 3 digits, we can’t accurately represent 4 digits decimal. It’s the same issue. Not printing - no place in memory!
–Loren
Thanks Loren. I understand that maybe “…the “extra” digits are truncated for the single variable (as fewer bits for the mantissa)” whereas they may appear in the double. That’s fine. What puzzles me is that the *conversion* of single->double yielding this i.e. it happens on extending the single (already truncated) into double. I take I wrongly assume promotion single->double involves only extending the extra mantissa bits with zeros (which I take can not add extra digits)?
Ljubomir
Ljubomir, using format hex will show you what’s really going on:
>> format hex
>> .4
ans =
3fd999999999999a
>> qsingle = single(0.4)
qsingle =
3ecccccd
>> qdouble = double(qsingle)
qdouble =
3fd99999a0000000
.4 as a double is represented by a 0 sign bit, 01111111101 (1023-biased) exponent, and [1]10011001…10011010 mantissa, where the [1] is implied. Cast that to single and you get 0 sign bit, 01111101 (127-biased) exponent, and [1]10011001100110011001101 mantissa, which is the 52 bit d.p. mantissa, rounded to 23 bits. Cast that back to double and you get the exact same mantissa, but zero extended.
Everything else is just an artifact of printing a binary number in base 10.
Sorry for the delay, Ljubomir. I’ve been away. My decimal example was formulated to show the essential principle without the distractions of binary floating point. At any rate, I agree with Peter’s exposition. You may round up or down when losing mantissa bits in the conversion from double to single. When converting back to double, the mantissa, such as it is, is padded with zeros. The phenomenon is obfuscated, however, by the conversion of the binary form to decimal form for display.
Thanks guys, most helpful, all clear now wrt my example, format hex separates presentation from content nicely. :-)
Ljubomir