ACOS ARGUMENT CLOSE TO UNITY

[Anthony_Richards]

I came across the following awkward behaviour:

an argument being used in ACOS function displayed in the debugger shows a value 1.00000000000000 but it resulted an a Nan result.

When displayed as hexadecimal,  the argument shows a value #3FF0000000000001.

A PARAMETER initialised to 1.0d+00 displays as  1.00000000000000 and when displayed as hexadecimal shows #3FF0000000000000.

This means that every time before I use ACOS, I have to test the magnitude of the argument against an exact unity value and make a decision based on the size of the difference from unity. This could become an annoying overhead - which surely will also affect the ASIN function as well. Would it be more efficient to use ISNAN on the result and set the result to zero or PI, depending on the sign of the argument, than test the argument beforehand?

[TimP]

I suppose max(min(arg,1_real64),-1_real64) might be more efficient in some contexts, if simply clamping the range is sufficient in your usage.  You may not be able to see the difference if it's not a question of maintaining vectorization. 

[mecej4]

What would you have done had the argument to ACOS been, say, 1.57? Perhaps the underlying issue has to do with some calculation which, if done in a stable and accurate way, should never have yielded a value with absolute value larger than 1.0.

[jimdempseyatthecove]

IF(X .GT. 1_REAL64 .OR. X .LT. -1_REAL64) CALL YourAssertRoutine("Argument out of range")

It is the responsibility of the caller to assert validity of arguments.

Jim Dempsey

[FortranFan]

Anthony Richards wrote:

.. This means that every time before I use ACOS, I have to test the magnitude of the argument against an exact unity value and make a decision based on the size of the difference from unity. This could become an annoying overhead - which surely will also affect the ASIN function as well. Would it be more efficient to use ISNAN on the result and set the result to zero or PI, depending on the sign of the argument, than test the argument beforehand? ..

Given the IEEE facility in standard Fortran since 2003 revision, my take is that it is best to do any check post-calculation.  I expect IEEE_INEXACT flag will be turned on if any of the 4 inverse trigonometry functions are invoked outside the allowed domain of x.  Making use of this flag using IEEE standard intrinsic module facilities in perhaps a custom inverse function that normally simply invokes the straight intrinsic inverse trigonometry function but which does appropriate corrective post-processing action for invalid x, as suitable for one's needs, is worth a consideration.

[Simon_Geard]

We had this problem in our code where we need the angle between two vectors and the direction of a vector. In the latter case we found that even if the vector had been normalized it as still possible for the argument to acos to be slightly greater than 1. To fix the problem we switched to using atan2. You don't have to worry about edge cases and it is certainly faster.