1  /*
       2   * ====================================================
       3   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       4   *
       5   * Developed at SunPro, a Sun Microsystems, Inc. business.
       6   * Permission to use, copy, modify, and distribute this
       7   * software is freely granted, provided that this notice
       8   * is preserved.
       9   * ====================================================
      10   */
      11  
      12  /* __ieee754_acosh(x)
      13   * Method :
      14   *	Based on
      15   *		acosh(x) = log [ x + sqrt(x*x-1) ]
      16   *	we have
      17   *		acosh(x) := log(x)+ln2,	if x is large; else
      18   *		acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
      19   *		acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
      20   *
      21   * Special cases:
      22   *	acosh(x) is NaN with signal if x<1.
      23   *	acosh(NaN) is NaN without signal.
      24   */
      25  
      26  #include <math.h>
      27  #include <math_private.h>
      28  #include <libm-alias-finite.h>
      29  
      30  static const double
      31  one	= 1.0,
      32  ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
      33  
      34  double
      35  __ieee754_acosh (double x)
      36  {
      37    int64_t hx;
      38    EXTRACT_WORDS64 (hx, x);
      39  
      40    if (hx > INT64_C (0x4000000000000000))
      41      {
      42        if (__glibc_unlikely (hx >= INT64_C (0x41b0000000000000)))
      43  	{
      44  	  /* x > 2**28 */
      45  	  if (hx >= INT64_C (0x7ff0000000000000))
      46  	    /* x is inf of NaN */
      47  	    return x + x;
      48  	  else
      49  	    return __ieee754_log (x) + ln2;/* acosh(huge)=log(2x) */
      50  	}
      51  
      52        /* 2**28 > x > 2 */
      53        double t = x * x;
      54        return __ieee754_log (2.0 * x - one / (x + sqrt (t - one)));
      55      }
      56    else if (__glibc_likely (hx > INT64_C (0x3ff0000000000000)))
      57      {
      58        /* 1<x<2 */
      59        double t = x - one;
      60        return __log1p (t + sqrt (2.0 * t + t * t));
      61      }
      62    else if (__glibc_likely (hx == INT64_C (0x3ff0000000000000)))
      63      return 0.0;				/* acosh(1) = 0 */
      64    else					/* x < 1 */
      65      return (x - x) / (x - x);
      66  }
      67  libm_alias_finite (__ieee754_acosh, __acosh)