1//! `x^n` with fractional `n` approximation for a single-precision float.
2
3use super::F32;
4
5impl F32 {
6 /// Approximates a number raised to a floating point power.
7 pub fn powf(self, n: Self) -> Self {
8 // using x^n = exp(ln(x^n)) = exp(n*ln(x))
9 if self >= Self::ZERO {
10 (n * self.ln()).exp()
11 } else if !n.is_integer() {
12 Self::NAN
13 } else if n.is_even() {
14 // if n is even, then we know that the result will have no sign, so we can remove it
15 (n * self.without_sign().ln()).exp()
16 } else {
17 // if n isn't even, we need to multiply by -1.0 at the end.
18 -(n * self.without_sign().ln()).exp()
19 }
20 }
21}
22
23#[cfg(test)]
24mod tests {
25 use super::F32;
26
27 /// error builds up from both exp and ln approximation, so we double the error allowed.
28 pub(crate) const MAX_ERROR: f32 = 0.002;
29
30 /// powf(3,x) test vectors - `(input, output)`
31 pub(crate) const TEST_VECTORS_POW3: &[(f32, f32)] = &[
32 (-1e-20, 1.0),
33 (-1e-19, 1.0),
34 (-1e-18, 1.0),
35 (-1e-17, 1.0),
36 (-1e-16, 0.9999999999999999),
37 (-1e-15, 0.9999999999999989),
38 (-1e-14, 0.999999999999989),
39 (-1e-13, 0.9999999999998901),
40 (-1e-12, 0.9999999999989014),
41 (-1e-11, 0.9999999999890139),
42 (-1e-10, 0.9999999998901388),
43 (-1e-09, 0.9999999989013877),
44 (-1e-08, 0.9999999890138772),
45 (-1e-07, 0.999_999_9),
46 (-1e-06, 0.999_998_9),
47 (-1e-05, 0.999_989_03),
48 (-1e-04, 0.999_890_15),
49 (-0.001, 0.998_901_96),
50 (-0.01, 0.989_074),
51 (-0.1, 0.895_958_5),
52 (-1.0, 0.333_333_34),
53 (-10.0, 1.693_508_8e-5),
54 (-100.0, 0e0),
55 (-1000.0, 0.0),
56 (1e-20, 1.0),
57 (1e-19, 1.0),
58 (1e-18, 1.0),
59 (1e-17, 1.0),
60 (1e-16, 1.0),
61 (1e-15, 1.000000000000001),
62 (1e-14, 1.0000000000000109),
63 (1e-13, 1.00000000000011),
64 (1e-12, 1.0000000000010987),
65 (1e-11, 1.000000000010986),
66 (1e-10, 1.0000000001098612),
67 (1e-09, 1.0000000010986123),
68 (1e-08, 1.000000010986123),
69 (1e-07, 1.000_000_1),
70 (1e-06, 1.000_001_1),
71 (1e-05, 1.000_011),
72 (1e-04, 1.000_109_9),
73 (0.001, 1.001_099_2),
74 (0.01, 1.011_046_6),
75 (0.1, 1.116_123_2),
76 (1.0, 3.0),
77 (10.0, 59049.0),
78 ];
79
80 /// powf(150,x) test vectors - `(input, output)`
81 pub(crate) const TEST_VECTORS_POW150: &[(f32, f32)] = &[
82 (-1e-20, 1.0),
83 (-1e-19, 1.0),
84 (-1e-18, 1.0),
85 (-1e-17, 1.0),
86 (-1e-16, 0.9999999999999994),
87 (-1e-15, 0.999999999999995),
88 (-1e-14, 0.9999999999999499),
89 (-1e-13, 0.999999999999499),
90 (-1e-12, 0.9999999999949893),
91 (-1e-11, 0.9999999999498936),
92 (-1e-10, 0.9999999994989365),
93 (-1e-09, 0.9999999949893649),
94 (-1e-08, 0.999_999_94),
95 (-1e-07, 0.999_999_5),
96 (-1e-06, 0.999_995),
97 (-1e-05, 0.999_949_9),
98 (-1e-04, 0.999_499_1),
99 (-0.001, 0.995_001_9),
100 (-0.01, 0.951_128_24),
101 (-0.1, 0.605_885_9),
102 (-1.0, 0.006_666_667),
103 (-10.0, 1.734_153e-22),
104 (-100.0, 0e0),
105 (-1000.0, 0.0),
106 (-10000.0, 0.0),
107 (-100000.0, 0.0),
108 (-1000000.0, 0.0),
109 (-10000000.0, 0.0),
110 (-100000000.0, 0.0),
111 (-1000000000.0, 0.0),
112 (-10000000000.0, 0.0),
113 (-100000000000.0, 0.0),
114 (-1000000000000.0, 0.0),
115 (-10000000000000.0, 0.0),
116 (-100000000000000.0, 0.0),
117 (-1000000000000000.0, 0.0),
118 (-1e+16, 0.0),
119 (-1e+17, 0.0),
120 (-1e+18, 0.0),
121 (-1e+19, 0.0),
122 (1e-20, 1.0),
123 (1e-19, 1.0),
124 (1e-18, 1.0),
125 (1e-17, 1.0),
126 (1e-16, 1.0000000000000004),
127 (1e-15, 1.000000000000005),
128 (1e-14, 1.0000000000000502),
129 (1e-13, 1.0000000000005012),
130 (1e-12, 1.0000000000050107),
131 (1e-11, 1.0000000000501064),
132 (1e-10, 1.0000000005010636),
133 (1e-09, 1.0000000050106352),
134 (1e-08, 1.000000050106354),
135 (1e-07, 1.000_000_5),
136 (1e-06, 1.000_005),
137 (1e-05, 1.000_050_1),
138 (1e-04, 1.000_501_2),
139 (0.001, 1.005_023_2),
140 (0.01, 1.051_382_9),
141 (0.1, 1.650_475_6),
142 (1.0, 150.0),
143 (10.0, 5.766_504e21),
144 ];
145
146 /// misc powf(x,n) test vectors - `(base_input, power_input, output)`
147 pub(crate) const TEST_VECTORS_MISC: &[(f32, f32, f32)] = &[
148 (-0.5881598, 2.0, 0.345_931_95),
149 (-0.5881598, 3.2, f32::NAN),
150 (-0.5881598, 3.0, -0.203_463_27),
151 (-1000000.0, 4.0, 1e+24),
152 ];
153
154 fn calc_relative_error(experimental: F32, expected: f32) -> F32 {
155 if experimental.is_nan() && expected.is_nan() {
156 F32::ZERO
157 } else if expected != 0.0 {
158 (experimental - expected) / expected
159 } else {
160 (experimental - expected) / (expected + 1.0e-20)
161 }
162 }
163
164 #[test]
165 fn sanity_check() {
166 for &(x, expected) in TEST_VECTORS_POW3 {
167 let exp_x = F32(3.0).powf(F32(x));
168 let relative_error = calc_relative_error(exp_x, expected);
169
170 assert!(
171 relative_error <= MAX_ERROR,
172 "relative_error {} too large for input {} : {} vs {}",
173 relative_error,
174 x,
175 exp_x,
176 expected
177 );
178 }
179
180 for &(x, expected) in TEST_VECTORS_POW150 {
181 let exp_x = F32(150.0).powf(F32(x));
182 let relative_error = calc_relative_error(exp_x, expected);
183
184 assert!(
185 relative_error <= MAX_ERROR,
186 "relative_error {} too large for input {} : {} vs {}",
187 relative_error,
188 x,
189 exp_x,
190 expected
191 );
192 }
193
194 for &(base_input, power_input, expected) in TEST_VECTORS_MISC {
195 let exp_x = F32(base_input).powf(F32(power_input));
196 let relative_error = calc_relative_error(exp_x, expected);
197
198 assert!(
199 relative_error <= MAX_ERROR,
200 "relative_error {} too large for input {}.powf({}) : {} vs {}",
201 relative_error,
202 base_input,
203 power_input,
204 exp_x,
205 expected
206 );
207 }
208 }
209}
210