sysdeps/ieee754/flt-32/e_asinf.c

   1 /* e_asinf.c -- float version of e_asin.c.
   2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   3  */
   4
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 /*
  17   Modifications for single precision expansion are
  18   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  19   and are incorporated herein by permission of the author.  The author
  20   reserves the right to distribute this material elsewhere under different
  21   copying permissions.  These modifications are distributed here under
  22   the following terms:
  23
  24     This library is free software; you can redistribute it and/or
  25     modify it under the terms of the GNU Lesser General Public
  26     License as published by the Free Software Foundation; either
  27     version 2.1 of the License, or (at your option) any later version.
  28
  29     This library is distributed in the hope that it will be useful,
  30     but WITHOUT ANY WARRANTY; without even the implied warranty of
  31     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  32     Lesser General Public License for more details.
  33
  34     You should have received a copy of the GNU Lesser General Public
  35     License along with this library; if not, write to the Free Software
  36     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  37
  38 #if defined(LIBM_SCCS) && !defined(lint)
  39 static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $";
  40 #endif
  41
  42 #include "math.h"
  43 #include "math_private.h"
  44
  45 #ifdef __STDC__
  46 static const float
  47 #else
  48 static float
  49 #endif
  50 one =  1.0000000000e+00, /* 0x3F800000 */
  51 huge =  1.000e+30,
  52
  53 pio2_hi = 1.57079637050628662109375f,
  54 pio2_lo = -4.37113900018624283e-8f,
  55 pio4_hi = 0.785398185253143310546875f,
  56
  57 /* asin x = x + x^3 p(x^2)
  58    -0.5 <= x <= 0.5;
  59    Peak relative error 4.8e-9 */
  60 p0 = 1.666675248e-1f,
  61 p1 = 7.495297643e-2f,
  62 p2 = 4.547037598e-2f,
  63 p3 = 2.417951451e-2f,
  64 p4 = 4.216630880e-2f;
  65
  66 #ifdef __STDC__
  67         float __ieee754_asinf(float x)
  68 #else
  69         float __ieee754_asinf(x)
  70         float x;
  71 #endif
  72 {
  73         float t,w,p,q,c,r,s;
  74         int32_t hx,ix;
  75         GET_FLOAT_WORD(hx,x);
  76         ix = hx&0x7fffffff;
  77         if(ix==0x3f800000) {
  78                 /* asin(1)=+-pi/2 with inexact */
  79             return x*pio2_hi+x*pio2_lo;
  80         } else if(ix> 0x3f800000) {     /* |x|>= 1 */
  81             return (x-x)/(x-x);         /* asin(|x|>1) is NaN */
  82         } else if (ix<0x3f000000) {     /* |x|<0.5 */
  83             if(ix<0x32000000) {         /* if |x| < 2**-27 */
  84                 if(huge+x>one) return x;/* return x with inexact if x!=0*/
  85             } else {
  86                 t = x*x;
  87                 w = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
  88                 return x+x*w;
  89             }
  90         }
  91         /* 1> |x|>= 0.5 */
  92         w = one-fabsf(x);
  93         t = w*0.5f;
  94         p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
  95         s = __ieee754_sqrtf(t);
  96         if(ix>=0x3F79999A) {    /* if |x| > 0.975 */
  97             t = pio2_hi-(2.0f*(s+s*p)-pio2_lo);
  98         } else {
  99             int32_t iw;
 100             w  = s;
 101             GET_FLOAT_WORD(iw,w);
 102             SET_FLOAT_WORD(w,iw&0xfffff000);
 103             c  = (t-w*w)/(s+w);
 104             r  = p;
 105             p  = 2.0f*s*r-(pio2_lo-2.0f*c);
 106             q  = pio4_hi-2.0f*w;
 107             t  = pio4_hi-(p-q);
 108         }
 109         if(hx>0) return t; else return -t;
 110 }