lib/msun/src/e_acosl.c

   1
   2 /* @(#)e_acos.c 1.3 95/01/18 */
   3 /* FreeBSD: head/lib/msun/src/e_acos.c 176451 2008-02-22 02:30:36Z 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 #include <sys/cdefs.h>
  16 /*
  17  * See comments in e_acos.c.
  18  * Converted to long double by David Schultz <das@FreeBSD.ORG>.
  19  */
  20
  21 #include <float.h>
  22
  23 #include "invtrig.h"
  24 #include "math.h"
  25 #include "math_private.h"
  26
  27 static const long double
  28 one=  1.00000000000000000000e+00;
  29
  30 #ifdef __i386__
  31 /* XXX Work around the fact that gcc truncates long double constants on i386 */
  32 static volatile double
  33 pi1 =  3.14159265358979311600e+00,      /*  0x1.921fb54442d18p+1  */
  34 pi2 =  1.22514845490862001043e-16;      /*  0x1.1a80000000000p-53 */
  35 #define pi      ((long double)pi1 + pi2)
  36 #else
  37 static const long double
  38 pi =  3.14159265358979323846264338327950280e+00L;
  39 #endif
  40
  41 long double
  42 acosl(long double x)
  43 {
  44         union IEEEl2bits u;
  45         long double z,p,q,r,w,s,c,df;
  46         int16_t expsign, expt;
  47         u.e = x;
  48         expsign = u.xbits.expsign;
  49         expt = expsign & 0x7fff;
  50         if(expt >= BIAS) {      /* |x| >= 1 */
  51             if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0) {
  52                 if (expsign>0) return 0.0;      /* acos(1) = 0  */
  53                 else return pi+2.0*pio2_lo;     /* acos(-1)= pi */
  54             }
  55             return (x-x)/(x-x);         /* acos(|x|>1) is NaN */
  56         }
  57         if(expt<BIAS-1) {       /* |x| < 0.5 */
  58             if(expt<ACOS_CONST) return pio2_hi+pio2_lo;/*x tiny: acosl=pi/2*/
  59             z = x*x;
  60             p = P(z);
  61             q = Q(z);
  62             r = p/q;
  63             return pio2_hi - (x - (pio2_lo-x*r));
  64         } else  if (expsign<0) {        /* x < -0.5 */
  65             z = (one+x)*0.5;
  66             p = P(z);
  67             q = Q(z);
  68             s = sqrtl(z);
  69             r = p/q;
  70             w = r*s-pio2_lo;
  71             return pi - 2.0*(s+w);
  72         } else {                        /* x > 0.5 */
  73             z = (one-x)*0.5;
  74             s = sqrtl(z);
  75             u.e = s;
  76             u.bits.manl = 0;
  77             df = u.e;
  78             c  = (z-df*df)/(s+df);
  79             p = P(z);
  80             q = Q(z);
  81             r = p/q;
  82             w = r*s+c;
  83             return 2.0*(df+w);
  84         }
  85 }