POK
|
00001 /* 00002 * POK header 00003 * 00004 * The following file is a part of the POK project. Any modification should 00005 * made according to the POK licence. You CANNOT use this file or a part of 00006 * this file is this part of a file for your own project 00007 * 00008 * For more information on the POK licence, please see our LICENCE FILE 00009 * 00010 * Please follow the coding guidelines described in doc/CODING_GUIDELINES 00011 * 00012 * Copyright (c) 2007-2009 POK team 00013 * 00014 * Created by julien on Fri Jan 30 14:41:34 2009 00015 */ 00016 00017 /* e_sinhf.c -- float version of e_sinh.c. 00018 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 00019 */ 00020 00021 /* 00022 * ==================================================== 00023 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 00024 * 00025 * Developed at SunPro, a Sun Microsystems, Inc. business. 00026 * Permission to use, copy, modify, and distribute this 00027 * software is freely granted, provided that this notice 00028 * is preserved. 00029 * ==================================================== 00030 */ 00031 00032 #ifdef POK_NEEDS_LIBMATH 00033 00034 #include <libm.h> 00035 #include "math_private.h" 00036 00037 static const float one = 1.0, shuge = 1.0e37; 00038 00039 float 00040 __ieee754_sinhf(float x) 00041 { 00042 float t,w,h; 00043 int32_t ix,jx; 00044 00045 GET_FLOAT_WORD(jx,x); 00046 ix = jx&0x7fffffff; 00047 00048 /* x is INF or NaN */ 00049 if(ix>=0x7f800000) return x+x; 00050 00051 h = 0.5; 00052 if (jx<0) h = -h; 00053 /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ 00054 if (ix < 0x41b00000) { /* |x|<22 */ 00055 if (ix<0x31800000) /* |x|<2**-28 */ 00056 if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ 00057 t = expm1f(fabsf(x)); 00058 if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one)); 00059 return h*(t+t/(t+one)); 00060 } 00061 00062 /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ 00063 if (ix < 0x42b17180) return h*__ieee754_expf(fabsf(x)); 00064 00065 /* |x| in [log(maxdouble), overflowthresold] */ 00066 if (ix<=0x42b2d4fc) { 00067 w = __ieee754_expf((float)0.5*fabsf(x)); 00068 t = h*w; 00069 return t*w; 00070 } 00071 00072 /* |x| > overflowthresold, sinh(x) overflow */ 00073 return x*shuge; 00074 } 00075 #endif 00076