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