ff14451b4a
Comparisons with absolute epsilons are usually useful when comparing numbers to zero, for non-zero numbers it's advised to switch to relative epsilons instead to obtain meaningful results (check [1] for more details). The new API introduces approxEqAbs and approxEqRel, where the former aliases and deprecated the old `approxEq`, allowing the user to pick the right tool for the job. The documentation is meant to guide the user in the choice of the correct alternative. [1] https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
144 lines
3.4 KiB
Zig
144 lines
3.4 KiB
Zig
// SPDX-License-Identifier: MIT
|
|
// Copyright (c) 2015-2020 Zig Contributors
|
|
// This file is part of [zig](https://ziglang.org/), which is MIT licensed.
|
|
// The MIT license requires this copyright notice to be included in all copies
|
|
// and substantial portions of the software.
|
|
// Ported from musl, which is licensed under the MIT license:
|
|
// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
|
|
//
|
|
// https://git.musl-libc.org/cgit/musl/tree/src/complex/catanf.c
|
|
// https://git.musl-libc.org/cgit/musl/tree/src/complex/catan.c
|
|
|
|
const std = @import("../../std.zig");
|
|
const builtin = @import("builtin");
|
|
const testing = std.testing;
|
|
const math = std.math;
|
|
const cmath = math.complex;
|
|
const Complex = cmath.Complex;
|
|
|
|
/// Returns the arc-tangent of z.
|
|
pub fn atan(z: anytype) @TypeOf(z) {
|
|
const T = @TypeOf(z.re);
|
|
return switch (T) {
|
|
f32 => atan32(z),
|
|
f64 => atan64(z),
|
|
else => @compileError("atan not implemented for " ++ @typeName(z)),
|
|
};
|
|
}
|
|
|
|
fn redupif32(x: f32) f32 {
|
|
const DP1 = 3.140625;
|
|
const DP2 = 9.67502593994140625e-4;
|
|
const DP3 = 1.509957990978376432e-7;
|
|
|
|
var t = x / math.pi;
|
|
if (t >= 0.0) {
|
|
t += 0.5;
|
|
} else {
|
|
t -= 0.5;
|
|
}
|
|
|
|
const u = @intToFloat(f32, @floatToInt(i32, t));
|
|
return ((x - u * DP1) - u * DP2) - t * DP3;
|
|
}
|
|
|
|
fn atan32(z: Complex(f32)) Complex(f32) {
|
|
const maxnum = 1.0e38;
|
|
|
|
const x = z.re;
|
|
const y = z.im;
|
|
|
|
if ((x == 0.0) and (y > 1.0)) {
|
|
// overflow
|
|
return Complex(f32).new(maxnum, maxnum);
|
|
}
|
|
|
|
const x2 = x * x;
|
|
var a = 1.0 - x2 - (y * y);
|
|
if (a == 0.0) {
|
|
// overflow
|
|
return Complex(f32).new(maxnum, maxnum);
|
|
}
|
|
|
|
var t = 0.5 * math.atan2(f32, 2.0 * x, a);
|
|
var w = redupif32(t);
|
|
|
|
t = y - 1.0;
|
|
a = x2 + t * t;
|
|
if (a == 0.0) {
|
|
// overflow
|
|
return Complex(f32).new(maxnum, maxnum);
|
|
}
|
|
|
|
t = y + 1.0;
|
|
a = (x2 + (t * t)) / a;
|
|
return Complex(f32).new(w, 0.25 * math.ln(a));
|
|
}
|
|
|
|
fn redupif64(x: f64) f64 {
|
|
const DP1 = 3.14159265160560607910;
|
|
const DP2 = 1.98418714791870343106e-9;
|
|
const DP3 = 1.14423774522196636802e-17;
|
|
|
|
var t = x / math.pi;
|
|
if (t >= 0.0) {
|
|
t += 0.5;
|
|
} else {
|
|
t -= 0.5;
|
|
}
|
|
|
|
const u = @intToFloat(f64, @floatToInt(i64, t));
|
|
return ((x - u * DP1) - u * DP2) - t * DP3;
|
|
}
|
|
|
|
fn atan64(z: Complex(f64)) Complex(f64) {
|
|
const maxnum = 1.0e308;
|
|
|
|
const x = z.re;
|
|
const y = z.im;
|
|
|
|
if ((x == 0.0) and (y > 1.0)) {
|
|
// overflow
|
|
return Complex(f64).new(maxnum, maxnum);
|
|
}
|
|
|
|
const x2 = x * x;
|
|
var a = 1.0 - x2 - (y * y);
|
|
if (a == 0.0) {
|
|
// overflow
|
|
return Complex(f64).new(maxnum, maxnum);
|
|
}
|
|
|
|
var t = 0.5 * math.atan2(f64, 2.0 * x, a);
|
|
var w = redupif64(t);
|
|
|
|
t = y - 1.0;
|
|
a = x2 + t * t;
|
|
if (a == 0.0) {
|
|
// overflow
|
|
return Complex(f64).new(maxnum, maxnum);
|
|
}
|
|
|
|
t = y + 1.0;
|
|
a = (x2 + (t * t)) / a;
|
|
return Complex(f64).new(w, 0.25 * math.ln(a));
|
|
}
|
|
|
|
const epsilon = 0.0001;
|
|
|
|
test "complex.catan32" {
|
|
const a = Complex(f32).new(5, 3);
|
|
const c = atan(a);
|
|
|
|
testing.expect(math.approxEqAbs(f32, c.re, 1.423679, epsilon));
|
|
testing.expect(math.approxEqAbs(f32, c.im, 0.086569, epsilon));
|
|
}
|
|
|
|
test "complex.catan64" {
|
|
const a = Complex(f64).new(5, 3);
|
|
const c = atan(a);
|
|
|
|
testing.expect(math.approxEqAbs(f64, c.re, 1.423679, epsilon));
|
|
testing.expect(math.approxEqAbs(f64, c.im, 0.086569, epsilon));
|
|
}
|