RTAI API: base/math/e_acosh.c Source File

00001 /* @(#)e_acosh.c 5.1 93/09/24 */ 00002 /* 00003 * ==================================================== 00004 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 00005 * 00006 * Developed at SunPro, a Sun Microsystems, Inc. business. 00007 * Permission to use, copy, modify, and distribute this 00008 * software is freely granted, provided that this notice 00009 * is preserved. 00010 * ==================================================== 00011 */ 00012 00013 #if defined(LIBM_SCCS) && !defined(lint) 00014 static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $"; 00015 #endif 00016 00017 /* __ieee754_acosh(x) 00018 * Method : 00019 * Based on 00020 * acosh(x) = log [ x + sqrt(x*x-1) ] 00021 * we have 00022 * acosh(x) := log(x)+ln2, if x is large; else 00023 * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else 00024 * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. 00025 * 00026 * Special cases: 00027 * acosh(x) is NaN with signal if x<1. 00028 * acosh(NaN) is NaN without signal. 00029 */ 00030 00031 #include "math.h" 00032 #include "mathP.h" 00033 00034 #ifdef __STDC__ 00035 static const double 00036 #else 00037 static double 00038 #endif 00039 one = 1.0, 00040 ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ 00041 00042 #ifdef __STDC__ 00043 double __ieee754_acosh(double x) 00044 #else 00045 double __ieee754_acosh(x) 00046 double x; 00047 #endif 00048 { 00049 double t; 00050 int32_t hx; 00051 u_int32_t lx; 00052 EXTRACT_WORDS(hx,lx,x); 00053 if(hx<0x3ff00000) { /* x < 1 */ 00054 return (x-x)/(x-x); 00055 } else if(hx >=0x41b00000) { /* x > 2**28 */ 00056 if(hx >=0x7ff00000) { /* x is inf of NaN */ 00057 return x+x; 00058 } else 00059 return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */ 00060 } else if(((hx-0x3ff00000)|lx)==0) { 00061 return 0.0; /* acosh(1) = 0 */ 00062 } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ 00063 t=x*x; 00064 return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one))); 00065 } else { /* 1<x<2 */ 00066 t = x-one; 00067 return log1p(t+sqrt(2.0*t+t*t)); 00068 } 00069 }