missing/acosh.c

DEFINITIONS

1. acosh
2. asinh
3. atanh

   1  /**********************************************************************
   2  
   3    acosh.c -
   4  
   5    $Author: nobu $
   6    $Date: 2002/04/12 03:24:41 $
   7    created at: Fri Apr 12 00:34:17 JST 2002
   8  
   9    public domain rewrite of acosh(3), asinh(3) and atanh(3)
  10  
  11  **********************************************************************/
  12  
  13  #include <errno.h>
  14  #include <float.h>
  15  #include <math.h>
  16  
  17  /* DBL_MANT_DIG must be less than 4 times of bits of int */
  18  #ifdef DBL_MANT_DIG
  19  #define DBL_MANT_DIG 53         /* in this case, at least 12 digit precision */
  20  #endif
  21  #define BIG_CRITERIA_BIT (1<<DBL_MANT_DIG/2)
  22  #if BIG_CRITERIA_BIT > 0
  23  #define BIG_CRITERIA (1.0*BIG_CRITERIA_BIT)
  24  #else
  25  #define BIG_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/2+1-DBL_MANT_DIG/4)))
  26  #endif
  27  #define SMALL_CRITERIA_BIT (1<<(DBL_MANT_DIG/3))
  28  #if SMALL_CRITERIA_BIT > 0
  29  #define SMALL_CRITERIA (1.0/SMALL_CRITERIA_BIT)
  30  #else
  31  #define SMALL_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/3+1-DBL_MANT_DIG/4)))
  32  #endif
  33  
  34  #ifndef HAVE_ACOSH
  35  double
acosh(x)
  37      double x;
  38  {
  39      if (x < 1)
  40          x = -1;                 /* NaN */
  41      else if (x == 1)
  42          return 0;
  43      else if (x > BIG_CRITERIA)
  44          x += x;
  45      else
  46          x += sqrt((x + 1) * (x - 1));
  47      return log(x);
  48  }
  49  #endif
  50  
  51  #ifndef HAVE_ASINH
  52  double
asinh(x)
  54      double x;
  55  {
  56      int neg = x < 0;
  57      double z = fabs(x);
  58  
  59      if (z < SMALL_CRITERIA) return x;
  60      if (z < (1.0/(1<<DBL_MANT_DIG/5))) {
  61          double x2 = z * z;
  62          z *= 1 + x2 * (-1.0/6.0 + x2 * 3.0/40.0);
  63      }
  64      else if (z > BIG_CRITERIA) {
  65          z = log(z + z);
  66      }
  67      else {
  68          z = log(z + sqrt(z * z + 1));
  69      }
  70      if (neg) z = -z;
  71      return z;
  72  }
  73  #endif
  74  
  75  #ifndef HAVE_ATANH
  76  double
atanh(x)
  78      double x;
  79  {
  80      int neg = x < 0;
  81      double z = fabs(x);
  82  
  83      if (z < SMALL_CRITERIA) return x;
  84      z = log(z > 1 ? -1 : (1 + z) / (1 - z)) / 2;
  85      if (neg) z = -z;
  86      return z;
  87  }
  88  #endif