diff options
Diffstat (limited to 'library/compiler-builtins/libm/src/math/atan2.rs')
| -rw-r--r-- | library/compiler-builtins/libm/src/math/atan2.rs | 131 |
1 files changed, 131 insertions, 0 deletions
diff --git a/library/compiler-builtins/libm/src/math/atan2.rs b/library/compiler-builtins/libm/src/math/atan2.rs new file mode 100644 index 00000000000..c668731cf37 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/atan2.rs @@ -0,0 +1,131 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2.c */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ +/* atan2(y,x) + * Method : + * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). + * 2. Reduce x to positive by (if x and y are unexceptional): + * ARG (x+iy) = arctan(y/x) ... if x > 0, + * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, + * + * Special cases: + * + * ATAN2((anything), NaN ) is NaN; + * ATAN2(NAN , (anything) ) is NaN; + * ATAN2(+-0, +(anything but NaN)) is +-0 ; + * ATAN2(+-0, -(anything but NaN)) is +-pi ; + * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; + * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; + * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; + * ATAN2(+-INF,+INF ) is +-pi/4 ; + * ATAN2(+-INF,-INF ) is +-3pi/4; + * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +use super::{atan, fabs}; + +const PI: f64 = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */ +const PI_LO: f64 = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ + +/// Arctangent of y/x (f64) +/// +/// Computes the inverse tangent (arc tangent) of `y/x`. +/// Produces the correct result even for angles near pi/2 or -pi/2 (that is, when `x` is near 0). +/// Returns a value in radians, in the range of -pi to pi. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] +pub fn atan2(y: f64, x: f64) -> f64 { + if x.is_nan() || y.is_nan() { + return x + y; + } + let mut ix = (x.to_bits() >> 32) as u32; + let lx = x.to_bits() as u32; + let mut iy = (y.to_bits() >> 32) as u32; + let ly = y.to_bits() as u32; + if ((ix.wrapping_sub(0x3ff00000)) | lx) == 0 { + /* x = 1.0 */ + return atan(y); + } + let m = ((iy >> 31) & 1) | ((ix >> 30) & 2); /* 2*sign(x)+sign(y) */ + ix &= 0x7fffffff; + iy &= 0x7fffffff; + + /* when y = 0 */ + if (iy | ly) == 0 { + return match m { + 0 | 1 => y, /* atan(+-0,+anything)=+-0 */ + 2 => PI, /* atan(+0,-anything) = PI */ + _ => -PI, /* atan(-0,-anything) =-PI */ + }; + } + /* when x = 0 */ + if (ix | lx) == 0 { + return if m & 1 != 0 { -PI / 2.0 } else { PI / 2.0 }; + } + /* when x is INF */ + if ix == 0x7ff00000 { + if iy == 0x7ff00000 { + return match m { + 0 => PI / 4.0, /* atan(+INF,+INF) */ + 1 => -PI / 4.0, /* atan(-INF,+INF) */ + 2 => 3.0 * PI / 4.0, /* atan(+INF,-INF) */ + _ => -3.0 * PI / 4.0, /* atan(-INF,-INF) */ + }; + } else { + return match m { + 0 => 0.0, /* atan(+...,+INF) */ + 1 => -0.0, /* atan(-...,+INF) */ + 2 => PI, /* atan(+...,-INF) */ + _ => -PI, /* atan(-...,-INF) */ + }; + } + } + /* |y/x| > 0x1p64 */ + if ix.wrapping_add(64 << 20) < iy || iy == 0x7ff00000 { + return if m & 1 != 0 { -PI / 2.0 } else { PI / 2.0 }; + } + + /* z = atan(|y/x|) without spurious underflow */ + let z = if (m & 2 != 0) && iy.wrapping_add(64 << 20) < ix { + /* |y/x| < 0x1p-64, x<0 */ + 0.0 + } else { + atan(fabs(y / x)) + }; + match m { + 0 => z, /* atan(+,+) */ + 1 => -z, /* atan(-,+) */ + 2 => PI - (z - PI_LO), /* atan(+,-) */ + _ => (z - PI_LO) - PI, /* atan(-,-) */ + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + #[cfg_attr(x86_no_sse, ignore = "FIXME(i586): possible incorrect rounding")] + fn sanity_check() { + assert_eq!(atan2(0.0, 1.0), 0.0); + assert_eq!(atan2(0.0, -1.0), PI); + assert_eq!(atan2(-0.0, -1.0), -PI); + assert_eq!(atan2(3.0, 2.0), atan(3.0 / 2.0)); + assert_eq!(atan2(2.0, -1.0), atan(2.0 / -1.0) + PI); + assert_eq!(atan2(-2.0, -1.0), atan(-2.0 / -1.0) - PI); + } +} |