# How to properly decide what is the appropriate rounding off error for a mathematical function?

In the space of 3D geometry, I need to compute the magnitude of the cross product of two `Line3D` because I want to check whether these two Lines are "parallel" to each other.

`Line3D` is consisted of 4 fields, `x1`, `y1`, `x2`, `y2` denoting the two points. the `x` and `y` are of the type `double`

In theory, two lines are parallel if the cross product of their unit vectors are equal to zero.

But in practice, C# double type has rounding errors in the arithmetic operation. So I need to prescribe a tolerance for it. How should I proceed to decide what is correct magnitude of the tolerance, for this function specifically, or or any other types of arithmetic generally?

Any rule of thumbs or even better, mathematically way to decide the tolerance?

• given that the x and y coordinate will have the same rounding problems surely you can never have parallel lines except in trivial cases. You need to decide how parallel is parallel for your purposes – Ewan Sep 21 '18 at 9:38

## 4 Answers

Your space does not go to infinity, so you'll be able to calculate an angle between to lines that will result in a difference in your space that you'll still call parallel.

If you'll display the lines on a screen you'll might go for a difference that will result in two parallel lines only differ by a few pixels.

Your space could be 25,000 by 25,000 by 25,000 so the longest distance is the distance from one corner of the cube to the opposite corner.

What you then do is to create two lines starting from the one corner and end at the other corner with one line missing by say 10. Then you calculate the cross product of those lines and the result is your target for lines being parallel.

(a) The correct amount of rounding depends on your problem domain. E.g. you might decide that lines within 1° of each other should be considered parallel.

(b) The floating-point representation puts a lower bound on the necessary rounding. Determining error propagation is a large subject in the field of numerical analysis.

All finite floating point numbers have a predecessor and successor number. This “step size” between numbers or “machine epsilon” depends on the magnitude of the number, i.e. is smallest around zero and grows towards positive and negative infinity. But conveniently, this is a constant relative error. For a 64-bit double there are 53 significant bits, leading to a relative step size of around 2^-52 = 2.22E-16 (the machine epsilon is usually defined as half that).

Assuming that your input values are as exact as possible (i.e. their relative error is the machine epsilon) you can trace your calculations to see how these errors propagate. E.g. for an addition or subtraction, you must add the absolute errors. Note that e.g. subtracting two similar values might result in a value that is on the order of its absolute error (catastrophic cancellation). It is sometimes possible to re-arrange calculations to avoid this. For a multiplication or division, you can add the relative errors.

You can use the results of your analysis to calculate suitable absolute and relative bounds at runtime.

As a super simple example, consider a test for equality of two floating point numbers a and b, under the assumption that their error is the machine epsilon. One simple approach would be to compare the next and previous float, i.e. `prev(a) <= b && b <= next(a)`. If the two numbers have arbitrary errors this is more difficult but the idea is the same: derive absolute error bounds and compare if the bounded area overlaps. This can be done conveniently by comparing their difference, i.e. `abs(a - b) <= absolute_error`.

For a thorough treatment of floating point comparisons, read https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/. That page also discusses when it's more appropriate to count ULPs (number of discrete floats between two numbers) than to deal with machine epsilons.

Choosing the right tolerance can only be done with knowledge of the problem domain.

One needs to know what the "lines" are used for, what they represent, when they are "parallel enough" for the particular use case, and what the parallelity test is used for. Then one can try to define a metrics for this use case and use it for determining a tolerance value for the cross product.

For example, if the use case involves just graphical lines on a screen and "parallel enough" is defined by "when the lines look parallel", then you will need quite a different tolerance than, for example, for lines representing boundaries of some object on a map, measured in meters, kilometers or millimeters, or if the lines represent some abstract mathematical object.

If you are going to write a general purpose library with such a test, which shall be reusable for many different use cases, the only sensible way is to provide a tolerance parameter for any function where it is required. Note that it may turn out that the distance of the cross product from zero is not necessarily the best metrics for the specific case.

The best way it to create your own function

``````double CrossProductMagnitude( Line3D l1, Line3d l2 )
``````

and test it with a variety of inputs, until you decide for yourself how small a result is appropriate for the lines to be considered parallel.

For example

if l1 is (0,0) to (1000,1000) and l2 is (1,0) to (x,1000), are the lines considered parallel if x = 1002 / 1001.5 / 1001.01 / 1001.000001?

Only you can decide that.