A sample library function

  #include <Numerics.h>
  #include <FPEnvironment.h>
  
  
#pragma fenv_access on
  double SampleArcTanh(double x) {
      TFPEnvironmentSaveAndUpdate environment;
  

      double_t tmp = Abs(x);
      if (IsLess(tmp, 0x1.0p-27)) { // threshold 2^-27
          if (tmp > 0.0)
              environment.SetFlags(kFPInexact);
      } else
          tmp = 0.5 * Log1Plus(2.0 * tmp / (1.0 - tmp));
      return CopySign(tmp, x);
  }

The temporary variable tmp has the type double_t, the platform's most efficient type at least as wide as double.

An important feature of robust numerical functions is their management of the execution environment--the rounding method employed and the reporting of exceptions. The computation of SampleArcTanh is bracketed by the constructor and destructor of the TFPEnvironmentSaveAndUpdate environment object. After the incoming environment is saved by the constructor, the call to Reset forces rounding to nearest and clears all exceptions. (By convention, functions are called with rounding to nearest unless otherwise specified. SampleArcTanh works only with rounding to nearest.) When the function is complete, the destructor restores the incoming state and raises the flags accumulated by SampleArcTanh. In this way, SampleArcTanh executes in the environment of its choosing and signals precisely the exceptions that pertain to its result.

Computing ArcTanh

The mathematical function arctanh is an odd function (that is, ), so for convenience the computation proceeds with tmp, the absolute value of the argument. The argument's sign is copied back onto the result later. The approximation splits into two cases.

1. If tmp is so small that its square would round off completely when added to 1, then mathematically to within a rounding error.
   The threshold for this test, 0x1.0p27, exhibits the hexadecimal floating constant syntax. Because the double format has 53 significant bits, is a reasonable choice for the threshold (even if double_t is wider than double). The fragment 0x1.0 of the constant is a fixed-point string of hexadecimal digits with a radix point separating its integer and fraction parts. Its value is 1. The fragment p27 signifies the scale factor .
   The function IsLess performs the test to filter small operands. Unlike its built-in counterpart <, IsLess does not raise the invalid operation flag when either operand is a NaN. Once the argument is known to be small, it remains to distinguish the special case 0, because exactly. At this stage the argument is known not to be a NaN (which compares unordered with, rather than less than, equal to, or greater than, any quantity), so the ordinary C++ relational operator suffices. The call to the member function SetFlags records the rounding error.
   

1. For magnitudes above the threshold but no greater than 1, the formula
   
   applies. A more numerically stable form of the expression is
   
   which appears in the code. (See "Rewriting mathematical expressions" on page 225 for a comparison of the two formulas.) Note the use of the function Log1Plus, which computes the natural logarithm of its argument as if added to 1 but without performing the addition--in this case without rounding off significant bits of
   .
   This branch of the code incorporates several special cases as well:
   
   • If the magnitude of the argument is 1, then the formula produces the floating-point value and raises the divide-by-zero flag (in the spirit of the underlying arithmetic).
   • If the magnitude is greater than 1, the formula produces a NaN and raises the invalid flag.
   • If the argument is a NaN, that NaN simply passes through the computation without exception.

Although this computation boils down to a few tests and one brief formula, it exploits the underlying arithmetic to handle every possible input. A more detailed analysis than perhaps warranted by this first example would show that tiny values, though inexact and subnormal, are correctly rounded to full double precision, so that SampleArcTanh need never raise underflow. Also, though the function arctanh approaches as its argument approaches 1, it approaches so steeply that SampleArcTanh never exceeds the value 19 for double numbers.