compiler/stdc/math/ld80/e_acoshl.c

   1 /* @(#)e_acosh.c 5.1 93/09/24 */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12
  13 /* acoshl(x)
  14  * Method :
  15  *      Based on
  16  *              acoshl(x) = logl [ x + sqrtl(x*x-1) ]
  17  *      we have
  18  *              acoshl(x) := logl(x)+ln2,       if x is large; else
  19  *              acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
  20  *              acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
  21  *
  22  * Special cases:
  23  *      acoshl(x) is NaN with signal if x<1.
  24  *      acoshl(NaN) is NaN without signal.
  25  */
  26
  27 #include "math.h"
  28
  29 #include "math_private.h"
  30
  31 static const long double
  32 one     = 1.0,
  33 ln2     = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
  34
  35 long double
  36 acoshl(long double x)
  37 {
  38         long double t;
  39         uint32_t se,i0,i1;
  40         GET_LDOUBLE_WORDS(se,i0,i1,x);
  41         if(se<0x3fff || se & 0x8000) {  /* x < 1 */
  42             return (x-x)/(x-x);
  43         } else if(se >=0x401d) {        /* x > 2**30 */
  44             if(se >=0x7fff) {           /* x is inf of NaN */
  45                 return x+x;
  46             } else
  47                 return logl(x)+ln2;     /* acoshl(huge)=logl(2x) */
  48         } else if(((se-0x3fff)|i0|i1)==0) {
  49             return 0.0;                 /* acosh(1) = 0 */
  50         } else if (se > 0x4000) {       /* 2**28 > x > 2 */
  51             t=x*x;
  52             return logl(2.0*x-one/(x+sqrtl(t-one)));
  53         } else {                        /* 1<x<2 */
  54             t = x-one;
  55             return log1pl(t+sqrtl(2.0*t+t*t));
  56         }
  57 }