[AUD-006][High] Avoid oversized shift in matrix exponential scaling

[gabrielsantosphilips]

Audit finding

• ID: AUD-006
• Status: Verified defect
• Severity: High
• Confidence: High
• Audited revision: 2f479320d805a1f9f35ebe4afaaeeded48913a94

Problem

Scaling and squaring computes 1 << s as a signed integer before converting to T.

Source:

numerical-toolbox-cpp/numerical/math/MatrixExponential.hpp

Lines 125 to 143 in 2f47932

	{
	const T lg = math::Log2(norm);
	s = static_cast<int>(math::Ceil(lg));
	if (s < 0)
	s = 0;
	}
	const T scale = T{ 1 } / static_cast<T>(1 << s);
	const auto as = a * scale;
	SquareMatrix<T, N> num{};
	SquareMatrix<T, N> den{};
	PadeNumeratorDenominator(as, num, den);
	auto r = SolvePade(den, num);
	for (int k = 0; k < s; ++k)
	r = r * r;

Reproduction

For the finite nilpotent matrix [[0, 1e20], [0, 0]], s=67. UBSan reports a shift exponent too large for 32-bit int, and the result is NaN. The exact result is I + A.

Acceptance criteria

• Replace integer shifting with a bounded floating scaling operation such as Ldexp through the math abstraction.
• Detect unsupported/non-finite scaling ranges and Padé solve failure.
• The large finite nilpotent case returns I + A within a justified float tolerance.
• Add large-norm and overflow-boundary TEST_F cases.