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