compiler/stdc/math/e_asinl.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 #ifndef lint
  13 static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_asinl.c,v 1.2 2008/08/03 17:49:05 das Exp $";
  14 #endif
  15
  16 /*
  17  * See comments in e_asin.c.
  18  * Converted to long double by David Schultz <das@FreeBSD.ORG>.
  19  */
  20
  21 #include <float.h>
  22 #include "math.h"
  23
  24 #include "invtrig.h"
  25 #include "math_private.h"
  26
  27 static const long double
  28 one =  1.00000000000000000000e+00,
  29 huge = 1.000e+300;
  30
  31 long double
  32 asinl(long double x)
  33 {
  34         union IEEEl2bits u;
  35         long double t=0.0,w,p,q,c,r,s;
  36         int16_t expsign, expt;
  37         u.e = x;
  38         expsign = u.xbits.expsign;
  39         expt = expsign & 0x7fff;
  40         if(expt >= BIAS) {              /* |x|>= 1 */
  41                 if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0)
  42                     /* asin(1)=+-pi/2 with inexact */
  43                     return x*pio2_hi+x*pio2_lo;
  44             return (x-x)/(x-x);         /* asin(|x|>1) is NaN */
  45         } else if (expt<BIAS-1) {       /* |x|<0.5 */
  46             if(expt<ASIN_LINEAR) {      /* if |x| is small, asinl(x)=x */
  47                 if(huge+x>one) return x;/* return x with inexact if x!=0*/
  48             }
  49             t = x*x;
  50             p = P(t);
  51             q = Q(t);
  52             w = p/q;
  53             return x+x*w;
  54         }
  55         /* 1> |x|>= 0.5 */
  56         w = one-fabsl(x);
  57         t = w*0.5;
  58         p = P(t);
  59         q = Q(t);
  60         s = sqrtl(t);
  61         if(u.bits.manh>=THRESH) {       /* if |x| is close to 1 */
  62             w = p/q;
  63             t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
  64         } else {
  65             u.e = s;
  66             u.bits.manl = 0;
  67             w = u.e;
  68             c  = (t-w*w)/(s+w);
  69             r  = p/q;
  70             p  = 2.0*s*r-(pio2_lo-2.0*c);
  71             q  = pio4_hi-2.0*w;
  72             t  = pio4_hi-(p-q);
  73         }
  74         if(expsign>0) return t; else return -t;
  75 }