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