fix: handle aliased operands in the in-place operators

[oberbichler]

Problem

x *= x and x /= x compute derivatives from data they have already overwritten, because b and *this are the same storage:

	mechanism
DDScalar::operator*=	the Hessian loop reads b.m_data[1 + j] after the gradient loop has replaced it
DDScalar::operator/=	binary() writes into m_data, which is also its b
SScalar::operator*=	the second loop reads b.m_d after the first has scaled it
SScalar::operator/=	same

Measured on x = 3 with dx = 1:

	result	expected
DDScalar x *= x	h00 = 12	h00 = 2	wrong
DDScalar x /= x	0	0	accidentally right
SScalar x *= x	dx = 12	dx = 6	wrong
SScalar x /= x	dx = 0.222	dx = 0	wrong

Two corrections to what I claimed when filing this:

• Only the second-order part of DDScalar::operator*= is affected. The gradient loop reads b.m_data[i] in the same expression that writes it, so it still sees the original value and produces the correct 2·f·g. That is why the first-order types come out right.
• DDScalar::operator/= is not observably broken. When both operands alias, the gradient loop zeroes every slot before the Hessian loop reads them — and the correct Hessian of x / x is zero as well, so the stale data happens to give the right answer. It is guarded anyway: that only holds because aliasing forces the trivial result, and a change to the loop order would break it silently.

+= and -= are correct in both classes and are left alone. They read and write the same element within one expression, which yields 2a and 0 as it should.

Fix

All four delegate to the binary operator, which builds a fresh result and therefore cannot read overwritten data:

// aliased operands would be read again after being overwritten below
if (this == &b) {
  return *this = *this * b;
}

The extra copy only happens for the aliased call. Cost of the pointer comparison on the hot path — 2M iterations of a *= b; a /= b, -O3 -DNDEBUG, static order 2 size 8, three alternating runs:

run	with guard	without
1	177.3 ns/op	178.5 ns/op
2	177.0 ns/op	175.8 ns/op
3	179.5 ns/op	174.2 ns/op

Noise. Dynamic size 8 and 32 likewise (152.4 vs 152.2, 754.0 vs 758.3), with identical accumulated results.

Tests

Written first. The oracle is the binary operator — the invariant being that in-place and binary agree — plus analytic anchors so the test cannot pass by both paths breaking the same way: x² at x = 3 with g0 = 1 has f = 9, g0 = 6, h00 = 2, and x / x is 1 with every derivative zero.

C++ — Self multiplication, Self division, SScalar Self multiplication, SScalar Self division. RED:

test.cpp:303: FATAL ERROR: REQUIRE( r1.h(0, 0) == doctest::Approx(r2.h(0, 0)) ) is NOT correct!
test.cpp:810: FATAL ERROR: REQUIRE( r1.d("x") == doctest::Approx(r2.d("x")) ) is NOT correct!
test.cpp:826: FATAL ERROR: REQUIRE( r1.d("x") == doctest::Approx(0.0) ) is NOT correct!

Python — test_self_multiplication and test_self_division over all four {order 1, 2} × {static, dynamic} contexts, plus the two SScalar cases. RED failed on exactly the four combinations predicted above, which is what confirmed the two corrections: test_self_multiplication failed for 2nd order only, and test_self_division passed for all four DDScalar variants.

After: C++ 46/46 test cases, 456 assertions · Python 176 passed · clang-format --Werror and ruff format --check clean.

Notes

Based on main with #73 and #74 merged. Independent of #75 (still open) — and once #75 lands, this code is covered by the sanitizer job as well.

Last item from the review's P0 list:

• resize() silently scrambles the Hessian layout for order 2 (pad_right() does it correctly). fix: validate indices and sizes instead of relying on asserts #73 only stopped it from accepting a negative size.