1
2 /* @(#)z_pow.c 1.0 98/08/13 */
3
4 /*
5 FUNCTION
6 <<pow>>, <<powf>>---x to the power y
7 INDEX
8 pow
9 INDEX
10 powf
11
12
13 ANSI_SYNOPSIS
14 #include <math.h>
15 double pow(double <[x]>, double <[y]>);
16 float pow(float <[x]>, float <[y]>);
17
18 TRAD_SYNOPSIS
19 #include <math.h>
20 double pow(<[x]>, <[y]>);
21 double <[x]>, <[y]>;
22
23 float pow(<[x]>, <[y]>);
24 float <[x]>, <[y]>;
25
26 DESCRIPTION
27 <<pow>> and <<powf>> calculate <[x]> raised to the exponent <[y]>.
28 @tex
29 (That is, $x^y$.)
30 @end tex
31
32 RETURNS
33 On success, <<pow>> and <<powf>> return the value calculated.
34
35 When the argument values would produce overflow, <<pow>>
36 returns <<HUGE_VAL>> and set <<errno>> to <<ERANGE>>. If the
37 argument <[x]> passed to <<pow>> or <<powf>> is a negative
38 noninteger, and <[y]> is also not an integer, then <<errno>>
39 is set to <<EDOM>>. If <[x]> and <[y]> are both 0, then
40 <<pow>> and <<powf>> return <<1>>.
41
42 You can modify error handling for these functions using <<matherr>>.
43
44 PORTABILITY
45 <<pow>> is ANSI C. <<powf>> is an extension. */
46
47 #include <float.h>
48 #include "fdlibm.h"
49 #include "zmath.h"
50
51 #ifndef _DOUBLE_IS_32BITS
52
pow(double x,double y)53 double pow (double x, double y)
54 {
55 double d, k, t, r = 1.0;
56 int n, sign, exponent_is_even_int = 0;
57 __uint32_t px;
58
59 GET_HIGH_WORD (px, x);
60
61 k = modf (y, &d);
62
63 if (k == 0.0)
64 {
65 /* Exponent y is an integer. */
66 if (modf (ldexp (y, -1), &t))
67 {
68 /* y is odd. */
69 exponent_is_even_int = 0;
70 }
71 else
72 {
73 /* y is even. */
74 exponent_is_even_int = 1;
75 }
76 }
77
78 if (x == 0.0)
79 {
80 if (y <= 0.0)
81 errno = EDOM;
82 }
83 else if ((t = y * log (fabs (x))) >= BIGX)
84 {
85 errno = ERANGE;
86 if (px & 0x80000000)
87 {
88 /* x is negative. */
89 if (k)
90 {
91 /* y is not an integer. */
92 errno = EDOM;
93 x = 0.0;
94 }
95 else if (exponent_is_even_int)
96 x = z_infinity.d;
97 else
98 x = -z_infinity.d;
99 }
100 else
101 {
102 x = z_infinity.d;
103 }
104 }
105 else if (t < SMALLX)
106 {
107 errno = ERANGE;
108 x = 0.0;
109 }
110 else
111 {
112 if ( !k && fabs(d) <= 32767 )
113 {
114 n = (int) d;
115
116 if ((sign = (n < 0)))
117 n = -n;
118
119 while ( n > 0 )
120 {
121 if ((unsigned int) n % 2)
122 r *= x;
123 x *= x;
124 n = (unsigned int) n / 2;
125 }
126
127 if (sign)
128 r = 1.0 / r;
129
130 return r;
131 }
132 else
133 {
134 if ( px & 0x80000000 )
135 {
136 /* x is negative. */
137 if ( k )
138 {
139 /* y is not an integer. */
140 errno = EDOM;
141 return 0.0;
142 }
143 }
144
145 x = exp (t);
146
147 if (!exponent_is_even_int)
148 {
149 if (px & 0x80000000)
150 {
151 /* y is an odd integer, and x is negative,
152 so the result is negative. */
153 GET_HIGH_WORD (px, x);
154 px |= 0x80000000;
155 SET_HIGH_WORD (x, px);
156 }
157 }
158 }
159 }
160
161 return x;
162 }
163
164 #endif _DOUBLE_IS_32BITS
165