1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
|
/* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "libm.h" // for GET_FLOAT_WORD, SET_FLOAT_WORD
#include <math.h> // for sqrtf, acosf, float_t
#include <stdint.h> // for uint32_t
static const float pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
pio2_lo = 7.5497894159e-08, /* 0x33a22168 */
pS0 = 1.6666586697e-01, pS1 = -4.2743422091e-02, pS2 = -8.6563630030e-03, qS1 = -7.0662963390e-01;
static float R(float z)
{
float_t p, q;
p = z * (pS0 + z * (pS1 + z * pS2));
q = 1.0f + z * qS1;
return p / q;
}
float acosf(float x)
{
float z, w, s, c, df;
uint32_t hx, ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
/* |x| >= 1 or nan */
if (ix >= 0x3f800000) {
if (ix == 0x3f800000) {
if (hx >> 31)
return 2 * pio2_hi + 0x1p-120f;
return 0;
}
return 0 / (x - x);
}
/* |x| < 0.5 */
if (ix < 0x3f000000) {
if (ix <= 0x32800000) /* |x| < 2**-26 */
return pio2_hi + 0x1p-120f;
return pio2_hi - (x - (pio2_lo - x * R(x * x)));
}
/* x < -0.5 */
if (hx >> 31) {
z = (1 + x) * 0.5f;
s = sqrtf(z);
w = R(z) * s - pio2_lo;
return 2 * (pio2_hi - (s + w));
}
/* x > 0.5 */
z = (1 - x) * 0.5f;
s = sqrtf(z);
GET_FLOAT_WORD(hx, s);
SET_FLOAT_WORD(df, hx & 0xfffff000);
c = (z - df * df) / (s + df);
w = R(z) * s + c;
return 2 * (df + w);
}
|
