Compute Pow and Exp beyond the range of a double [patch]

SignificantNumber.Test/SignificantNumberRangeTests.cs

@@ -0,0 +1,72 @@
+// Copyright (c) 2023-2026 ktsu-dev contributors
+
+namespace SignificantNumber.Test;
+
+using System.Globalization;
+using ktsu.PreciseNumber;
+using ktsu.SignificantNumber;
+
+/// <summary>
+/// Covers <see cref="SignificantNumber.Pow"/> and <see cref="SignificantNumber.Exp"/> for results outside the range
+/// of a <see cref="double"/>.
+/// </summary>
+[TestClass]
+public class SignificantNumberRangeTests
+{
+	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);
+
+	[TestMethod]
+	[DataRow("10", "2000", "1E2000")]
+	[DataRow("1.5", "2000", "1.5E352")]
+	[DataRow("1E400", "0.5", "1E200")]
+	[DataRow("-10", "2001", "-1E2001")]
+	public void Pow_ResultAboveDoubleRange_ReturnsMagnitude(string baseNumber, string power, string expected) =>
+		Assert.AreEqual(Parse(expected), Parse(baseNumber).Pow(Precise(power)));
+
+	[TestMethod]
+	[DataRow("2", "-1100", "7E-332")]
+	[DataRow("1E-400", "0.5", "1E-200")]
+	[DataRow("1E400", "-0.5", "1E-200")]
+	[DataRow("-2", "-1101", "-4E-332")]
+	public void Pow_ResultBelowDoubleRange_ReturnsNonzeroMagnitude(string baseNumber, string power, string expected) =>
+		Assert.AreEqual(Parse(expected), Parse(baseNumber).Pow(Precise(power)));
+
+	[TestMethod]
+	public void Exp_ResultAboveDoubleRange_ReturnsMagnitude()
+	{
+		// e^710 = 2.2339...E308 and e^1000 = 1.9700711140170469...E434
+		SignificantNumber exp710 = SignificantNumber.Exp(Precise("710"));
+		SignificantNumber exp1000 = SignificantNumber.Exp(Precise("1000"));
+
+		Assert.AreEqual(308, exp710.Exponent + exp710.SignificantDigits - 1);
+		Assert.AreEqual(434, exp1000.Exponent + exp1000.SignificantDigits - 1);
+		Assert.AreEqual(Parse("1.97007111402E434"), exp1000);
+	}
+
+	[TestMethod]
+	public void Exp_ResultBelowDoubleRange_ReturnsNonzeroMagnitude()
+	{
+		// e^-1000 = 5.0759588975494567...E-435
+		SignificantNumber result = SignificantNumber.Exp(Precise("-1000"));
+
+		Assert.AreEqual(Parse("5.07595889755E-435"), result);
+	}
+
+	[TestMethod]
+	public void Pow_PowerOfTenBeyondExponentRange_ThrowsOverflowException()
+	{
+		Assert.ThrowsExactly<OverflowException>(() => Parse("10").Pow(Precise("1E300")));
+		Assert.ThrowsExactly<OverflowException>(() => Parse("10").Pow(Precise("-1E300")));
+	}
+
+	[TestMethod]
+	public void Exp_PowerBeyondExponentRange_ThrowsOverflowException()
+	{
+		Assert.ThrowsExactly<OverflowException>(() => SignificantNumber.Exp(Precise("1E300")));
+		Assert.ThrowsExactly<OverflowException>(() => SignificantNumber.Exp(Precise("-1E300")));
+	}
+}

SignificantNumber/SignificantNumber.cs

@@ -563,6 +563,7 @@ internal static bool DoesImplementGenericInterface(Type type, Type genericInterf
 	/// <returns>The current number raised to <paramref name="power"/>, rounded to the fewest significant digits of the two.</returns>
 	/// <exception cref="DivideByZeroException">Thrown when the current number is zero and <paramref name="power"/> is negative.</exception>
 	/// <exception cref="ArgumentException">Thrown when the current number is negative and <paramref name="power"/> is not an integer, as the result is not a real number.</exception>
+	/// <exception cref="OverflowException">Thrown when the result's power of ten is too large or too small to represent.</exception>
 	public SignificantNumber Pow(PreciseNumber power)
 	{
 		if (Equal(power, Zero))
@@ -603,19 +604,41 @@ public SignificantNumber Pow(PreciseNumber power)
 			? int.Min(Value.SignificantDigits, DoubleSignificantDigits)
 			: LowestSignificantDigits(this, power);
 
-		// Use logarithm and exponential to support decimal powers. This computes |x|^p, so the
-		// sign of a negative base is restored for odd integer powers.
-		double logValue = Math.Log(Math.Abs(Value.To<double>()));
+		// This computes |x|^p, so the sign of a negative base is restored for odd integer powers.
+		bool isNegativeResult = isNegativeBase && PreciseNumber.IsOddInteger(power);
+		return PowByLogarithm(Value, power, isNegativeResult, significantDigits);
+	}
+
+	/// <summary>
+	/// Raises the magnitude of a nonzero number to a power by logarithm and exponential, which supports decimal powers.
+	/// </summary>
+	/// <param name="value">The nonzero base.</param>
+	/// <param name="power">The power to raise the magnitude of <paramref name="value"/> to.</param>
+	/// <param name="isNegative">Whether the result is negative.</param>
+	/// <param name="significantDigits">The number of significant digits the result should have.</param>
+	/// <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)
+	{
+		double logValue = Math.Log(Math.Abs(value.To<double>()));
 		double magnitude = Math.Exp(logValue * power.To<double>());
-		double result = isNegativeBase && PreciseNumber.IsOddInteger(power) ? -magnitude : magnitude;
-		return result.ToSignificantNumber(significantDigits);
+		if (double.IsNormal(magnitude))
+		{
+			double result = isNegative ? -magnitude : magnitude;
+			return result.ToSignificantNumber(significantDigits);
+		}
+
+		// The base or the result is outside the range of a double, so work in base-10 logarithms
+		// built from the base's significand and exponent, which never pass through a double.
+		double log10Magnitude = Log10OfMagnitude(value) * power.To<double>();
+		return FromLog10OfMagnitude(log10Magnitude, isNegative, significantDigits);
 	}
 
 	/// <summary>
 	/// Returns the result of raising e to the specified power.
 	/// </summary>
 	/// <param name="power">The power to raise e to.</param>
 	/// <returns>e raised to <paramref name="power"/>, rounded to the fewest significant digits of the two.</returns>
+	/// <exception cref="OverflowException">Thrown when the result's power of ten is too large or too small to represent.</exception>
 	public static SignificantNumber Exp(PreciseNumber power)
 	{
 		if (Equal(power, Zero))
@@ -633,8 +656,60 @@ public static SignificantNumber Exp(PreciseNumber power)
 			? int.Min(E.SignificantDigits, DoubleSignificantDigits)
 			: LowestSignificantDigits(E, power);
 
-		return Math.Exp(power.To<double>())
-			.ToSignificantNumber(significantDigits);
+		double result = Math.Exp(power.To<double>());
+		if (double.IsNormal(result))
+		{
+			return result.ToSignificantNumber(significantDigits);
+		}
+
+		// The result is outside the range of a double: e^p = 10^(p × log10(e)).
+		return FromLog10OfMagnitude(power.To<double>() * Math.Log10(Math.E), isNegative: false, significantDigits);
+	}
+
+	/// <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"/>.
+	/// </summary>
+	/// <param name="value">A nonzero number.</param>
+	/// <returns>log10(|<paramref name="value"/>|).</returns>
+	private static double Log10OfMagnitude(PreciseNumber value) =>
+		BigInteger.Log10(BigInteger.Abs(value.Significand)) + value.Exponent;
+
+	/// <summary>
+	/// Builds a number from the base-10 logarithm of its magnitude, keeping the power of ten as an integer so that
+	/// the result can lie outside the range of a <see cref="double"/>.
+	/// </summary>
+	/// <param name="log10Magnitude">log10 of the magnitude of the result.</param>
+	/// <param name="isNegative">Whether the result is negative.</param>
+	/// <param name="significantDigits">The number of significant digits the result should have.</param>
+	/// <returns>
+	/// The result, rounded to <paramref name="significantDigits"/>, or to fewer digits when the logarithm cannot
+	/// determine that many: a double carries about 15 significant digits, and the integer part of the logarithm
+	/// uses some of them.
+	/// </returns>
+	/// <exception cref="OverflowException">Thrown when the result's power of ten is too large or too small to represent.</exception>
+	private static SignificantNumber FromLog10OfMagnitude(double log10Magnitude, bool isNegative, int significantDigits)
+	{
+		double wholePart = Math.Floor(log10Magnitude);
+		if (!double.IsFinite(wholePart) || Math.Abs(wholePart) > MaxLog10Magnitude)
+		{
+			throw new OverflowException(log10Magnitude > 0
+				? "The result is too large to represent."
+				: "The result is too small to represent.");
+		}
+
+		// 10^fraction is in [1, 10), so scaling it by 10^16 gives a 17-digit significand that a double holds exactly.
+		double mantissa = Math.Pow(10, log10Magnitude - wholePart);
+		BigInteger significand = new(Math.Round(mantissa * 1e16));
+		int powerOfTen = (int)wholePart;
+		int exponent = powerOfTen - 16;
+
+		int integerDigits = int.Abs(powerOfTen).ToString(CultureInfo.InvariantCulture).Length;
+		int determinedDigits = int.Max(1, DoubleSignificantDigits - integerDigits);
+
+		return CreateFromComponents(exponent, isNegative ? -significand : significand)
+			.Value
+			.ToSignificantNumber(int.Min(significantDigits, determinedDigits));
 	}
 
 	/// <summary>
@@ -648,6 +723,12 @@ public static SignificantNumber Exp(PreciseNumber power)
 	/// </summary>
 	private const int MaxExactIntegerPower = 1024;
 
+	/// <summary>
+	/// The largest power of ten, in either direction, that <see cref="Pow"/> and <see cref="Exp"/> produce. It keeps the
+	/// exponent of the result well inside the range of an <see cref="int"/>.
+	/// </summary>
+	private const int MaxLog10Magnitude = 1_000_000_000;
+
 	/// <summary>
 	/// Gets the value of an integer <see cref="PreciseNumber"/> as a <see cref="BigInteger"/>.
 	/// </summary>