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