Subversion Repositories Kolibri OS

Rev

Rev 4872 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  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 */
  108.