1
2 /* @(#)z_powf.c 1.0 98/08/13 */
3 #include <float.h>
4 #include "fdlibm.h"
5 #include "zmath.h"
6
powf(float x,float y)7 float powf (float x, float y)
8 {
9 float d, k, t, r = 1.0;
10 int n, sign, exponent_is_even_int = 0;
11 __int32_t px;
12
13 GET_FLOAT_WORD (px, x);
14
15 k = modff (y, &d);
16
17 if (k == 0.0)
18 {
19 /* Exponent y is an integer. */
20 if (modff (ldexpf (y, -1), &t))
21 {
22 /* y is odd. */
23 exponent_is_even_int = 0;
24 }
25 else
26 {
27 /* y is even. */
28 exponent_is_even_int = 1;
29 }
30 }
31
32 if (x == 0.0 && y <= 0.0)
33 {
34 errno = EDOM;
35 }
36 else if ((t = y * log (fabsf (x))) >= BIGX)
37 {
38 errno = ERANGE;
39 if (px & 0x80000000)
40 {
41 /* x is negative. */
42 if (k)
43 {
44 /* y is not an integer. */
45 errno = EDOM;
46 x = 0.0;
47 }
48 else if (exponent_is_even_int)
49 x = z_infinity_f.f;
50 else
51 x = -z_infinity_f.f;
52 }
53 else
54 {
55 x = z_infinity_f.f;
56 }
57 }
58 else if (t < SMALLX)
59 {
60 errno = ERANGE;
61 x = 0.0;
62 }
63 else
64 {
65 if ( !k && fabsf (d) <= 32767 )
66 {
67 n = (int) d;
68
69 if ((sign = (n < 0)))
70 n = -n;
71
72 while ( n > 0 )
73 {
74 if ((unsigned int) n % 2)
75 r *= x;
76 x *= x;
77 n = (unsigned int) n / 2;
78 }
79
80 if (sign)
81 r = 1.0 / r;
82
83 return r;
84 }
85 else
86 {
87 if ( px & 0x80000000 )
88 {
89 /* x is negative. */
90 if (k)
91 {
92 /* y is not an integer. */
93 errno = EDOM;
94 return 0.0;
95 }
96 }
97
98 x = exp (t);
99
100 if (!exponent_is_even_int)
101 {
102 if (px & 0x80000000)
103 {
104 /* y is an odd integer, and x is negative,
105 so the result is negative. */
106 GET_FLOAT_WORD (px, x);
107 px |= 0x80000000;
108 SET_FLOAT_WORD (x, px);
109 }
110 }
111 }
112 }
113
114 return x;
115 }
116
117 #ifdef _DOUBLE_IS_32BITS
118
pow(double x,double y)119 double pow (double x, double y)
120 {
121 return (double) powf ((float) x, (float) y);
122 }
123
124 #endif /* defined(_DOUBLE_IS_32BITS) */
125