fix: handle aliased operands in the in-place operators

include/hyperjet/hyperjet.h

@@ -951,6 +951,11 @@ class DDScalar {
   Type &operator*=(const Type &b) {
     check_equal_size(size(), b.size());
 
+    // aliased operands would be read again after being overwritten below
+    if (this == &b) {
+      return *this = *this * b;
+    }
+
     const Data a_m_data = m_data;
 
     const Scalar da = b.f();
@@ -1029,6 +1034,11 @@ class DDScalar {
 
     check_equal_size(size(), b.size());
 
+    // aliased operands would be read again after being overwritten below
+    if (this == &b) {
+      return *this = *this / b;
+    }
+
     const Data a_m_data = m_data;
 
     const auto f = a_m_data[0] / b.f();
@@ -1947,6 +1957,11 @@ template <typename TScalar> class SScalar {
   friend Type operator*(const Scalar a, const Type &b) { return b * a; }
 
   Type &operator*=(const Type &b) {
+    // aliased operands would be read again after being overwritten below
+    if (this == &b) {
+      return *this = *this * b;
+    }
+
     const Scalar da = b.m_f;
     const Scalar db = m_f;
 
@@ -1998,6 +2013,11 @@ template <typename TScalar> class SScalar {
   Type &operator/=(const Type &b) {
     using std::pow;
 
+    // aliased operands would be read again after being overwritten below
+    if (this == &b) {
+      return *this = *this / b;
+    }
+
     const auto da = 1 / b.m_f;
     const auto db = -m_f / pow(b.m_f, 2);
 

python/tests/test_DDScalar.py

@@ -711,6 +711,25 @@ def test_pad_left(ctx):
             u.pad_left(new_size=5)
 
 
+@pytest.mark.parametrize("ctx", **test_data)
+def test_self_multiplication(ctx):
+    r = ctx.u1
+    r *= r
+
+    assert_allclose(r.data, (ctx.u1 * ctx.u1).data, atol=1e-16)
+
+
+@pytest.mark.parametrize("ctx", **test_data)
+def test_self_division(ctx):
+    r = ctx.u1
+    r /= r
+
+    expected = np.zeros(len(r.data))
+    expected[0] = 1
+
+    assert_allclose(r.data, expected, atol=1e-15)
+
+
 # serialization
 
 

python/tests/test_SScalar.py

@@ -40,3 +40,28 @@ def test_mul():
     assert_allclose(r.d("x"), 1)
     assert_allclose(r.d("y"), 2)
     assert_allclose(r.d("z"), 0)
+
+
+def test_self_multiplication():
+    r = hj.SScalar(f=3, d={"x": 1, "y": 6})
+    r *= r
+
+    u = hj.SScalar(f=3, d={"x": 1, "y": 6})
+    e = u * hj.SScalar(f=3, d={"x": 1, "y": 6})
+
+    assert_allclose(r.f, e.f)
+    assert_allclose(r.d("x"), e.d("x"))
+    assert_allclose(r.d("y"), e.d("y"))
+
+    assert_allclose(r.f, 9)
+    assert_allclose(r.d("x"), 6)
+    assert_allclose(r.d("y"), 36)
+
+
+def test_self_division():
+    r = hj.SScalar(f=3, d={"x": 1, "y": 6})
+    r /= r
+
+    assert_allclose(r.f, 1)
+    assert_allclose(r.d("x"), 0, atol=1e-15)
+    assert_allclose(r.d("y"), 0, atol=1e-15)

test/src/test.cpp

@@ -287,6 +287,49 @@ TEST_CASE("IncDiv") {
   REQUIRE(r1.h(2, 2) == doctest::Approx(1.43750));
 }
 
+// An in-place operator has to give the same result as its binary counterpart,
+// including when both operands are the same object.
+
+TEST_CASE("Self multiplication") {
+  auto r1 = dd1;
+  r1 *= r1;
+
+  const auto r2 = dd1 * dd1;
+
+  REQUIRE(r1.f() == doctest::Approx(r2.f()));
+  REQUIRE(r1.g(0) == doctest::Approx(r2.g(0)));
+  REQUIRE(r1.g(1) == doctest::Approx(r2.g(1)));
+  REQUIRE(r1.g(2) == doctest::Approx(r2.g(2)));
+  REQUIRE(r1.h(0, 0) == doctest::Approx(r2.h(0, 0)));
+  REQUIRE(r1.h(0, 1) == doctest::Approx(r2.h(0, 1)));
+  REQUIRE(r1.h(0, 2) == doctest::Approx(r2.h(0, 2)));
+  REQUIRE(r1.h(1, 1) == doctest::Approx(r2.h(1, 1)));
+  REQUIRE(r1.h(1, 2) == doctest::Approx(r2.h(1, 2)));
+  REQUIRE(r1.h(2, 2) == doctest::Approx(r2.h(2, 2)));
+
+  // dd1 is 3 with g0 = 1 and h00 = 0, so x^2 gives 9, 2*x*g0 and 2*g0^2
+  REQUIRE(r1.f() == doctest::Approx(9.0));
+  REQUIRE(r1.g(0) == doctest::Approx(6.0));
+  REQUIRE(r1.h(0, 0) == doctest::Approx(2.0));
+}
+
+TEST_CASE("Self division") {
+  auto r1 = dd1;
+  r1 /= r1;
+
+  // x / x is 1 and all of its derivatives vanish
+  REQUIRE(r1.f() == doctest::Approx(1.0));
+  REQUIRE(r1.g(0) == doctest::Approx(0.0));
+  REQUIRE(r1.g(1) == doctest::Approx(0.0));
+  REQUIRE(r1.g(2) == doctest::Approx(0.0));
+  REQUIRE(r1.h(0, 0) == doctest::Approx(0.0));
+  REQUIRE(r1.h(0, 1) == doctest::Approx(0.0));
+  REQUIRE(r1.h(0, 2) == doctest::Approx(0.0));
+  REQUIRE(r1.h(1, 1) == doctest::Approx(0.0));
+  REQUIRE(r1.h(1, 2) == doctest::Approx(0.0));
+  REQUIRE(r1.h(2, 2) == doctest::Approx(0.0));
+}
+
 TEST_CASE("Pow") {
   using std::pow;
 
@@ -757,6 +800,34 @@ TEST_CASE("SScalar IncDiv") {
   REQUIRE(r1.d("z") == doctest::Approx(1.00000));
 }
 
+TEST_CASE("SScalar Self multiplication") {
+  auto r1 = s1;
+  r1 *= r1;
+
+  const auto r2 = s1 * s1;
+
+  REQUIRE(r1.f() == doctest::Approx(r2.f()));
+  REQUIRE(r1.d("x") == doctest::Approx(r2.d("x")));
+  REQUIRE(r1.d("y") == doctest::Approx(r2.d("y")));
+  REQUIRE(r1.d("z") == doctest::Approx(r2.d("z")));
+
+  // s1 is 3 with dx = 1, dy = 6 and dz = 4
+  REQUIRE(r1.f() == doctest::Approx(9.0));
+  REQUIRE(r1.d("x") == doctest::Approx(6.0));
+  REQUIRE(r1.d("y") == doctest::Approx(36.0));
+  REQUIRE(r1.d("z") == doctest::Approx(24.0));
+}
+
+TEST_CASE("SScalar Self division") {
+  auto r1 = s1;
+  r1 /= r1;
+
+  REQUIRE(r1.f() == doctest::Approx(1.0));
+  REQUIRE(r1.d("x") == doctest::Approx(0.0));
+  REQUIRE(r1.d("y") == doctest::Approx(0.0));
+  REQUIRE(r1.d("z") == doctest::Approx(0.0));
+}
+
 TEST_CASE("SScalar Pow") {
   using std::pow;