Fix recurrence order in dawson's exponential-sum branch

src/dawson.rs

@@ -40,11 +40,13 @@ pub fn dawson(x: f64) -> f64 {
 
         for i in 0..NMAX {
             let cval = f64::exp(-((2. * i as f64 + 1.) * h).powi(2));
+            sum += cval * (e1 / d1 + 1.0 / (d2 * e1));
+
+            // Advance the recurrence after accumulating, as in the
+            // Numerical Recipes for-loop post-clause.
             d1 += 2.;
             d2 -= 2.;
             e1 *= e2;
-
-            sum += cval * (e1 / d1 + 1.0 / (d2 * e1));
         }
 
         ans = (1. / PI.sqrt()) * sign(f64::exp(-xp * xp), x) * sum;

tests/dawson_test.rs

@@ -1,9 +1,76 @@
 use puruspe::dawson;
 
+// Accuracy tolerance. `dawson` documents ~2e-7; the differential grid below
+// agrees with mpmath to better than 6e-8.
+const TOL: f64 = 1e-7;
+
 #[test]
 fn dawson_test() {
     let expected = 0.13818492867352312; // From Wolfram|Alpha
     let res = dawson(0.14);
     dbg!(expected, res);
-    assert!((expected - res).abs() < 1e-7); // Verify accuracy.
+    assert!((expected - res).abs() < TOL); // Verify accuracy.
+}
+
+// Series branch (|x| < 0.2): the Maclaurin approximation must stay correct.
+// Oracle: mpmath D(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x) at dps 40.
+#[test]
+fn dawson_series_branch_test() {
+    let cases = [
+        (0.0, 0.0),
+        (0.05, 0.049916749940476193),
+        (0.1, 0.099335992397852861),
+        (0.14, 0.13818492867352313),
+        (0.19, 0.18549268702269875),
+        (0.199, 0.19382855517484729),
+        (-0.1, -0.099335992397852861),
+    ];
+    for (x, expected) in cases {
+        let res = dawson(x);
+        println!("dawson({x}) = {res:.17e}  expected {expected:.17e}");
+        assert!((res - expected).abs() < TOL, "series branch wrong at x={}", x);
+    }
+}
+
+// Exponential-sum branch (|x| >= 0.2), which was previously wrong for every
+// argument (wrong sign near x=0.5, orders of magnitude off at x=10). Covers the
+// branch boundary, the sign-flip region, the asymptotic tail, and negatives.
+#[test]
+fn dawson_exp_sum_branch_test() {
+    let cases = [
+        (0.2, 0.19475103336802805),
+        (0.201, 0.19567275321979942),
+        (0.25, 0.23983916356289821),
+        (0.3, 0.28263166502131193),
+        (0.5, 0.4244363835020223),
+        (0.7, 0.51050405755923178),
+        (1.0, 0.53807950691276842),
+        (1.5, 0.42824907108539863),
+        (2.0, 0.30134038892379197),
+        (3.0, 0.17827103061055829),
+        (5.0, 0.10213407442427684),
+        (8.0, 0.063000198707553388),
+        (10.0, 0.050253847187598528),
+        (-0.3, -0.28263166502131193),
+        (-0.5, -0.4244363835020223),
+        (-1.0, -0.53807950691276842),
+        (-2.0, -0.30134038892379197),
+        (-5.0, -0.10213407442427684),
+    ];
+    for (x, expected) in cases {
+        let res = dawson(x);
+        println!("dawson({x}) = {res:.17e}  expected {expected:.17e}");
+        assert!((res - expected).abs() < TOL, "exp-sum branch wrong at x={}", x);
+    }
+}
+
+// Dawson's integral is odd: D(-x) = -D(x). Exercises both branches together.
+#[test]
+fn dawson_odd_symmetry_test() {
+    for x in [0.05, 0.14, 0.2, 0.5, 1.0, 3.0, 7.5] {
+        let a = dawson(x);
+        let b = dawson(-x);
+        println!("dawson({x}) = {a:.17e}, dawson({}) = {b:.17e}", -x);
+        assert!((a + b).abs() < 1e-15, "not odd at x={}", x);
+    }
 }

tests/faddeeva_test.rs

@@ -1043,3 +1043,31 @@ const W_OF_Z_TABLE: [(f64, f64); 1000] = [
     (0.02830791446808592, 0.02816645639318190),
     (0.02827946745423294, 0.02813843327633737),
 ];
+
+// Real-axis case: for y = 0 (and y just inside the |y| < 1e-8*|x| branch), the
+// imaginary part is Im w(x,0) = 2/sqrt(pi) * D(x), so it delegates to `dawson`.
+// The C-table test above only samples Re(z) == Im(z), never this branch.
+// Oracle: mpmath 2/sqrt(pi) * D(x) at dps 40.
+#[test]
+fn faddeeva_real_axis_im_test() {
+    let cases = [
+        (0.2, 0.21975300882280587),
+        (0.3, 0.31891568277156587),
+        (0.5, 0.47892517290104347),
+        (1.0, 0.60715770584139373),
+        (2.0, 0.3400262170660662),
+        (5.0, 0.11524596183093659),
+    ];
+    for (x, im_expected) in cases {
+        // Exactly on the real axis: Re w(x,0) = exp(-x^2).
+        let (re, im) = faddeeva(x, 0.0);
+        let re_expected = (-x * x as f64).exp();
+        println!("faddeeva({x},0) = ({re:.17e}, {im:.17e})  im_expected {im_expected:.17e}");
+        assert!((re - re_expected).abs() < 1e-15, "re wrong at x={}", x);
+        assert!((im - im_expected).abs() < 1e-7, "im wrong at x={}", x);
+
+        // Just inside the near-real-axis branch: y = 1e-12 * x.
+        let (_, im2) = faddeeva(x, 1e-12 * x);
+        assert!((im2 - im_expected).abs() < 1e-7, "near-axis im wrong at x={}", x);
+    }
+}