lib/libm/src/e_acosl.c

   1
   2 /* @(#)e_acos.c 1.3 95/01/18 */
   3 /* $FreeBSD: head/lib/msun/src/e_acosl.c 181152 2008-08-02 03:56:22Z das $ */
   4 /*
   5  * ====================================================
   6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   7  *
   8  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   9  * Permission to use, copy, modify, and distribute this
  10  * software is freely granted, provided that this notice
  11  * is preserved.
  12  * ====================================================
  13  */
  14
  15 /*
  16  * See comments in e_acos.c.
  17  * Converted to long double by David Schultz <das@FreeBSD.ORG>.
  18  */
  19
  20 #include <float.h>
  21
  22 #include "invtrig.h"
  23 #include "math.h"
  24 #include "math_private.h"
  25
  26 static const long double
  27 one=  1.00000000000000000000e+00;
  28
  29 #ifdef __i386__
  30 /* XXX Work around the fact that gcc truncates long double constants on i386 */
  31 static volatile double
  32 pi1 =  3.14159265358979311600e+00,      /*  0x1.921fb54442d18p+1  */
  33 pi2 =  1.22514845490862001043e-16;      /*  0x1.1a80000000000p-53 */
  34 #define pi      ((long double)pi1 + pi2)
  35 #else
  36 static const long double
  37 pi =  3.14159265358979323846264338327950280e+00L;
  38 #endif
  39
  40 long double
  41 acosl(long double x)
  42 {
  43         union IEEEl2bits u;
  44         long double z,p,q,r,w,s,c,df;
  45         int16_t expsign, expt;
  46         u.e = x;
  47         expsign = u.xbits.expsign;
  48         expt = expsign & 0x7fff;
  49         if(expt >= BIAS) {      /* |x| >= 1 */
  50             if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0) {
  51                 if (expsign>0) return 0.0;      /* acos(1) = 0  */
  52                 else return pi+2.0*pio2_lo;     /* acos(-1)= pi */
  53             }
  54             return (x-x)/(x-x);         /* acos(|x|>1) is NaN */
  55         }
  56         if(expt<BIAS-1) {       /* |x| < 0.5 */
  57             if(expt<ACOS_CONST) return pio2_hi+pio2_lo;/*x tiny: acosl=pi/2*/
  58             z = x*x;
  59             p = P(z);
  60             q = Q(z);
  61             r = p/q;
  62             return pio2_hi - (x - (pio2_lo-x*r));
  63         } else  if (expsign<0) {        /* x < -0.5 */
  64             z = (one+x)*0.5;
  65             p = P(z);
  66             q = Q(z);
  67             s = sqrtl(z);
  68             r = p/q;
  69             w = r*s-pio2_lo;
  70             return pi - 2.0*(s+w);
  71         } else {                        /* x > 0.5 */
  72             z = (one-x)*0.5;
  73             s = sqrtl(z);
  74             u.e = s;
  75             u.bits.manl = 0;
  76             df = u.e;
  77             c  = (z-df*df)/(s+df);
  78             p = P(z);
  79             q = Q(z);
  80             r = p/q;
  81             w = r*s+c;
  82             return 2.0*(df+w);
  83         }
  84 }