newlib/libm/mathfp/e_acosh.c

   1
   2 /* @(#)e_acosh.c 5.1 93/09/24 */
   3
   4 /*
   5 FUNCTION
   6 <<acosh>>, <<acoshf>>---inverse hyperbolic cosine
   7
   8 INDEX
   9 acosh
  10 INDEX
  11 acoshf
  12
  13 SYNOPSIS
  14         #include <math.h>
  15         double acosh(double <[x]>);
  16         float acoshf(float <[x]>);
  17
  18 DESCRIPTION
  19 <<acosh>> calculates the inverse hyperbolic cosine of <[x]>.
  20 <<acosh>> is defined as
  21 @ifnottex
  22 . log(<[x]> + sqrt(<[x]>*<[x]>-1))
  23 @end ifnottex
  24 @tex
  25 $$ln\Bigl(x + \sqrt{x^2-1}\Bigr)$$
  26 @end tex
  27
  28 <[x]> must be a number greater than or equal to 1.
  29
  30 <<acoshf>> is identical, other than taking and returning floats.
  31
  32 RETURNS
  33 <<acosh>> and <<acoshf>> return the calculated value.  If <[x]>
  34 less than 1, the return value is NaN and <<errno>> is set to <<EDOM>>.
  35
  36 PORTABILITY
  37 Neither <<acosh>> nor <<acoshf>> are ANSI C.  They are not recommended
  38 for portable programs.
  39
  40
  41 QUICKREF
  42  ansi posix rentrant
  43  acos    n,n,m
  44  acosf   n,n,m
  45
  46 MATHREF
  47  acosh, NAN,   arg,DOMAIN,EDOM
  48  acosh, < 1.0, NAN,DOMAIN,EDOM
  49  acosh, >=1.0, acosh(arg),,,
  50
  51 MATHREF
  52  acoshf, NAN,   arg,DOMAIN,EDOM
  53  acoshf, < 1.0, NAN,DOMAIN,EDOM
  54  acoshf, >=1.0, acosh(arg),,,
  55
  56 */
  57
  58 /*
  59  * ====================================================
  60  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  61  *
  62  * Developed at SunPro, a Sun Microsystems, Inc. business.
  63  * Permission to use, copy, modify, and distribute this
  64  * software is freely granted, provided that this notice
  65  * is preserved.
  66  * ====================================================
  67  *
  68  */
  69
  70 /* acosh(x)
  71  * Method :
  72  *      Based on
  73  *              acosh(x) = log [ x + sqrt(x*x-1) ]
  74  *      we have
  75  *              acosh(x) := log(x)+ln2, if x is large; else
  76  *              acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
  77  *              acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
  78  *
  79  * Special cases:
  80  *      acosh(x) is NaN with signal if x<1.
  81  *      acosh(NaN) is NaN without signal.
  82  */
  83
  84 #include "fdlibm.h"
  85
  86 #ifndef _DOUBLE_IS_32BITS
  87
  88 #ifdef __STDC__
  89 static const double
  90 #else
  91 static double
  92 #endif
  93 one     = 1.0,
  94 ln2     = 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
  95
  96 #ifdef __STDC__
  97         double acosh(double x)
  98 #else
  99         double acosh(x)
 100         double x;
 101 #endif
 102 {
 103         double t;
 104         __int32_t hx;
 105         __uint32_t lx;
 106         EXTRACT_WORDS(hx,lx,x);
 107         if(hx<0x3ff00000) {             /* x < 1 */
 108             return (x-x)/(x-x);
 109         } else if(hx >=0x41b00000) {    /* x > 2**28 */
 110             if(hx >=0x7ff00000) {       /* x is inf of NaN */
 111                 return x+x;
 112             } else
 113                 return log(x)+ln2;      /* acosh(huge)=log(2x) */
 114         } else if(((hx-0x3ff00000)|lx)==0) {
 115             return 0.0;                 /* acosh(1) = 0 */
 116         } else if (hx > 0x40000000) {   /* 2**28 > x > 2 */
 117             t=x*x;
 118             return log(2.0*x-one/(x+sqrt(t-one)));
 119         } else {                        /* 1<x<2 */
 120             t = x-one;
 121             return log1p(t+sqrt(2.0*t+t*t));
 122         }
 123 }
 124
 125 #endif /* defined(_DOUBLE_IS_32BITS) */