

Title  Bug: Watch for floating point errors 
Description  This bug alert points out some possible sources of error in floating point calculations. 
Keywords  floating point errors, single, double, overflow, underflow, Visual Basic 6, VB 6, VB.NET, VB 2005 
Categories  Algorithms, Bug Alerts 


Ronny Zulaikha pointed out a potential problem with the example Find the distance between a point and a line segment. The code performs the following check:
If dx = 0 And dy = 0 Then
If this test is false (so dx or dy is not 0), then it performs this calculation:
t = ((px  X1) * dx + (py  Y1) * dy) _
/ (dx * dx + dy * dy)
As Ronny pointed out, dx and dy might not be zero but dx * dx + dy * dy might be zero. When you multiply two small numbers together, the result is even smaller. If dx is small enough, then dx * dx might be so small that a Single might not be able to represent it other than as zero. Then the calculation of t will result in a divide by zero error.
In fact, the problem is even worse than this. If dx * dx + dy * dy is small enough and the numerator is large enough, the resulting calculation might lead to a floating point overflow even if the denominator is not zero.
In the example program, if you enter a value for X and click the Calculate button, the program calculates X * X and 1 / (X * X). If X is 1E19, then both calculations produce valid results. If X is 1E20, then X * X is a very small nonzero number but 1 / (X * X) produces an overflow error. If X = 1E23, then X * X equals 0.
The following code shows how the program performs its calculations.


Private Sub cmdCalculate_Click()
Dim X As Single
Dim result1 As Single
Dim result2 As Single
X = CSng(txtX.Text)
On Error Resume Next
result1 = X * X
If Err.Number = 0 Then
lblResult1.Caption = Format$(result1)
Else
lblResult1.Caption = Err.Description
End If
result2 = 1 / (X * X)
If Err.Number = 0 Then
lblResult2.Caption = Format$(result2)
Else
lblResult2.Caption = Err.Description
End If
End Sub


Unfortunately there isn't a simple way to protect a program against these sorts of errors. In some applications, you may be able to perform some sanity checks. For example, if you know that a fraction's numerator and denominator must be between 10 and 100, then you should be in good shape. Otherwise you should perform whatever checks you can and then use error handling to protect against unexpected errors.
Bruce Deam wrote the following expanded discussion:
Loss of Significance with Numerical Subtraction
Subtracting two numbers that are similar in magnitude introduces significant error into the result and should therefore be avoided wherever possible. Unfortunately subtraction is a common operation in geometric calculations (e.g. computer graphics) and finite difference arithmetic (interpolation).
Consider the apparently simple calculation of the distance between a point and a line segment. The following algorithm by Stephens (2007) calculates the distance using a parameter, t, which is used to ascertain whether one end or the other is nearest the point or whether the nearest point is somewhere within the segment.


' Calculate the distance between the point and the segment.
Private Function DistToSegment( _
ByVal px As Single, ByVal py As Single, _
ByVal X1 As Single, ByVal Y1 As Single, _
ByVal X2 As Single, ByVal Y2 As Single, _
ByRef near_x As Single, _
ByRef near_y As Single) As Single
Dim dx As Single
Dim dy As Single
Dim t As Single
dx = X2  X1
dy = Y2  Y1
If dx = 0 And dy = 0 Then
' It's a point not a line segment.
dx = px  X1
dy = py  Y1
near_x = X1
near_y = Y1
DistToSegment = Sqr(dx * dx + dy * dy)
Exit Function
End If
' Calculate the t that minimizes the distance.
t = ((px  X1) * dx + (py  Y1) * dy) / _
(dx * dx + dy * dy)
' See if this represents one of the segment's
' end points or a point in the middle.
If t < 0 Then
dx = px  X1
dy = py  Y1
near_x = X1
near_y = Y1
ElseIf t > 1 Then
dx = px  X2
dy = py  Y2
near_x = X2
near_y = Y2
Else
near_x = X1 + t * dx
near_y = Y1 + t * dy
dx = px  near_x
dy = py  near_y
End If
DistToSegment = Sqr(dx * dx + dy * dy)
End Function


In this implementation, the calculation of the parameter t can cause trouble when the x and y components of the line segment's length (ie dx and dy respectively) are both close to zero because the denominator of the equation also becomes close to zero. This is accentuated by the values of the two components being squared. In practice, this is likely to be a common situation because many computer graphics calculations only use single precision arithmetic.
A simple improvement to the calculation of parameter t would be to divide both the top and bottom lines of the calculation by (say) dx, and introducing a second parameter, s = dy/dx. This of course would require comparison between dy and dx and use of a similar expression in terms of dx/dy to avoid the same problem reoccurring. E.g.


' Calculate the t that minimizes the distance.
If Abs(dx) > Abs(dy) Then
s = dy / dx
t = ((px  X1) + (py  Y1) * s) / dx / (1 + s * s)
Else
s = dx / dy
t = ((px  X1) * s + (py  Y1)) / dy / (1 + s * s)
End If


This local improvement, however, fails to recognise the whole problem (as well as accommodating 45 degree line segments where dx = dy). Once a problem is identified, it is usually better to consider both how the values used in the calculation are themselves calculated and how the result of the calculation is subsequently used. The segment length components are calculated by subtracting two (possibly similar magnitude) numbers. We will investigate whether this can be avoided later.
The first comparison following the calculation of t (ie whether t is smaller than zero) could make use of the fact that the denominator is always positive and just use the numerator to make this comparison, particularly as t is not required for any subsequent calculation if the comparison is true. E.g.:
If (px  X1) * dx + (py  Y1) * dy < 0 Then
The test for it being closest to the other end of the line segment can compare the numerator of t with its denominator to avoid dividing by the denominator. E.g.
ElseIf (px  X2) * dx + (py  Y2) * dy >= dx * dx + dy * dy Then
Use of these two comparisons allows the value of t to be rearranged to be calculated more accurately using a modified form of the simple improvement suggested above. Note that the following expressions have only one subtraction and only include dx and dy as a ratio, which will always be smaller than 1:
ElseIf Abs(dx) > Abs(dy)
s = dy / dx
near_x = (px + s * (py  Y1 + s * X1)) / (1 + s * s)
near_y = (Y1 + s * (px  X1 + s * py)) / (1 + s * s)
Finally, the original algorithm used a simple comparison at the beginning that would provide the best efficiency when the function is often called with two points rather than a point and a line segment. If this was an infrequent situation, the special case can be included in the main calculation, giving the following implementation:


' Calculate the distance between the point and the segment.
Private Function DistToSegment( _
ByVal px As Single, ByVal py As Single, _
ByVal X1 As Single, ByVal Y1 As Single, _
ByVal X2 As Single, ByVal Y2 As Single, _
ByRef near_x As Single, ByRef near_y As Single) _
As Single
Dim dx As Single
Dim dy As Single
Dim test As Single
Dim s As Single
dx = X2  X1
dy = Y2  Y1
' See if this represents one of the segment's
' end points or a point in the middle.
test = (px  X1) * dx + (py  Y1) * dy
If test <= 0 Then
near_x = X1
near_y = Y1
ElseIf test >= dx * dx + dy * dy Then
near_x = X2
near_y = Y2
ElseIf Abs(dx) > Abs(dy)
s = dy / dx
near_x = (px + s * (py  Y1 + s * X1)) / (1 + s * s)
near_y = (Y1 + s * (px  X1 + s * py)) / (1 + s * s)
ElseIf Abs(dx) < Abs(dy)
s = dx / dy
near_x = (X1 + s * (py  Y1 + s * px)) / (1 + s * s)
near_y = (py + s * (px  X1 + s * Y1)) / (1 + s * s)
Else
s = Sgn(dx * dy)
near_x = (X1 + s * (py  Y1 + s * px)) / 2
near_y = (py + s * (px  X1 + s * Y1)) / 2
End If
dx = px  near_x
dy = py  near_y
DistToSegment = Sqr(dx * dx + dy * dy)
End Function


A further improvement in efficiency may be obtained by precalculating values that are used in several other calculations; however this makes the code even harder to read and may not provide much improvement on a modern processor if it needed to be stored outside the fast internal cache memory.
A final aspect of accommodating the potential loss of precision introduced by subtracting numbers is that the simple concept embodied in the original algorithm has now become much more obscure. It also requires much more elaborate documentation to describe why it was implemented in that manner.
Kahn (1997) provides good notes on a similar problem for triangles. There are a number of similar articles available from his home page http://www.cs.berkeley.edu/~wkahan/
[Note that the original example program lets the user pick coordinates from the form in native .NET coordinates, which are pixels. Because the user cannot click on fractional pixel values, the points on the line segment are either identical (so dx = dy = 0), or either dx or dy is at least 1. That means this should not be an issue in this particular application.  Rod]





