authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2017-06-17 20:39:45-04:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2017-06-17 20:39:45-04:00
log62323eeb75b94829f0941cf97067780f0999573f
treecc167f7fa556d5b618723f8f53d046233647c471
parent4c16f9a3c35b23b9917f2a27b91ba8cd20e6fd82

std: refactor pow to be generic

4 files changed, 36 insertions(+), 190 deletions(-)

std/math/exp.zig+1-1
...@@ -11,7 +11,7 @@ pub fn exp(x: var) -> @typeOf(x) {...@@ -11,7 +11,7 @@ pub fn exp(x: var) -> @typeOf(x) {
11}11}
1212
13fn exp32(x_: f32) -> f32 {13fn exp32(x_: f32) -> f32 {
14 const half = []const f32 { 0.5, -0.5 };14 const half = []f32 { 0.5, -0.5 };
15 const ln2hi = 6.9314575195e-1;15 const ln2hi = 6.9314575195e-1;
16 const ln2lo = 1.4286067653e-6;16 const ln2lo = 1.4286067653e-6;
17 const invln2 = 1.4426950216e+0;17 const invln2 = 1.4426950216e+0;
std/math/ln.zig+3-3
...@@ -120,12 +120,12 @@ fn lnd(x_: f64) -> f64 {...@@ -120,12 +120,12 @@ fn lnd(x_: f64) -> f64 {
120 s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi120 s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi
121}121}
122122
123test "log" {123test "math.ln" {
124 assert(ln(f32(0.2)) == lnf(0.2));124 assert(ln(f32(0.2)) == lnf(0.2));
125 assert(ln(f64(0.2)) == lnd(0.2));125 assert(ln(f64(0.2)) == lnd(0.2));
126}126}
127127
128test "logf" {128test "math.ln32" {
129 const epsilon = 0.000001;129 const epsilon = 0.000001;
130130
131 assert(math.approxEq(f32, lnf(0.2), -1.609438, epsilon));131 assert(math.approxEq(f32, lnf(0.2), -1.609438, epsilon));
...@@ -136,7 +136,7 @@ test "logf" {...@@ -136,7 +136,7 @@ test "logf" {
136 assert(math.approxEq(f32, lnf(123123.234375), 11.720941, epsilon));136 assert(math.approxEq(f32, lnf(123123.234375), 11.720941, epsilon));
137}137}
138138
139test "logd" {139test "math.ln64" {
140 const epsilon = 0.000001;140 const epsilon = 0.000001;
141141
142 assert(math.approxEq(f64, lnd(0.2), -1.609438, epsilon));142 assert(math.approxEq(f64, lnd(0.2), -1.609438, epsilon));
std/math/pow.zig+28-177
...@@ -1,21 +1,12 @@...@@ -1,21 +1,12 @@
1const math = @import("index.zig");1const math = @import("index.zig");
2const assert = @import("../debug.zig").assert;2const assert = @import("../debug.zig").assert;
33
4// This implementation is taken from the go stlib, musl is a bit more complex.
4pub fn pow(comptime T: type, x: T, y: T) -> T {5pub fn pow(comptime T: type, x: T, y: T) -> T {
5 switch (T) {6 if (T != f32 and T != f64) {
6 f32 => @inlineCall(pow32, x, y),7 @compileError("pow not implemented for " ++ @typeName(T));
7 f64 => @inlineCall(pow64, x, y),
8 else => @compileError("pow not implemented for " ++ @typeName(T)),
9 }8 }
10}
119
12fn isOddInteger(x: f64) -> bool {
13 const r = math.modf(x);
14 r.fpart == 0.0 and i64(r.ipart) & 1 == 1
15}
16
17// This implementation is taken from the go stlib, musl is a bit more complex.
18fn pow32(x: f32, y: f32) -> f32 {
19 // pow(x, +-0) = 1 for all x10 // pow(x, +-0) = 1 for all x
20 // pow(1, y) = 1 for all y11 // pow(1, y) = 1 for all y
21 if (y == 0 or x == 1) {12 if (y == 0 or x == 1) {
...@@ -25,7 +16,7 @@ fn pow32(x: f32, y: f32) -> f32 {...@@ -25,7 +16,7 @@ fn pow32(x: f32, y: f32) -> f32 {
25 // pow(nan, y) = nan for all y16 // pow(nan, y) = nan for all y
26 // pow(x, nan) = nan for all x17 // pow(x, nan) = nan for all x
27 if (math.isNan(x) or math.isNan(y)) {18 if (math.isNan(x) or math.isNan(y)) {
28 return math.nan(f32);19 return math.nan(T);
29 }20 }
3021
31 // pow(x, 1) = x for all x22 // pow(x, 1) = x for all x
...@@ -46,11 +37,11 @@ fn pow32(x: f32, y: f32) -> f32 {...@@ -46,11 +37,11 @@ fn pow32(x: f32, y: f32) -> f32 {
46 if (y < 0) {37 if (y < 0) {
47 // pow(+-0, y) = +- 0 for y an odd integer38 // pow(+-0, y) = +- 0 for y an odd integer
48 if (isOddInteger(y)) {39 if (isOddInteger(y)) {
49 return math.copysign(f32, math.inf(f32), x);40 return math.copysign(T, math.inf(T), x);
50 }41 }
51 // pow(+-0, y) = +inf for y an even integer42 // pow(+-0, y) = +inf for y an even integer
52 else {43 else {
53 return math.inf(f32);44 return math.inf(T);
54 }45 }
55 } else {46 } else {
56 if (isOddInteger(y)) {47 if (isOddInteger(y)) {
...@@ -74,13 +65,13 @@ fn pow32(x: f32, y: f32) -> f32 {...@@ -74,13 +65,13 @@ fn pow32(x: f32, y: f32) -> f32 {
74 // pow(x, -inf) = +inf for |x| < 165 // pow(x, -inf) = +inf for |x| < 1
75 // pow(x, +inf) = +inf for |x| > 166 // pow(x, +inf) = +inf for |x| > 1
76 else {67 else {
77 return math.inf(f32);68 return math.inf(T);
78 }69 }
79 }70 }
8071
81 if (math.isInf(x)) {72 if (math.isInf(x)) {
82 if (math.isNegativeInf(x)) {73 if (math.isNegativeInf(x)) {
83 return pow32(1 / x, -y);74 return pow(T, 1 / x, -y);
84 }75 }
85 // pow(+inf, y) = +0 for y < 076 // pow(+inf, y) = +0 for y < 0
86 else if (y < 0) {77 else if (y < 0) {
...@@ -88,7 +79,7 @@ fn pow32(x: f32, y: f32) -> f32 {...@@ -88,7 +79,7 @@ fn pow32(x: f32, y: f32) -> f32 {
88 }79 }
89 // pow(+inf, y) = +0 for y > 080 // pow(+inf, y) = +0 for y > 0
90 else if (y > 0) {81 else if (y > 0) {
91 return math.inf(f32);82 return math.inf(T);
92 }83 }
93 }84 }
9485
...@@ -104,14 +95,14 @@ fn pow32(x: f32, y: f32) -> f32 {...@@ -104,14 +95,14 @@ fn pow32(x: f32, y: f32) -> f32 {
104 var yf = r1.fpart;95 var yf = r1.fpart;
10596
106 if (yf != 0 and x < 0) {97 if (yf != 0 and x < 0) {
107 return math.nan(f32);98 return math.nan(T);
108 }99 }
109 if (yi >= 1 << 31) {100 if (yi >= 1 << (T.bit_count - 1)) {
110 return math.exp(y * math.ln(x));101 return math.exp(y * math.ln(x));
111 }102 }
112103
113 // a = a1 * 2^ae104 // a = a1 * 2^ae
114 var a1: f32 = 1.0;105 var a1: T = 1.0;
115 var ae: i32 = 0;106 var ae: i32 = 0;
116107
117 // a *= x^yf108 // a *= x^yf
...@@ -151,166 +142,26 @@ fn pow32(x: f32, y: f32) -> f32 {...@@ -151,166 +142,26 @@ fn pow32(x: f32, y: f32) -> f32 {
151 math.scalbn(a1, ae)142 math.scalbn(a1, ae)
152}143}
153144
154// This implementation is taken from the go stlib, musl is a bit more complex.145fn isOddInteger(x: f64) -> bool {
155fn pow64(x: f64, y: f64) -> f64 {146 const r = math.modf(x);
156 // pow(x, +-0) = 1 for all x147 r.fpart == 0.0 and i64(r.ipart) & 1 == 1
157 // pow(1, y) = 1 for all y
158 if (y == 0 or x == 1) {
159 return 1;
160 }
161
162 // pow(nan, y) = nan for all y
163 // pow(x, nan) = nan for all x
164 if (math.isNan(x) or math.isNan(y)) {
165 return math.nan(f64);
166 }
167
168 // pow(x, 1) = x for all x
169 if (y == 1) {
170 return x;
171 }
172
173 // special case sqrt
174 if (y == 0.5) {
175 return math.sqrt(x);
176 }
177
178 if (y == -0.5) {
179 return 1 / math.sqrt(x);
180 }
181
182 if (x == 0) {
183 if (y < 0) {
184 // pow(+-0, y) = +- 0 for y an odd integer
185 if (isOddInteger(y)) {
186 return math.copysign(f64, math.inf(f64), x);
187 }
188 // pow(+-0, y) = +inf for y an even integer
189 else {
190 return math.inf(f64);
191 }
192 } else {
193 if (isOddInteger(y)) {
194 return x;
195 } else {
196 return 0;
197 }
198 }
199 }
200
201 if (math.isInf(y)) {
202 // pow(-1, inf) = -1 for all x
203 if (x == -1) {
204 return -1;
205 }
206 // pow(x, +inf) = +0 for |x| < 1
207 // pow(x, -inf) = +0 for |x| > 1
208 else if ((math.fabs(x) < 1) == math.isInf(y)) {
209 return 0;
210 }
211 // pow(x, -inf) = +inf for |x| < 1
212 // pow(x, +inf) = +inf for |x| > 1
213 else {
214 return math.inf(f64);
215 }
216 }
217
218 if (math.isInf(x)) {
219 if (math.isInf(x)) {
220 return pow64(1 / x, -y);
221 }
222 // pow(+inf, y) = +0 for y < 0
223 else if (y < 0) {
224 return 0;
225 }
226 // pow(+inf, y) = +0 for y > 0
227 else if (y > 0) {
228 return math.inf(f64);
229 }
230 }
231
232 var ay = y;
233 var flip = false;
234 if (ay < 0) {
235 ay = -ay;
236 flip = true;
237 }
238
239 const r1 = math.modf(ay);
240 var yi = r1.ipart;
241 var yf = r1.fpart;
242
243 if (yf != 0 and x < 0) {
244 return math.nan(f64);
245 }
246 if (yi >= 1 << 63) {
247 return math.exp(y * math.ln(x));
248 }
249
250 // a = a1 * 2^ae
251 var a1: f64 = 1.0;
252 var ae: i32 = 0;
253
254 // a *= x^yf
255 if (yf != 0) {
256 if (yf > 0.5) {
257 yf -= 1;
258 yi += 1;
259 }
260 a1 = math.exp(yf * math.ln(x));
261 }
262
263 // a *= x^yi
264 const r2 = math.frexp(x);
265 var xe = r2.exponent;
266 var x1 = r2.significand;
267
268 var i = i64(yi);
269 while (i != 0) : (i >>= 1) {
270 if (i & 1 == 1) {
271 a1 *= x1;
272 ae += xe;
273 }
274 x1 *= x1;
275 xe <<= 1;
276 if (x1 < 0.5) {
277 x1 += x1;
278 xe -= 1;
279 }
280 }
281
282 // a *= a1 * 2^ae
283 if (flip) {
284 a1 = 1 / a1;
285 ae = -ae;
286 }
287
288 math.scalbn(a1, ae)
289}
290
291test "pow" {
292 assert(pow(f32, 0.2, 3.3) == pow32(0.2, 3.3));
293 assert(pow(f64, 0.2, 3.3) == pow64(0.2, 3.3));
294}148}
295149
296test "pow32" {150test "math.pow" {
297 const epsilon = 0.000001;151 const epsilon = 0.000001;
298152
299 // assert(math.approxEq(f32, pow32(0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero153 // assert(math.approxEq(f32, pow(f32, 0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero
300 assert(math.approxEq(f32, pow32(0.8923, 3.3), 0.686572, epsilon));154 assert(math.approxEq(f32, pow(f32, 0.8923, 3.3), 0.686572, epsilon));
301 assert(math.approxEq(f32, pow32(0.2, 3.3), 0.004936, epsilon));155 assert(math.approxEq(f32, pow(f32, 0.2, 3.3), 0.004936, epsilon));
302 assert(math.approxEq(f32, pow32(1.5, 3.3), 3.811546, epsilon));156 assert(math.approxEq(f32, pow(f32, 1.5, 3.3), 3.811546, epsilon));
303 assert(math.approxEq(f32, pow32(37.45, 3.3), 155736.703125, epsilon));157 assert(math.approxEq(f32, pow(f32, 37.45, 3.3), 155736.703125, epsilon));
304 assert(math.approxEq(f32, pow32(89.123, 3.3), 2722489.5, epsilon));158 assert(math.approxEq(f32, pow(f32, 89.123, 3.3), 2722489.5, epsilon));
305}
306159
307test "pow64" {
308 const epsilon = 0.000001;
309160
310 // assert(math.approxEq(f32, pow32(0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero161 // assert(math.approxEq(f32, pow(f64, 0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero
311 assert(math.approxEq(f64, pow64(0.8923, 3.3), 0.686572, epsilon));162 assert(math.approxEq(f64, pow(f64, 0.8923, 3.3), 0.686572, epsilon));
312 assert(math.approxEq(f64, pow64(0.2, 3.3), 0.004936, epsilon));163 assert(math.approxEq(f64, pow(f64, 0.2, 3.3), 0.004936, epsilon));
313 assert(math.approxEq(f64, pow64(1.5, 3.3), 3.811546, epsilon));164 assert(math.approxEq(f64, pow(f64, 1.5, 3.3), 3.811546, epsilon));
314 assert(math.approxEq(f64, pow64(37.45, 3.3), 155736.7160616, epsilon));165 assert(math.approxEq(f64, pow(f64, 37.45, 3.3), 155736.7160616, epsilon));
315 assert(math.approxEq(f64, pow64(89.123, 3.3), 2722490.231436, epsilon));166 assert(math.approxEq(f64, pow(f64, 89.123, 3.3), 2722490.231436, epsilon));
316}167}
std/special/builtin.zig+4-9
...@@ -60,7 +60,7 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {...@@ -60,7 +60,7 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {
60 if (ex == 0) {60 if (ex == 0) {
61 i = ux <<% exp_bits;61 i = ux <<% exp_bits;
62 while (i >> bits_minus_1 == 0) : ({ex -= 1; i <<%= 1}) {}62 while (i >> bits_minus_1 == 0) : ({ex -= 1; i <<%= 1}) {}
63 ux <<%= twosComplementCast(uint, -ex + 1);63 ux <<%= @bitCast(u32, -ex + 1);
64 } else {64 } else {
65 ux &= @maxValue(uint) >> exp_bits;65 ux &= @maxValue(uint) >> exp_bits;
66 ux |= 1 <<% digits;66 ux |= 1 <<% digits;
...@@ -68,7 +68,7 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {...@@ -68,7 +68,7 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {
68 if (ey == 0) {68 if (ey == 0) {
69 i = uy <<% exp_bits;69 i = uy <<% exp_bits;
70 while (i >> bits_minus_1 == 0) : ({ey -= 1; i <<%= 1}) {}70 while (i >> bits_minus_1 == 0) : ({ey -= 1; i <<%= 1}) {}
71 uy <<= twosComplementCast(uint, -ey + 1);71 uy <<= @bitCast(u32, -ey + 1);
72 } else {72 } else {
73 uy &= @maxValue(uint) >> exp_bits;73 uy &= @maxValue(uint) >> exp_bits;
74 uy |= 1 <<% digits;74 uy |= 1 <<% digits;
...@@ -95,9 +95,9 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {...@@ -95,9 +95,9 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {
95 // scale result up95 // scale result up
96 if (ex > 0) {96 if (ex > 0) {
97 ux -%= 1 <<% digits;97 ux -%= 1 <<% digits;
98 ux |= twosComplementCast(uint, ex) <<% digits;98 ux |= @bitCast(u32, ex) <<% digits;
99 } else {99 } else {
100 ux >>= twosComplementCast(uint, -ex + 1);100 ux >>= @bitCast(u32, -ex + 1);
101 }101 }
102 if (T == f32) {102 if (T == f32) {
103 ux |= sx;103 ux |= sx;
...@@ -116,8 +116,3 @@ fn isNan(comptime T: type, bits: T) -> bool {...@@ -116,8 +116,3 @@ fn isNan(comptime T: type, bits: T) -> bool {
116 unreachable;116 unreachable;
117 }117 }
118}118}
119
120// TODO this should be a builtin function and it shouldn't do a ptr cast
121fn twosComplementCast(comptime T: type, src: var) -> T {
122 return *@ptrCast(&const @IntType(T.is_signed, @typeOf(src).bit_count), &src);
123}