| 1 |
|
|---|
| 2 | /* @(#)e_atanh.c 1.3 95/01/18 */
|
|---|
| 3 | /*
|
|---|
| 4 | * ====================================================
|
|---|
| 5 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|---|
| 6 | *
|
|---|
| 7 | * Developed at SunSoft, a Sun Microsystems, Inc. business.
|
|---|
| 8 | * Permission to use, copy, modify, and distribute this
|
|---|
| 9 | * software is freely granted, provided that this notice
|
|---|
| 10 | * is preserved.
|
|---|
| 11 | * ====================================================
|
|---|
| 12 | *
|
|---|
| 13 | */
|
|---|
| 14 |
|
|---|
| 15 | #ifndef lint
|
|---|
| 16 | static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_atanh.c,v 1.7 2005/02/04 18:26:05 das Exp $";
|
|---|
| 17 | #endif
|
|---|
| 18 |
|
|---|
| 19 | /* __ieee754_atanh(x)
|
|---|
| 20 | * Method :
|
|---|
| 21 | * 1.Reduced x to positive by atanh(-x) = -atanh(x)
|
|---|
| 22 | * 2.For x>=0.5
|
|---|
| 23 | * 1 2x x
|
|---|
| 24 | * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
|
|---|
| 25 | * 2 1 - x 1 - x
|
|---|
| 26 | *
|
|---|
| 27 | * For x<0.5
|
|---|
| 28 | * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
|
|---|
| 29 | *
|
|---|
| 30 | * Special cases:
|
|---|
| 31 | * atanh(x) is NaN if |x| > 1 with signal;
|
|---|
| 32 | * atanh(NaN) is that NaN with no signal;
|
|---|
| 33 | * atanh(+-1) is +-INF with signal.
|
|---|
| 34 | *
|
|---|
| 35 | */
|
|---|
| 36 |
|
|---|
| 37 | #include "math.h"
|
|---|
| 38 | #include "math_private.h"
|
|---|
| 39 |
|
|---|
| 40 | static const double one = 1.0, huge = 1e300;
|
|---|
| 41 | static const double zero = 0.0;
|
|---|
| 42 |
|
|---|
| 43 | double
|
|---|
| 44 | __ieee754_atanh(double x)
|
|---|
| 45 | {
|
|---|
| 46 | double t;
|
|---|
| 47 | int32_t hx,ix;
|
|---|
| 48 | u_int32_t lx;
|
|---|
| 49 | EXTRACT_WORDS(hx,lx,x);
|
|---|
| 50 | ix = hx&0x7fffffff;
|
|---|
| 51 | if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
|
|---|
| 52 | return (x-x)/(x-x);
|
|---|
| 53 | if(ix==0x3ff00000)
|
|---|
| 54 | return x/zero;
|
|---|
| 55 | if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */
|
|---|
| 56 | SET_HIGH_WORD(x,ix);
|
|---|
| 57 | if(ix<0x3fe00000) { /* x < 0.5 */
|
|---|
| 58 | t = x+x;
|
|---|
| 59 | t = 0.5*log1p(t+t*x/(one-x));
|
|---|
| 60 | } else
|
|---|
| 61 | t = 0.5*log1p((x+x)/(one-x));
|
|---|
| 62 | if(hx>=0) return t; else return -t;
|
|---|
| 63 | }
|
|---|
