Cap Pow and Exp digits by the precision of the double exponent [patch]

SignificantNumber.Test/SignificantNumberPowExpPrecisionTests.cs

@@ -0,0 +1,52 @@
+// Copyright (c) 2023-2026 ktsu-dev contributors
+
+namespace SignificantNumber.Test;
+
+using System.Globalization;
+using ktsu.PreciseNumber;
+using ktsu.SignificantNumber;
+
+/// <summary>
+/// Covers the significant digits <see cref="SignificantNumber.Pow"/> and <see cref="SignificantNumber.Exp"/> report
+/// when the result is computed as e^y in a <see cref="double"/> and |y| is large enough to cost digits.
+/// </summary>
+[TestClass]
+public class SignificantNumberPowExpPrecisionTests
+{
+	private static SignificantNumber Parse(string text) =>
+		SignificantNumber.Parse(text, NumberStyles.Float, CultureInfo.InvariantCulture);
+
+	private static PreciseNumber Precise(string text) =>
+		PreciseNumber.Parse(text, NumberStyles.Float, CultureInfo.InvariantCulture);
+
+	/// <summary>
+	/// Asserts that every digit <paramref name="result"/> reports as significant matches <paramref name="trueValue"/>.
+	/// </summary>
+	private static void AssertReportedDigitsAreCorrect(string trueValue, SignificantNumber result) =>
+		Assert.AreEqual(
+			Precise(trueValue).ToSignificantNumber(result.SignificantDigits),
+			result,
+			$"{result} reports {result.SignificantDigits} significant digits of {trueValue}");
+
+	// Reference values from Python decimal at 50 digits.
+	[TestMethod]
+	[DataRow("700.12345678901234567890", "1.1475032771153835242760515951099859734230596995731E+304")]
+	[DataRow("-700.12345678901234567890", "8.7145720621715330952462542732455202987551213258088E-305")]
+	public void Exp_LargeNonIntegerPower_ReportsOnlyCorrectDigits(string power, string trueValue) =>
+		AssertReportedDigitsAreCorrect(trueValue, SignificantNumber.Exp(Precise(power)));
+
+	[TestMethod]
+	[DataRow("3.0000000000000000001", "200.50000000000000000001", "4.6005692393404968725867994056284596535772404109326E+95")]
+	[DataRow("3.0000000000000000001", "2.5000000000000000001", "15.588457268119895644758622250842932011355302787297")]
+	[DataRow("1.2345678901234567890", "2000", "1.0714068665496473862896950588742712767980475529713E+183")]
+	public void Pow_LargeLogarithmOfResult_ReportsOnlyCorrectDigits(string baseNumber, string power, string trueValue) =>
+		AssertReportedDigitsAreCorrect(trueValue, Parse(baseNumber).Pow(Precise(power)));
+
+	[TestMethod]
+	public void Exp_PowerNearSevenHundred_ReportsTwelveSignificantDigits() =>
+		Assert.AreEqual(12, SignificantNumber.Exp(Precise("700.12345678901234567890")).SignificantDigits);
+
+	[TestMethod]
+	public void Pow_SmallLogarithmOfResult_KeepsTheDoublesPrecision() =>
+		Assert.AreEqual(15, Parse("1.2345678901234567890").Pow(Precise("0.5000000000000000001")).SignificantDigits);
+}

SignificantNumber/SignificantNumber.cs

@@ -644,6 +644,9 @@
 	/// <returns>|<paramref name="value"/>|^<paramref name="power"/>, negated when <paramref name="isNegative"/> is set.</returns>
 	private static SignificantNumber PowByLogarithm(PreciseNumber value, PreciseNumber power, bool isNegative, int significantDigits)
 	{
+		// |x|^p = e^(p × ln|x|), so the digits are bounded by how precisely e^(p × ln|x|) is known.
+		significantDigits = CapToExponentPrecision(significantDigits, Log10OfMagnitude(value) * Math.Log(10) * power.To<double>());
+
 		double logValue = Math.Log(Math.Abs(value.To<double>()));
 		double magnitude = Math.Exp(logValue * power.To<double>());
 		if (double.IsNormal(magnitude))
@@ -680,6 +683,7 @@
 		int significantDigits = PreciseNumber.IsInteger(power)
 			? int.Min(E.SignificantDigits, DoubleSignificantDigits)
 			: LowestSignificantDigits(E, power);
+		significantDigits = CapToExponentPrecision(significantDigits, power.To<double>());
 
 		double result = Math.Exp(power.To<double>());
 		if (double.IsNormal(result))
@@ -691,6 +695,25 @@
 		return FromLog10OfMagnitude(power.To<double>() * Math.Log10(Math.E), isNegative: false, significantDigits);
 	}
 
+	/// <summary>
+	/// Limits a significant digit count to the digits that e^<paramref name="exponent"/> computed in a
+	/// <see cref="double"/> can be trusted to. The exponent carries a relative error of about 2⁻⁵³, which
+	/// becomes an absolute error of |exponent| × 2⁻⁵³ in it, and so a relative error of that size in the result:
+	/// every power of ten in |exponent| costs one of the double's significant digits.
+	/// </summary>
+	/// <param name="significantDigits">The significant digit count to limit.</param>
+	/// <param name="exponent">The natural logarithm of the magnitude of the result.</param>
+	/// <returns>The lower of <paramref name="significantDigits"/> and the digits the double result determines, and at least 1.</returns>
+	private static int CapToExponentPrecision(int significantDigits, double exponent)
+	{
+		double magnitude = Math.Abs(exponent);
+		int lostDigits = !double.IsFinite(magnitude)
+			? DoubleSignificantDigits
+			: magnitude > 1 ? (int)Math.Ceiling(Math.Log10(magnitude)) : 0;
+
+		return int.Max(1, int.Min(significantDigits, DoubleSignificantDigits - lostDigits));
+	}
+
 	/// <summary>
 	/// Computes the base-10 logarithm of the magnitude of a nonzero number from its significand and exponent, so
 	/// that it works for magnitudes outside the range of a <see cref="double"/>.