toolchain/gcc/newlib/libm/math/s_asinh.c

   1
   2 /* @(#)s_asinh.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, 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 FUNCTION
  16         <<asinh>>, <<asinhf>>---inverse hyperbolic sine
  17
  18 INDEX
  19         asinh
  20 INDEX
  21         asinhf
  22
  23 ANSI_SYNOPSIS
  24         #include <math.h>
  25         double asinh(double <[x]>);
  26         float asinhf(float <[x]>);
  27
  28 TRAD_SYNOPSIS
  29         #include <math.h>
  30         double asinh(<[x]>)
  31         double <[x]>;
  32
  33         float asinhf(<[x]>)
  34         float <[x]>;
  35
  36 DESCRIPTION
  37 <<asinh>> calculates the inverse hyperbolic sine of <[x]>.
  38 <<asinh>> is defined as
  39 @ifnottex
  40 . sgn(<[x]>) * log(abs(<[x]>) + sqrt(1+<[x]>*<[x]>))
  41 @end ifnottex
  42 @tex
  43 $$sign(x) \times ln\Bigl(|x| + \sqrt{1+x^2}\Bigr)$$
  44 @end tex
  45
  46 <<asinhf>> is identical, other than taking and returning floats.
  47
  48 RETURNS
  49 <<asinh>> and <<asinhf>> return the calculated value.
  50
  51 PORTABILITY
  52 Neither <<asinh>> nor <<asinhf>> are ANSI C.
  53
  54 */
  55
  56 /* asinh(x)
  57  * Method :
  58  *      Based on
  59  *              asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
  60  *      we have
  61  *      asinh(x) := x  if  1+x*x=1,
  62  *               := sign(x)*(log(x)+ln2)) for large |x|, else
  63  *               := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
  64  *               := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
  65  */
  66
  67 #include "fdlibm.h"
  68
  69 #ifndef _DOUBLE_IS_32BITS
  70
  71 #ifdef __STDC__
  72 static const double
  73 #else
  74 static double
  75 #endif
  76 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
  77 ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
  78 huge=  1.00000000000000000000e+300;
  79
  80 #ifdef __STDC__
  81         double asinh(double x)
  82 #else
  83         double asinh(x)
  84         double x;
  85 #endif
  86 {
  87         double t,w;
  88         __int32_t hx,ix;
  89         GET_HIGH_WORD(hx,x);
  90         ix = hx&0x7fffffff;
  91         if(ix>=0x7ff00000) return x+x;  /* x is inf or NaN */
  92         if(ix< 0x3e300000) {    /* |x|<2**-28 */
  93             if(huge+x>one) return x;    /* return x inexact except 0 */
  94         }
  95         if(ix>0x41b00000) {     /* |x| > 2**28 */
  96             w = __ieee754_log(fabs(x))+ln2;
  97         } else if (ix>0x40000000) {     /* 2**28 > |x| > 2.0 */
  98             t = fabs(x);
  99             w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
 100         } else {                /* 2.0 > |x| > 2**-28 */
 101             t = x*x;
 102             w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
 103         }
 104         if(hx>0) return w; else return -w;
 105 }
 106
 107 #endif /* _DOUBLE_IS_32BITS */