Back to index

glibc  2.9
e_asinf.c
Go to the documentation of this file.
00001 /* e_asinf.c -- float version of e_asin.c.
00002  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 /*
00017   Modifications for single precision expansion are 
00018   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00019   and are incorporated herein by permission of the author.  The author 
00020   reserves the right to distribute this material elsewhere under different
00021   copying permissions.  These modifications are distributed here under 
00022   the following terms:
00023 
00024     This library is free software; you can redistribute it and/or
00025     modify it under the terms of the GNU Lesser General Public
00026     License as published by the Free Software Foundation; either
00027     version 2.1 of the License, or (at your option) any later version.
00028 
00029     This library is distributed in the hope that it will be useful,
00030     but WITHOUT ANY WARRANTY; without even the implied warranty of
00031     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00032     Lesser General Public License for more details.
00033 
00034     You should have received a copy of the GNU Lesser General Public
00035     License along with this library; if not, write to the Free Software
00036     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00037 
00038 #if defined(LIBM_SCCS) && !defined(lint)
00039 static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $";
00040 #endif
00041 
00042 #include "math.h"
00043 #include "math_private.h"
00044 
00045 #ifdef __STDC__
00046 static const float
00047 #else
00048 static float
00049 #endif
00050 one =  1.0000000000e+00, /* 0x3F800000 */
00051 huge =  1.000e+30,
00052 
00053 pio2_hi = 1.57079637050628662109375f,
00054 pio2_lo = -4.37113900018624283e-8f,
00055 pio4_hi = 0.785398185253143310546875f,
00056 
00057 /* asin x = x + x^3 p(x^2)
00058    -0.5 <= x <= 0.5;
00059    Peak relative error 4.8e-9 */
00060 p0 = 1.666675248e-1f,
00061 p1 = 7.495297643e-2f,
00062 p2 = 4.547037598e-2f,
00063 p3 = 2.417951451e-2f,
00064 p4 = 4.216630880e-2f;
00065 
00066 #ifdef __STDC__
00067        float __ieee754_asinf(float x)
00068 #else
00069        float __ieee754_asinf(x)
00070        float x;
00071 #endif
00072 {
00073        float t,w,p,q,c,r,s;
00074        int32_t hx,ix;
00075        GET_FLOAT_WORD(hx,x);
00076        ix = hx&0x7fffffff;
00077        if(ix==0x3f800000) {
00078               /* asin(1)=+-pi/2 with inexact */
00079            return x*pio2_hi+x*pio2_lo;
00080        } else if(ix> 0x3f800000) { /* |x|>= 1 */
00081            return (x-x)/(x-x);            /* asin(|x|>1) is NaN */
00082        } else if (ix<0x3f000000) { /* |x|<0.5 */
00083            if(ix<0x32000000) {            /* if |x| < 2**-27 */
00084               if(huge+x>one) return x;/* return x with inexact if x!=0*/
00085            } else {
00086               t = x*x;
00087               w = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
00088               return x+x*w;
00089            }
00090        }
00091        /* 1> |x|>= 0.5 */
00092        w = one-fabsf(x);
00093        t = w*0.5f;
00094        p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
00095        s = __ieee754_sqrtf(t);
00096        if(ix>=0x3F79999A) {        /* if |x| > 0.975 */
00097            t = pio2_hi-(2.0f*(s+s*p)-pio2_lo);
00098        } else {
00099            int32_t iw;
00100            w  = s;
00101            GET_FLOAT_WORD(iw,w);
00102            SET_FLOAT_WORD(w,iw&0xfffff000);
00103            c  = (t-w*w)/(s+w);
00104            r  = p;
00105            p  = 2.0f*s*r-(pio2_lo-2.0f*c);
00106            q  = pio4_hi-2.0f*w;
00107            t  = pio4_hi-(p-q);
00108        }
00109        if(hx>0) return t; else return -t;
00110 }