fix atan2 so unresolved sign is not treated as proven zero.

Author: timschmidt

src/computable/node.rs

@@ -24,6 +24,7 @@ use std::{
 };
 
 pub type Precision = i32;
+const ATAN2_SIGN_REFINEMENT_FLOOR: Precision = -4096;
 
 #[derive(Clone, Debug, PartialEq, Default)]
 pub(crate) enum Cache {
@@ -419,6 +420,14 @@ fn public_sign(sign: Sign) -> RealSign {
     }
 }
 
+fn private_sign(sign: RealSign) -> Sign {
+    match sign {
+        RealSign::Negative => Sign::Minus,
+        RealSign::Zero => Sign::NoSign,
+        RealSign::Positive => Sign::Plus,
+    }
+}
+
 use std::sync::Arc;
 use std::sync::atomic::AtomicBool;
 
@@ -3217,28 +3226,71 @@ impl Computable {
     /// - axes return exact constants: `pi/2`, `-pi/2`, `pi`, or zero.
     /// - the origin `(0, 0)` returns zero, matching `f64::atan2`.
     pub fn atan2(self, x: Computable) -> Computable {
-        let y_sign = self.sign();
-        let x_sign = x.sign();
+        let y_sign = self.structural_facts().sign.map(private_sign);
+        let x_sign = x.structural_facts().sign.map(private_sign);
         match (y_sign, x_sign) {
-            (Sign::NoSign, Sign::NoSign) | (Sign::NoSign, Sign::Plus) => {
+            (Some(Sign::NoSign), Some(Sign::NoSign)) | (Some(Sign::NoSign), Some(Sign::Plus)) => {
                 crate::trace_dispatch!("computable", "atan2", "axis-zero-y");
                 return Self::zero();
             }
-            (Sign::NoSign, Sign::Minus) => {
+            (Some(Sign::NoSign), Some(Sign::Minus)) => {
                 crate::trace_dispatch!("computable", "atan2", "axis-negative-x");
                 return Self::pi();
             }
-            (Sign::Plus, Sign::NoSign) => {
+            (Some(Sign::Plus), Some(Sign::NoSign)) => {
                 crate::trace_dispatch!("computable", "atan2", "axis-positive-y");
                 return Self::pi().shift_right(1);
             }
-            (Sign::Minus, Sign::NoSign) => {
+            (Some(Sign::Minus), Some(Sign::NoSign)) => {
                 crate::trace_dispatch!("computable", "atan2", "axis-negative-y");
                 return Self::pi().shift_right(1).negate();
             }
             _ => {}
         }
-        let base = self.multiply(x.clone().inverse()).atan();
+        match (
+            y_sign.unwrap_or_else(|| self.sign()),
+            x_sign.unwrap_or_else(|| x.sign()),
+        ) {
+            (Sign::NoSign, Sign::Plus) => {
+                crate::trace_dispatch!("computable", "atan2", "quadrant-right");
+                return self.multiply(x.inverse()).atan();
+            }
+            (Sign::NoSign, Sign::NoSign) => {
+                crate::trace_dispatch!("computable", "atan2", "unresolved-origin");
+                return Self::zero();
+            }
+            (Sign::NoSign, Sign::Minus) => {
+                let y_sign = self
+                    .sign_until(ATAN2_SIGN_REFINEMENT_FLOOR)
+                    .map(private_sign);
+                return match y_sign {
+                    Some(Sign::Minus) => {
+                        crate::trace_dispatch!("computable", "atan2", "quadrant-lower-left");
+                        self.multiply(x.inverse()).atan().add(Self::pi().negate())
+                    }
+                    Some(Sign::Plus) => {
+                        crate::trace_dispatch!("computable", "atan2", "quadrant-upper-left");
+                        self.multiply(x.inverse()).atan().add(Self::pi())
+                    }
+                    _ => {
+                        crate::trace_dispatch!("computable", "atan2", "axis-negative-x");
+                        Self::pi()
+                    }
+                };
+            }
+            (Sign::Plus, Sign::NoSign) => {
+                crate::trace_dispatch!("computable", "atan2", "half-angle-positive-y");
+                return Self::atan2_half_angle(self, x);
+            }
+            (Sign::Minus, Sign::NoSign) => {
+                crate::trace_dispatch!("computable", "atan2", "half-angle-negative-y");
+                return Self::atan2_half_angle(self, x);
+            }
+            _ => {}
+        }
+        let x_sign = x_sign.unwrap_or_else(|| x.sign());
+        let y_sign = y_sign.unwrap_or_else(|| self.sign());
+        let base = self.multiply(x.inverse()).atan();
         if x_sign == Sign::Plus {
             crate::trace_dispatch!("computable", "atan2", "quadrant-right");
             base
@@ -3251,6 +3303,11 @@ impl Computable {
         }
     }
 
+    fn atan2_half_angle(y: Computable, x: Computable) -> Computable {
+        let radius = x.clone().square().add(y.clone().square()).sqrt();
+        y.multiply(radius.add(x).inverse()).atan().shift_left(1)
+    }
+
     /// Inverse sine of this number.
     pub fn asin(self) -> Computable {
         if let Some(rational) = self.exact_rational() {

src/real/arithmetic.rs

@@ -4032,37 +4032,37 @@ impl Real {
         // more expensive than reading already-derivable structural facts. Only
         // descend to the refinement path when structural inspection cannot
         // decide for one of the inputs.
-        let y_sign = match self.structural_facts().sign {
-            Some(RealSign::Zero) => Sign::NoSign,
-            Some(RealSign::Positive) => Sign::Plus,
-            Some(RealSign::Negative) => Sign::Minus,
-            None => self.best_sign(),
-        };
-        let x_sign = match x.structural_facts().sign {
-            Some(RealSign::Zero) => Sign::NoSign,
-            Some(RealSign::Positive) => Sign::Plus,
-            Some(RealSign::Negative) => Sign::Minus,
-            None => x.best_sign(),
-        };
+        let y_sign = self.structural_facts().sign.map(num_sign_from_real);
+        let x_sign = x.structural_facts().sign.map(num_sign_from_real);
         match (y_sign, x_sign) {
-            (Sign::NoSign, Sign::NoSign) | (Sign::NoSign, Sign::Plus) => {
+            (Some(Sign::NoSign), Some(Sign::NoSign)) | (Some(Sign::NoSign), Some(Sign::Plus)) => {
                 crate::trace_dispatch!("real", "atan2", "axis-zero-y");
                 return Self::zero();
             }
-            (Sign::NoSign, Sign::Minus) => {
+            (Some(Sign::NoSign), Some(Sign::Minus)) => {
                 crate::trace_dispatch!("real", "atan2", "axis-negative-x");
                 return Self::pi();
             }
-            (Sign::Plus, Sign::NoSign) => {
+            (Some(Sign::Plus), Some(Sign::NoSign)) => {
                 crate::trace_dispatch!("real", "atan2", "axis-positive-y");
                 return Self::pi_fraction(1, 2);
             }
-            (Sign::Minus, Sign::NoSign) => {
+            (Some(Sign::Minus), Some(Sign::NoSign)) => {
                 crate::trace_dispatch!("real", "atan2", "axis-negative-y");
                 return Self::pi_fraction(-1, 2);
             }
             _ => {}
         }
+        let y_sign = y_sign.unwrap_or_else(|| self.best_sign());
+        let x_sign = x_sign.unwrap_or_else(|| x.best_sign());
+        if y_sign == Sign::NoSign && x_sign != Sign::Plus {
+            crate::trace_dispatch!("real", "atan2", "generic-computable");
+            return Self::irrational_from_computable(self.fold().atan2(x.fold()));
+        }
+        if x_sign == Sign::NoSign {
+            crate::trace_dispatch!("real", "atan2", "generic-computable");
+            return Self::irrational_from_computable(self.fold().atan2(x.fold()));
+        }
         let ratio = (self / &x).expect("nonzero x rules out divide-by-zero");
         let base = ratio.atan().expect("Real::atan is total");
         if x_sign == Sign::Plus {
@@ -4140,7 +4140,7 @@ impl Real {
         crate::trace_dispatch!("real", "tanh", "generic-exp-identity");
         let positive = self.clone().exp()?;
         let negative = self.neg().exp()?;
-        (positive.clone() - negative.clone()) / (positive + negative)
+        (&positive - &negative) / (positive + negative)
     }
 
     /// The inverse hyperbolic sine of this Real.
@@ -4711,6 +4711,14 @@ fn real_sign_from_num(sign: Sign) -> RealSign {
     }
 }
 
+fn num_sign_from_real(sign: RealSign) -> Sign {
+    match sign {
+        RealSign::Negative => Sign::Minus,
+        RealSign::Zero => Sign::NoSign,
+        RealSign::Positive => Sign::Plus,
+    }
+}
+
 fn multiply_public_sign(left: Option<RealSign>, right: Option<RealSign>) -> Option<RealSign> {
     match (left?, right?) {
         (RealSign::Zero, _) | (_, RealSign::Zero) => Some(RealSign::Zero),

src/real/tests.rs

@@ -1913,6 +1913,24 @@ mod tests {
         }
     }
 
+    #[test]
+    fn atan2_unresolved_positive_y_does_not_collapse_to_axis() {
+        let tiny = Real::new(
+            Rational::from_bigint_fraction(
+                num::BigInt::from(1_u8),
+                num::BigUint::from(1_u8) << 2500,
+            )
+            .unwrap(),
+        );
+        let y = (Real::pi() + tiny.clone()) - Real::pi();
+        assert_eq!(y.structural_facts().sign, None);
+
+        let got = y.atan2(Real::one());
+        let expected = tiny.atan2(Real::one());
+        assert_ne!(got, Real::zero());
+        assert_eq!(got.to_f64_approx(), expected.to_f64_approx());
+    }
+
     #[test]
     fn atan2_is_consistent_under_uniform_positive_scaling() {
         // atan2(ky, kx) = atan2(y, x) for k > 0. Pick coords whose |y/x|
@@ -1990,6 +2008,51 @@ mod tests {
         }
     }
 
+    #[test]
+    fn computable_atan2_unresolved_positive_y_does_not_collapse_to_axis() {
+        use crate::Computable;
+        use num::Zero;
+
+        let tiny = Rational::from_bigint_fraction(
+            num::BigInt::from(1_u8),
+            num::BigUint::from(1_u8) << 2500,
+        )
+        .unwrap();
+        let y = Computable::pi()
+            .add(Computable::rational(tiny.clone()))
+            .add(Computable::pi().negate());
+        assert_eq!(y.sign(), num::bigint::Sign::NoSign);
+
+        let got = y.atan2(Computable::one()).approx(-2600);
+        let expected = Computable::rational(tiny)
+            .atan2(Computable::one())
+            .approx(-2600);
+        assert!(!got.is_zero());
+        assert_eq!(got, expected);
+    }
+
+    #[test]
+    fn computable_atan2_unresolved_negative_y_on_negative_x_keeps_lower_branch() {
+        use crate::Computable;
+
+        let tiny = Rational::from_bigint_fraction(
+            num::BigInt::from(1_u8),
+            num::BigUint::from(1_u8) << 2500,
+        )
+        .unwrap();
+        let y = Computable::pi()
+            .add(Computable::rational(tiny.clone()).negate())
+            .add(Computable::pi().negate());
+        assert_eq!(y.sign(), num::bigint::Sign::NoSign);
+
+        let got = y.atan2(Computable::one().negate()).approx(-2600);
+        let expected = Computable::rational(tiny)
+            .negate()
+            .atan2(Computable::one().negate())
+            .approx(-2600);
+        assert_eq!(got, expected);
+    }
+
     #[test]
     fn dot2_refs_all_zero_returns_zero() {
         let left = [Real::zero(), Real::zero()];