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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
/* @(#)z_powf.c 1.0 98/08/13 */
#include <float.h>
#include "fdlibm.h"
#include "zmath.h"
float powf (float x, float y)
{
float d, k, t, r = 1.0;
int n, sign, exponent_is_even_int = 0;
__int32_t px;
GET_FLOAT_WORD (px, x);
k = modff (y, &d);
if (k == 0.0)
{
/* Exponent y is an integer. */
if (modff (ldexpf (y, -1), &t))
{
/* y is odd. */
exponent_is_even_int = 0;
}
else
{
/* y is even. */
exponent_is_even_int = 1;
}
}
if (x == 0.0 && y <= 0.0)
{
errno = EDOM;
}
else if ((t = y * log (fabsf (x))) >= BIGX)
{
errno = ERANGE;
if (px & 0x80000000)
{
/* x is negative. */
if (k)
{
/* y is not an integer. */
errno = EDOM;
x = 0.0;
}
else if (exponent_is_even_int)
x = z_infinity_f.f;
else
x = -z_infinity_f.f;
}
else
{
x = z_infinity_f.f;
}
}
else if (t < SMALLX)
{
errno = ERANGE;
x = 0.0;
}
else
{
if ( !k && fabsf (d) <= 32767 )
{
n = (int) d;
if ((sign = (n < 0)))
n = -n;
while ( n > 0 )
{
if ((unsigned int) n % 2)
r *= x;
x *= x;
n = (unsigned int) n / 2;
}
if (sign)
r = 1.0 / r;
return r;
}
else
{
if ( px & 0x80000000 )
{
/* x is negative. */
if (k)
{
/* y is not an integer. */
errno = EDOM;
return 0.0;
}
}
x = exp (t);
if (!exponent_is_even_int)
{
if (px & 0x80000000)
{
/* y is an odd integer, and x is negative,
so the result is negative. */
GET_FLOAT_WORD (px, x);
px |= 0x80000000;
SET_FLOAT_WORD (x, px);
}
}
}
}
return x;
}
#ifdef _DOUBLE_IS_32BITS
double pow (double x, double y)
{
return (double) powf ((float) x, (float) y);
}
#endif /* defined(_DOUBLE_IS_32BITS) */
|
