question about magnitude of floating point variance for double

[LMHmedchem]

Hello,

I have some code comparing the angle between vectors. I am checking things like if the angle between a pair of vectors is > 90 degrees, or if two angles sum to > 180 degrees, etc.

Angle units are based on radians, which are based on Pi. To check if an angle is < 90 degrees,

if(largest_angle < 1.5707963268) { ...do something; }

or something like,

if(largest_angle < (atan(1)*4)*0.5) { ...do something; }

The issue is that pi goes to infinite decimal places and the angle I have calculated in double has a value that starts to become uncertain at some point to the right of the decimal. A reasonable comparison of the two values involves an effective method for handling differences that arise from floating point variation. It seems as if I should define a tolerance and add that to the logic resulting in three outcomes instead of two.

double tolerance = some_number;

// check if angle is > 90
if(largest_angle > ((atan(1)*4)*0.5)+tolerance ) { 
   // pretty confident the angle is > 90
   ...do something;
}
// check if angle is < 90
elseif(largest_angle < ((atan(1)*4)*0.5)-tolerance ) { 
   // pretty confident the angle is < 90
   ...do something;
}
else {
   // have to proceed on the assumption the the angle ~= 90
   ...do something;
}

This creates an uncertain range around (atan(1)*4)*0.5 +/- the tolerance. I'm not sure if this approach is the best way to handle the issue and if there is some other method that is considered more standard. I know that in managing floating point variance you are often fighting the balance between reproducibility and "correctness" in the calculation. If the method above is reasonable, I have no idea what to use for the tolerance.

It doesn't make sense to use degrees instead of radians even though 90.0000 degrees appears more certain. To get degrees, you have to convert radians using pi for a second time, degrees=(radians*180)/(pi), which just exasperates the problem. (apologizes for the elementary math, I know that many users here would know that already. Since some might not off the top of their head, I thought I would mention it anyway)

As always, suggestions would be greatly appreciated.

LMHmedchem

[VictorN]

Why do you think you need all these trigonometry calculations?
Maybe the elementary vector algebra would be enough?
just begin from the https://en.wikipedia.org/wiki/Dot_product

[LMHmedchem]

Why do you think you need all these trigonometry calculations?
Maybe the elementary vector algebra would be enough?
just begin from the https://en.wikipedia.org/wiki/Dot_product

I will have to think about that for a bit. The dot product involves the magnitude of the vectors, which in this case have no clear pattern. The dot product also gets bigger with the dimensional of the space, which is also not constant across multiple problems. I guess the question is whether or not I can learn everything I need to know from the dot product or not. It certainly would tell me if the the angle is over 90 or not, but I'm not sure about the other questions I need to have answered.

It's an interesting thought though, thanks.

LMHmedchem

[2kaud]

and the angle I have calculated in double has a value that starts to become uncertain at some point to the right of the decimal.

To how many DPs are you working/need to work?

Note that for VS in math.h there is the following

#if defined _USE_MATH_DEFINES && !defined _MATH_DEFINES_DEFINED
    #define _MATH_DEFINES_DEFINED
    // Definitions of useful mathematical constants
    //
    // Define _USE_MATH_DEFINES before including <math.h> to expose these macro
    // definitions for common math constants.  These are placed under an #ifdef
    // since these commonly-defined names are not part of the C or C++ standards
    #define M_E        2.71828182845904523536   // e
    #define M_LOG2E    1.44269504088896340736   // log2(e)
    #define M_LOG10E   0.434294481903251827651  // log10(e)
    #define M_LN2      0.693147180559945309417  // ln(2)
    #define M_LN10     2.30258509299404568402   // ln(10)
    #define M_PI       3.14159265358979323846   // pi
    #define M_PI_2     1.57079632679489661923   // pi/2
    #define M_PI_4     0.785398163397448309616  // pi/4
    #define M_1_PI     0.318309886183790671538  // 1/pi
    #define M_2_PI     0.636619772367581343076  // 2/pi
    #define M_2_SQRTPI 1.12837916709551257390   // 2/sqrt(pi)
    #define M_SQRT2    1.41421356237309504880   // sqrt(2)
    #define M_SQRT1_2  0.707106781186547524401  // 1/sqrt(2)
#endif

[LMHmedchem]

To how many DPs are you working/need to work?

This is part of my dilemma, I don't really know the answer to that. There are two basic cases I need to evaluate. The first is whether or not the angle between two vectors is < 90, ==90, or > 90. This is difficult when there is no explicit definition of the 90 degree cutoff and the computed angle will also have uncertainty at some point to the right of the decimal. Any incremental value > 90 changes the outcome of the answer. If the value is exactly 90, that is a different case. I am assuming that I have to live with some region where there is uncertainty.

I think that making the == 90 bin a bit wider like,

if(largest_angle > 90.00000001) {
   // angle is > than 90
}
if(largest_angle < 89.99999999) {
   // angle is < than 90
}
else {
   // angle == 90
}

makes the most sense for what I am doing but I'm not sure about how wide to make the == 90 bin. I would like the bin to be as wide as necessary to make sure that the > 90 and < 90 bins are certain. It is easier to live with uncertainty in the == 90 bin. I don't want to go wider than necessary, so I don't know what is a reasonable number to use.

I can play around with real data and see at what point I start getting answers that are different than existing methods for the same calculation but I thought I would post and see if there is some standard practice for this. I know the issue comes up in gaming allot.

As VictorN suggests, for the 90 degree test, I can just test the sign of the dot product. The dot product is exactly 0 at 90 degrees, so that would eliminate the floating point variation in one of the two values being compared. I don't know if I can do everything I need with the dot product because I need to evaluate some combinations of angles that sum to > 180 degrees.

Note that for VS in math.h there is the following

#if defined _USE_MATH_DEFINES && !defined _MATH_DEFINES_DEFINED
    #define _MATH_DEFINES_DEFINED
    // Definitions of useful mathematical constants
    //
    // Define _USE_MATH_DEFINES before including <math.h> to expose these macro
    // definitions for common math constants.  These are placed under an #ifdef
    // since these commonly-defined names are not part of the C or C++ standards
    #define M_E        2.71828182845904523536   // e
    #define M_LOG2E    1.44269504088896340736   // log2(e)
    #define M_LOG10E   0.434294481903251827651  // log10(e)
    #define M_LN2      0.693147180559945309417  // ln(2)
    #define M_LN10     2.30258509299404568402   // ln(10)
    #define M_PI       3.14159265358979323846   // pi
    #define M_PI_2     1.57079632679489661923   // pi/2
    #define M_PI_4     0.785398163397448309616  // pi/4
    #define M_1_PI     0.318309886183790671538  // 1/pi
    #define M_2_PI     0.636619772367581343076  // 2/pi
    #define M_2_SQRTPI 1.12837916709551257390   // 2/sqrt(pi)
    #define M_SQRT2    1.41421356237309504880   // sqrt(2)
    #define M_SQRT1_2  0.707106781186547524401  // 1/sqrt(2)
#endif

I don't use VS, but I think that math.h is defined in g++ as well. I don't know if it makes sense to use a predefined constant or to calculate the value on the fly.

LMHmedchem

[VictorN]

...
As VictorN suggests, for the 90 degree test, I can just test the sign of the dot product. The dot product is exactly 0 at 90 degrees, so that would eliminate the floating point variation in one of the two values being compared. I don't know if I can do everything I need with the dot product because I need to evaluate some combinations of angles that sum to > 180 degrees.

What do you mean by " some combinations of angles that sum to > 180 degrees?
And what is your primary data: vectors os angles or something else?

[2kaud]

I don't know if it makes sense to use a predefined constant or to calculate the value on the fly.

IMO I would use predefined constants over on-fly calculation unless greater precision was required to that offered by a predefined constant (I don't think that applies here?). If calculation is required I would consider having them as constexpr statements.

A common equality test is similar to

const double epsilon = ....;
double d1 = ...;
double d2 = ...;

if (abs(d1 - d2) < epsilon)
 // consider d1 and d2 equal.

where epsilon has to be determined depending upon required degree of DP accuracy.

If an equality test is done first and the test fails, then a simple > or < comparison should suffice to test for greater or lesser.

[LMHmedchem]

What do you mean by " some combinations of angles that sum to > 180 degrees?
And what is your primary data: vectors os angles or something else?

The primary data are rows of coordinates, so each of these can be considered as a vector from the origin to the point defined by the coordinates. I am working in 3 dimensions at present, but the algorithm is intended to work in high dimensional n-space.

The algorithm I am working on goes more or less like this,
- take a point in a data set to evaluate: point B (each point is processed in series starting with the closest to the centroid)
- generate a vector from point B to each other point in the data set (call each other point, point A(i), A(j), A(k), etc)
- sum all of the vectors BA(i,n) to create a sum vector BC (call the coordinates of the sum vector head point C)
- calculate the angle between each vector BA(i) and the sum vector BC

If all angles A(i)BC are < 90, point B is classified in category 1.

If any angles A(i)BC are > 90, additional processing is necessary.
- if any angle A(i)BC is > 90, evaluate the sum of the angles for each pair A(i)BC + A(j)BC to look for sums that are > 180.
- if A(i)BC + A(j)BC > 180, evaluate the angle A(i) B A(j)
- if any angle A(i) B A(j) is > 90, point B is classified in category 2
- if no angles A(i) B A(j) are > 90, point B is classified in category 1

The cases where an angle is exactly 90 or exactly 180 are special cases I haven't implemented yet. The algorithm works in 2D, but so far it doesn't work in 3D. I'm not sure why yet.

In 3D, what I am generating reminds me of the cross product, but since the vectors AB can be of all different magnitudes, I'm not sure that I can find a pattern in the cross product values that will answer my questions. My apologizes if my math conventions are not right in the description above, it has been quite a while since I did anything in this area.

I can post an excel spreadsheet with some data and plots if my explanation above is not clear, please let me know.

LMHmedchem

[VictorN]

...
If any angles A(i)BC are > 90, additional processing is necessary.
- [хл]if any angle A(i)BC is > 90, evaluate the sum of the angles for each pair A(i)BC + A(j)BC to look for sums that are > 180[/хл].
...

IMHO it makes sense only in case of coplanarity of these pairs. Doesn't it?