Ruby 1.9.3p327(2012-11-10revision37606)
|
00001 /********************************************************************** 00002 00003 acosh.c - 00004 00005 $Author: akr $ 00006 created at: Fri Apr 12 00:34:17 JST 2002 00007 00008 public domain rewrite of acosh(3), asinh(3) and atanh(3) 00009 00010 **********************************************************************/ 00011 00012 #include <errno.h> 00013 #include <float.h> 00014 #include <math.h> 00015 #include "ruby.h" 00016 00017 /* DBL_MANT_DIG must be less than 4 times of bits of int */ 00018 #ifndef DBL_MANT_DIG 00019 #define DBL_MANT_DIG 53 /* in this case, at least 12 digit precision */ 00020 #endif 00021 #define BIG_CRITERIA_BIT (1<<DBL_MANT_DIG/2) 00022 #if BIG_CRITERIA_BIT > 0 00023 #define BIG_CRITERIA (1.0*BIG_CRITERIA_BIT) 00024 #else 00025 #define BIG_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/2+1-DBL_MANT_DIG/4))) 00026 #endif 00027 #define SMALL_CRITERIA_BIT (1<<(DBL_MANT_DIG/3)) 00028 #if SMALL_CRITERIA_BIT > 0 00029 #define SMALL_CRITERIA (1.0/SMALL_CRITERIA_BIT) 00030 #else 00031 #define SMALL_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/3+1-DBL_MANT_DIG/4))) 00032 #endif 00033 00034 #ifndef HAVE_ACOSH 00035 double 00036 acosh(double x) 00037 { 00038 if (x < 1) 00039 x = -1; /* NaN */ 00040 else if (x == 1) 00041 return 0; 00042 else if (x > BIG_CRITERIA) 00043 x += x; 00044 else 00045 x += sqrt((x + 1) * (x - 1)); 00046 return log(x); 00047 } 00048 #endif 00049 00050 #ifndef HAVE_ASINH 00051 double 00052 asinh(double x) 00053 { 00054 int neg = x < 0; 00055 double z = fabs(x); 00056 00057 if (z < SMALL_CRITERIA) return x; 00058 if (z < (1.0/(1<<DBL_MANT_DIG/5))) { 00059 double x2 = z * z; 00060 z *= 1 + x2 * (-1.0/6.0 + x2 * 3.0/40.0); 00061 } 00062 else if (z > BIG_CRITERIA) { 00063 z = log(z + z); 00064 } 00065 else { 00066 z = log(z + sqrt(z * z + 1)); 00067 } 00068 if (neg) z = -z; 00069 return z; 00070 } 00071 #endif 00072 00073 #ifndef HAVE_ATANH 00074 double 00075 atanh(double x) 00076 { 00077 int neg = x < 0; 00078 double z = fabs(x); 00079 00080 if (z < SMALL_CRITERIA) return x; 00081 z = log(z > 1 ? -1 : (1 + z) / (1 - z)) / 2; 00082 if (neg) z = -z; 00083 if (isinf(z)) 00084 #if defined(ERANGE) 00085 errno = ERANGE; 00086 #elif defined(EDOM) 00087 errno = EDOM; 00088 #else 00089 ; 00090 #endif 00091 return z; 00092 } 00093 #endif 00094