Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 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 | |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 7 | // The original C code, the long comment, and the constants |
| 8 | // below are from FreeBSD's /usr/src/lib/msun/src/s_asinh.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 | // asinh(x) |
| 23 | // Method : |
| 24 | // Based on |
| 25 | // asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] |
| 26 | // we have |
| 27 | // asinh(x) := x if 1+x*x=1, |
| 28 | // := sign(x)*(log(x)+ln2)) for large |x|, else |
| 29 | // := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 30 | // := sign(x)*log1p(|x| + x**2/(1 + sqrt(1+x**2))) |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 31 | // |
| 32 | |
Charles L. Dorian | f273487 | 2012-04-06 14:01:12 -0400 | [diff] [blame] | 33 | // Asinh returns the inverse hyperbolic sine of x. |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 34 | // |
| 35 | // Special cases are: |
Charles L. Dorian | 94b0342 | 2011-12-08 17:07:13 -0500 | [diff] [blame] | 36 | // Asinh(±0) = ±0 |
Charles L. Dorian | c8d2544 | 2011-11-28 13:04:52 -0800 | [diff] [blame] | 37 | // Asinh(±Inf) = ±Inf |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 38 | // Asinh(NaN) = NaN |
| 39 | func Asinh(x float64) float64 { |
| 40 | const ( |
| 41 | Ln2 = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 42 | NearZero = 1.0 / (1 << 28) // 2**-28 |
| 43 | Large = 1 << 28 // 2**28 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 44 | ) |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 45 | // special cases |
Luuk van Dijk | 8dd3de4 | 2012-02-01 16:08:31 +0100 | [diff] [blame] | 46 | if IsNaN(x) || IsInf(x, 0) { |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 47 | return x |
| 48 | } |
| 49 | sign := false |
| 50 | if x < 0 { |
| 51 | x = -x |
| 52 | sign = true |
| 53 | } |
| 54 | var temp float64 |
| 55 | switch { |
| 56 | case x > Large: |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 57 | temp = Log(x) + Ln2 // |x| > 2**28 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 58 | case x > 2: |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 59 | temp = Log(2*x + 1/(Sqrt(x*x+1)+x)) // 2**28 > |x| > 2.0 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 60 | case x < NearZero: |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 61 | temp = x // |x| < 2**-28 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 62 | default: |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 63 | temp = Log1p(x + x*x/(1+Sqrt(1+x*x))) // 2.0 > |x| > 2**-28 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 64 | } |
| 65 | if sign { |
| 66 | temp = -temp |
| 67 | } |
| 68 | return temp |
| 69 | } |