acosh.mx raw

   1  // Copyright 2010 The Go Authors. All rights reserved.
   2  // Use of this source code is governed by a BSD-style
   3  // license that can be found in the LICENSE file.
   4  
   5  package math
   6  
   7  // The original C code, the long comment, and the constants
   8  // below are from FreeBSD's /usr/src/lib/msun/src/e_acosh.c
   9  // and came with this notice. The go code is a simplified
  10  // version of the original C.
  11  //
  12  // ====================================================
  13  // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  14  //
  15  // Developed at SunPro, a Sun Microsystems, Inc. business.
  16  // Permission to use, copy, modify, and distribute this
  17  // software is freely granted, provided that this notice
  18  // is preserved.
  19  // ====================================================
  20  //
  21  //
  22  // __ieee754_acosh(x)
  23  // Method :
  24  //	Based on
  25  //	        acosh(x) = log [ x + sqrt(x*x-1) ]
  26  //	we have
  27  //	        acosh(x) := log(x)+ln2,	if x is large; else
  28  //	        acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
  29  //	        acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
  30  //
  31  // Special cases:
  32  //	acosh(x) is NaN with signal if x<1.
  33  //	acosh(NaN) is NaN without signal.
  34  //
  35  
  36  // Acosh returns the inverse hyperbolic cosine of x.
  37  //
  38  // Special cases are:
  39  //
  40  //	Acosh(+Inf) = +Inf
  41  //	Acosh(x) = NaN if x < 1
  42  //	Acosh(NaN) = NaN
  43  func Acosh(x float64) float64 {
  44  	if haveArchAcosh {
  45  		return archAcosh(x)
  46  	}
  47  	return acosh(x)
  48  }
  49  
  50  func acosh(x float64) float64 {
  51  	const Large = 1 << 28 // 2**28
  52  	// first case is special case
  53  	switch {
  54  	case x < 1 || IsNaN(x):
  55  		return NaN()
  56  	case x == 1:
  57  		return 0
  58  	case x >= Large:
  59  		return Log(x) + Ln2 // x > 2**28
  60  	case x > 2:
  61  		return Log(2*x - 1/(x+Sqrt(x*x-1))) // 2**28 > x > 2
  62  	}
  63  	t := x - 1
  64  	return Log1p(t + Sqrt(2*t+t*t)) // 2 >= x > 1
  65  }
  66