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

00001 /* @(#)s_asinh.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: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $";
00015 #endif
00016 
00017 /* asinh(x)
00018  * Method :
00019  *  Based on 
00020  *      asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
00021  *  we have
00022  *  asinh(x) := x  if  1+x*x=1,
00023  *       := sign(x)*(log(x)+ln2)) for large |x|, else
00024  *       := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
00025  *       := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))  
00026  */
00027 
00028 #include "math.h"
00029 #include "mathP.h"
00030 
00031 #ifdef __STDC__
00032 static const double 
00033 #else
00034 static double 
00035 #endif
00036 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
00037 ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
00038 huge=  1.00000000000000000000e+300; 
00039 
00040 #ifdef __STDC__
00041     double asinh(double x)
00042 #else
00043     double asinh(x)
00044     double x;
00045 #endif
00046 {   
00047     double t,w;
00048     int32_t hx,ix;
00049     GET_HIGH_WORD(hx,x);
00050     ix = hx&0x7fffffff;
00051     if(ix>=0x7ff00000) return x+x;  /* x is inf or NaN */
00052     if(ix< 0x3e300000) {    /* |x|<2**-28 */
00053         if(huge+x>one) return x;    /* return x inexact except 0 */
00054     } 
00055     if(ix>0x41b00000) { /* |x| > 2**28 */
00056         w = __ieee754_log(fabs(x))+ln2;
00057     } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
00058         t = fabs(x);
00059         w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
00060     } else {        /* 2.0 > |x| > 2**-28 */
00061         t = x*x;
00062         w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
00063     }
00064     if(hx>0) return w; else return -w;
00065 }