missing/acosh.c

   1 /**********************************************************************
   2
   3   acosh.c -
   4
   5   $Author$
   6   created at: Fri Apr 12 00:34:17 JST 2002
   7
   8   public domain rewrite of acosh(3), asinh(3) and atanh(3)
   9
  10 **********************************************************************/
  11
  12 #include <errno.h>
  13 #include <float.h>
  14 #include <math.h>
  15 #include "ruby.h"
  16
  17 /* DBL_MANT_DIG must be less than 4 times of bits of int */
  18 #ifndef 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
  36 acosh(double x)
  37 {
  38     if (x < 1)
  39         x = -1;                 /* NaN */
  40     else if (x == 1)
  41         return 0;
  42     else if (x > BIG_CRITERIA)
  43         x += x;
  44     else
  45         x += sqrt((x + 1) * (x - 1));
  46     return log(x);
  47 }
  48 #endif
  49
  50 #ifndef HAVE_ASINH
  51 double
  52 asinh(double x)
  53 {
  54     int neg = x < 0;
  55     double z = fabs(x);
  56
  57     if (z < SMALL_CRITERIA) return x;
  58     if (z < (1.0/(1<<DBL_MANT_DIG/5))) {
  59         double x2 = z * z;
  60         z *= 1 + x2 * (-1.0/6.0 + x2 * 3.0/40.0);
  61     }
  62     else if (z > BIG_CRITERIA) {
  63         z = log(z + z);
  64     }
  65     else {
  66         z = log(z + sqrt(z * z + 1));
  67     }
  68     if (neg) z = -z;
  69     return z;
  70 }
  71 #endif
  72
  73 #ifndef HAVE_ATANH
  74 double
  75 atanh(double x)
  76 {
  77     int neg = x < 0;
  78     double z = fabs(x);
  79
  80     if (z < SMALL_CRITERIA) return x;
  81     z = log(z > 1 ? -1 : (1 + z) / (1 - z)) / 2;
  82     if (neg) z = -z;
  83     if (isinf(z))
  84 #if defined(ERANGE)
  85         errno = ERANGE;
  86 #elif defined(EDOM)
  87         errno = EDOM;
  88 #else
  89         ;
  90 #endif
  91     return z;
  92 }
  93 #endif