small arithmetic optimization

2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.Format.cs

@@ -0,0 +1,278 @@
+// Printing and parsing of DoubleDoubleSample, by Marek Fiser.
+// Source: https://marekfiser.com/blog/double-double-arithmetic
+// Licensed under CC BY 4.0: https://creativecommons.org/licenses/by/4.0/
+// Changes for this benchmark: copied into this project, kept in the
+// Just.PreciseMath.PerformanceComparison.Reference namespace, and the file and type
+// renamed to DoubleDoubleSample.Format.cs and DoubleDoubleSample.
+// Printing 31 correct decimal digits requires computing them exactly, which is done with
+// BigInteger (part of the BCL, so still no NuGet dependencies).
+using System.Globalization;
+using System.Numerics;
+using System.Text;
+
+namespace Just.PreciseMath.PerformanceComparison.Reference;
+
+public readonly partial struct DoubleDoubleSample {
+
+	/// <summary>Significant digits <see cref="ToString()"/> prints: enough to round-trip any
+	/// value whose two components are contiguous (see <see cref="Parse"/>).</summary>
+	private const int DEFAULT_DIGITS = 36;
+
+
+	/// <summary>Formats with <see cref="DEFAULT_DIGITS"/> significant digits.</summary>
+	public override string ToString() {
+		return ToString(DEFAULT_DIGITS);
+	}
+
+	/// <summary>Formats with the given number of significant decimal digits, correctly rounded.</summary>
+	public string ToString(int significantDigits) {
+        ArgumentOutOfRangeException.ThrowIfLessThan(significantDigits, 1);
+
+        if (!double.IsFinite(m_hi)) {
+			return m_hi.ToString(CultureInfo.InvariantCulture);
+		}
+
+		(BigInteger numerator, BigInteger denominator) = toRational(this);
+		return formatRational(numerator, denominator, significantDigits);
+	}
+
+	/// <summary>
+	/// Formats the exact value. A double-double is a sum of two binary fractions, so it
+	/// always has a finite decimal expansion, but it can be hundreds of digits long when the
+	/// components are far apart.
+	/// </summary>
+	public string ToStringExact() {
+		if (!double.IsFinite(m_hi)) {
+			return m_hi.ToString(CultureInfo.InvariantCulture);
+		}
+
+		(BigInteger numerator, BigInteger denominator) = toRational(this);
+		if (numerator.IsZero) {
+			return "0";
+		}
+
+		// The denominator is a power of two, so multiplying by the matching power of five
+		// turns it into a power of ten, which gives the decimal digits directly.
+		int twos = (int)(denominator.GetBitLength() - 1);
+		BigInteger scaled = BigInteger.Abs(numerator) * BigInteger.Pow(5, twos);
+		string digits = scaled.ToString(CultureInfo.InvariantCulture).PadLeft(twos + 1, '0');
+		string sign = numerator.Sign < 0 ? "-" : "";
+		if (twos == 0) {
+			return sign + digits;
+		}
+
+		string result = sign + digits[..^twos] + "." + digits[^twos..];
+		return result.TrimEnd('0').TrimEnd('.');
+	}
+
+	/// <summary>
+	/// Parses a decimal string to the nearest double-double. The value is built exactly as a
+	/// rational, then rounded once into each component. Round-trips with
+	/// <see cref="ToString()"/> for every value whose components are contiguous;
+	/// a hand-built pair such as (1.0, 1e-300) holds more information than 36 digits can
+	/// represent.
+	/// </summary>
+	public static DoubleDoubleSample Parse(string text) {
+		(BigInteger numerator, BigInteger denominator) = parseRational(text);
+		double hi = toNearestDouble(numerator, denominator);
+		if (!double.IsFinite(hi)) {
+			return FromDouble(hi);
+		}
+
+		// Subtract the leading component exactly and round the remainder into the trailing one.
+		(BigInteger hiNumerator, BigInteger hiDenominator) = toRational(FromDouble(hi));
+		BigInteger restNumerator = numerator * hiDenominator - hiNumerator * denominator;
+		BigInteger restDenominator = denominator * hiDenominator;
+		double lo = toNearestDouble(restNumerator, restDenominator);
+		return FromTwoDoubles(hi, lo);
+	}
+
+	/// <summary>The exact value of a double-double as a rational (the denominator is a power of two).</summary>
+	private static (BigInteger Numerator, BigInteger Denominator) toRational(DoubleDoubleSample value) {
+		(BigInteger hiNumerator, BigInteger hiDenominator) = toRational(value.m_hi);
+		(BigInteger loNumerator, BigInteger loDenominator) = toRational(value.m_lo);
+		return (hiNumerator * loDenominator + loNumerator * hiDenominator, hiDenominator * loDenominator);
+	}
+
+	/// <summary>The exact value of a finite double as a rational, from its bits.</summary>
+	private static (BigInteger Numerator, BigInteger Denominator) toRational(double value) {
+		long bits = BitConverter.DoubleToInt64Bits(value);
+		int exponent = (int)((bits >> 52) & 0x7FF);
+		long mantissa = bits & 0xF_FFFF_FFFF_FFFFL;
+		if (exponent == 0) {
+			exponent = -1074;  // subnormal: no implicit leading bit
+		} else {
+			mantissa |= 1L << 52;
+			exponent -= 1075;
+		}
+
+		BigInteger numerator = bits < 0 ? -mantissa : mantissa;
+		return exponent >= 0
+			? (numerator << exponent, BigInteger.One)
+			: (numerator, BigInteger.One << -exponent);
+	}
+
+	/// <summary>Rounds an exact rational to the nearest double, ties to even.</summary>
+	private static double toNearestDouble(BigInteger numerator, BigInteger denominator) {
+		if (numerator.IsZero) {
+			return 0.0;
+		}
+
+		int sign = numerator.Sign;
+		numerator = BigInteger.Abs(numerator);
+
+		// Scale so that the quotient is in [2^52, 2^54). One integer division then gives all
+		// bits of the significand plus the remainder needed for the rounding.
+		long exponent = numerator.GetBitLength() - denominator.GetBitLength() - 53;
+		if (exponent > 0) {
+			denominator <<= (int)exponent;
+		} else {
+			numerator <<= (int)-exponent;
+		}
+
+		BigInteger quotient = BigInteger.DivRem(numerator, denominator, out BigInteger remainder);
+		if (quotient.GetBitLength() > 53) {
+			// One bit too many. Move it to the exponent and add the dropped bit to the
+			// remainder so that the rounding decision below stays exact.
+			if (!quotient.IsEven) {
+				remainder += denominator;
+			}
+			quotient >>= 1;
+			denominator <<= 1;
+			exponent++;
+		}
+
+		int comparison = (remainder << 1).CompareTo(denominator);
+		if (comparison > 0 || (comparison == 0 && !quotient.IsEven)) {
+			quotient++;
+			if (quotient.GetBitLength() > 53) {
+				quotient >>= 1;
+				exponent++;
+			}
+		}
+
+		return sign * Math.ScaleB((double)quotient, (int)exponent);
+	}
+
+	/// <summary>Parses a decimal string such as "-1.25e-7" into an exact rational.</summary>
+	private static (BigInteger Numerator, BigInteger Denominator) parseRational(string text) {
+		text = text.Trim();
+		if (text.Length == 0) {
+			throw new FormatException("Empty input.");
+		}
+
+		int index = 0;
+		bool negative = text[index] == '-';
+		if (negative || text[index] == '+') {
+			index++;
+		}
+
+		BigInteger mantissa = BigInteger.Zero;
+		int fractionDigits = 0;
+		bool seenDot = false;
+		bool seenDigit = false;
+
+		for (; index < text.Length; index++) {
+			char c = text[index];
+			if (c == '.') {
+				if (seenDot) {
+					throw new FormatException($"Two decimal points in '{text}'.");
+				}
+				seenDot = true;
+			} else if (c is >= '0' and <= '9') {
+				mantissa = mantissa * 10 + (c - '0');
+				seenDigit = true;
+				if (seenDot) {
+					fractionDigits++;
+				}
+			} else if (c is 'e' or 'E') {
+				break;
+			} else {
+				throw new FormatException($"Unexpected character '{c}' in '{text}'.");
+			}
+		}
+
+		if (!seenDigit) {
+			throw new FormatException($"No digits in '{text}'.");
+		}
+
+		int exponent = -fractionDigits;
+		if (index < text.Length) {
+			exponent += int.Parse(text[(index + 1)..], CultureInfo.InvariantCulture);
+		}
+
+		if (negative) {
+			mantissa = -mantissa;
+		}
+
+		return exponent >= 0
+			? (mantissa * BigInteger.Pow(10, exponent), BigInteger.One)
+			: (mantissa, BigInteger.Pow(10, -exponent));
+	}
+
+	/// <summary>Rounds an exact rational to N significant digits and formats them.</summary>
+	private static string formatRational(BigInteger numerator, BigInteger denominator, int digits) {
+		if (numerator.IsZero) {
+			return "0";
+		}
+
+		bool negative = numerator.Sign < 0;
+		numerator = BigInteger.Abs(numerator);
+
+		// Estimate of the decimal exponent of the leading digit from the bit lengths
+		// (log10(2) ~ 0.30103). The loop below corrects it when the estimate is off by one.
+		int exponent10 = (int)Math.Floor((numerator.GetBitLength() - denominator.GetBitLength()) * 0.30103);
+		string text;
+
+		while (true) {
+			int scale = digits - 1 - exponent10;
+			BigInteger scaledNumerator = numerator;
+			BigInteger scaledDenominator = denominator;
+			if (scale >= 0) {
+				scaledNumerator *= BigInteger.Pow(10, scale);
+			} else {
+				scaledDenominator *= BigInteger.Pow(10, -scale);
+			}
+			BigInteger rounded = roundToNearest(scaledNumerator, scaledDenominator);
+			text = rounded.ToString(CultureInfo.InvariantCulture);
+			if (text.Length == digits) {
+				break;
+			}
+			// Off by one either way (estimate too low, or rounding carried into a new digit).
+			exponent10 += text.Length - digits;
+		}
+
+		text = text.TrimEnd('0');
+		if (text.Length == 0) {
+			text = "0";
+		}
+
+		StringBuilder result = new(negative ? "-" : "");
+		if (exponent10 >= -5 && exponent10 < digits) {
+			// Positional notation for the range where it is readable.
+			if (exponent10 >= 0) {
+				text = text.PadRight(exponent10 + 1, '0');
+				result.Append(text[..(exponent10 + 1)]);
+				if (text.Length > exponent10 + 1) {
+					result.Append('.').Append(text[(exponent10 + 1)..]);
+				}
+			} else {
+				result.Append("0.").Append('0', -exponent10 - 1).Append(text);
+			}
+		} else {
+			result.Append(text[0]);
+			if (text.Length > 1) {
+				result.Append('.').Append(text[1..]);
+			}
+			result.Append('e').Append(exponent10.ToString(CultureInfo.InvariantCulture));
+		}
+		
+		return result.ToString();
+	}
+
+	/// <summary>Integer nearest to numerator/denominator, ties away from zero (both positive).</summary>
+	private static BigInteger roundToNearest(BigInteger numerator, BigInteger denominator) {
+		BigInteger quotient = BigInteger.DivRem(numerator, denominator, out BigInteger remainder);
+		return (remainder << 1) >= denominator ? quotient + 1 : quotient;
+	}
+}

2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.cs

@@ -0,0 +1,213 @@
+// Double-double arithmetic in about 200 lines, by Marek Fiser.
+// Source: https://marekfiser.com/blog/double-double-arithmetic
+// Licensed under CC BY 4.0: https://creativecommons.org/licenses/by/4.0/
+// Changes for this benchmark: copied into this project, kept in the
+// Just.PreciseMath.PerformanceComparison.Reference namespace, and the file and type
+// renamed to DoubleDoubleSample.
+// This file depends only on System so that it can be copied into any project. Printing and
+// parsing use BigInteger and are in DoubleDoubleSample.Format.cs.
+namespace Just.PreciseMath.PerformanceComparison.Reference;
+
+/// <summary>
+/// A value stored as the unevaluated sum of two doubles: ~106 significand bits
+/// (~31 decimal digits) at a few times the cost of a double, with no heap allocation.
+///
+/// Invariant: Hi == RoundToNearest(Hi + Lo), i.e. |Lo| is at most half an ulp of Hi.
+/// Every factory and operation below maintains it; comparisons and ToDouble rely on it.
+///
+/// Operations are not correctly rounded but stay within a few units of 2^-106 of the
+/// exact result. The exponent range is the same as a plain double's, and infinities are
+/// not special-cased: the first operation on one returns NaN (the error terms compute
+/// inf - inf).
+///
+/// By Marek Fiser, from marekfiser.com/blog/double-double-arithmetic/ (CC BY 4.0; the file
+/// header carries the full attribution and license notice).
+/// </summary>
+public readonly partial struct DoubleDoubleSample : IEquatable<DoubleDoubleSample>, IComparable<DoubleDoubleSample> {
+
+	private readonly double m_hi;
+	private readonly double m_lo;
+
+	/// <summary>The leading component; equals the value rounded to nearest double.</summary>
+	public double Hi => m_hi;
+
+	/// <summary>The trailing component (the rounding error of <see cref="Hi"/>).</summary>
+	public double Lo => m_lo;
+
+	public bool IsZero => m_hi == 0.0;
+
+	/// <summary>-1, 0 or +1 (0 also for NaN).</summary>
+	public int Sign => m_hi > 0.0 ? 1 : m_hi < 0.0 ? -1 : 0;
+
+	public static DoubleDoubleSample Zero { get; } = new DoubleDoubleSample(0.0, 0.0);
+
+	public static DoubleDoubleSample One { get; } = new DoubleDoubleSample(1.0, 0.0);
+
+
+	private DoubleDoubleSample(double hi, double lo) {
+		m_hi = hi;
+		m_lo = lo;
+	}
+
+
+	/// <summary>Exact conversion from a double.</summary>
+	public static DoubleDoubleSample FromDouble(double value) {
+		return new DoubleDoubleSample(value, 0.0);
+	}
+
+	/// <summary>The exact sum <paramref name="hi"/> + <paramref name="lo"/> of two arbitrary doubles.</summary>
+	public static DoubleDoubleSample FromTwoDoubles(double hi, double lo) {
+		// Renormalizes so that the invariant holds.
+		(double sum, double error) = twoSum(hi, lo);
+		return new DoubleDoubleSample(sum, error);
+	}
+
+	public DoubleDoubleSample Add(DoubleDoubleSample right) {
+		// Add the two components pairwise, then incorporate lowSum and lowError one by one.
+		// Adding lowError can round, so a second quickTwoSum renormalizes the result.
+		(double sum, double error) = twoSum(m_hi, right.m_hi);
+		(double lowSum, double lowError) = twoSum(m_lo, right.m_lo);
+		error += lowSum;
+		(sum, error) = quickTwoSum(sum, error);
+		error += lowError;
+		(sum, error) = quickTwoSum(sum, error);
+		return new DoubleDoubleSample(sum, error);
+	}
+
+	public DoubleDoubleSample Subtract(DoubleDoubleSample right) {
+		return Add(right.Negate());
+	}
+
+	public DoubleDoubleSample Multiply(DoubleDoubleSample right) {
+		// hi*hi exactly, then the three cross terms. Each cross term is at least 2^-53 smaller
+		// than the product, so plain double arithmetic is accurate enough for them.
+		(double product, double error) = twoProduct(m_hi, right.m_hi);
+		error += m_hi * right.m_lo + m_lo * right.m_hi + m_lo * right.m_lo;
+		(product, error) = quickTwoSum(product, error);
+		return new DoubleDoubleSample(product, error);
+	}
+
+	/// <summary>Multiplies by a plain double. Faster than the full double-double product.</summary>
+	public DoubleDoubleSample Multiply(double right) {
+		(double product, double error) = twoProduct(m_hi, right);
+		error += m_lo * right;
+		(product, error) = quickTwoSum(product, error);
+		return new DoubleDoubleSample(product, error);
+	}
+
+	/// <summary>Slower and slightly less accurate than the other operations.</summary>
+	public DoubleDoubleSample Divide(DoubleDoubleSample right) {
+		// Three steps of long division.
+		double quotient1 = m_hi / right.m_hi;
+		DoubleDoubleSample remainder = Subtract(right.Multiply(quotient1));
+		double quotient2 = remainder.m_hi / right.m_hi;
+		remainder = remainder.Subtract(right.Multiply(quotient2));
+		double quotient3 = remainder.m_hi / right.m_hi;
+
+		(double sum, double error) = quickTwoSum(quotient1, quotient2);
+		return new DoubleDoubleSample(sum, error).Add(FromDouble(quotient3));
+	}
+
+	/// <summary>Slightly faster than the general product.</summary>
+	public DoubleDoubleSample Square() {
+		(double product, double error) = twoProduct(m_hi, m_hi);
+		error += 2.0 * m_hi * m_lo + m_lo * m_lo;
+		(product, error) = quickTwoSum(product, error);
+		return new DoubleDoubleSample(product, error);
+	}
+
+	public DoubleDoubleSample Negate() {
+		// Negation never rounds, so this is exact.
+		return new DoubleDoubleSample(-m_hi, -m_lo);
+	}
+
+	public DoubleDoubleSample Abs() {
+		return Sign < 0 ? Negate() : this;
+	}
+
+	/// <summary>Slightly less accurate than the other operations. Negative and NaN inputs
+	/// return NaN, as Math.Sqrt does.</summary>
+	public DoubleDoubleSample Sqrt() {
+		if (!(m_hi > 0.0)) {           // Not "m_hi <= 0.0", because this form is also true for NaN.
+			return FromDouble(Math.Sqrt(m_hi));
+		}
+		// One Newton step from the double square root.
+		double approx = Math.Sqrt(m_hi);
+		DoubleDoubleSample residual = Subtract(FromDouble(approx).Square());
+		double correction = residual.m_hi / (approx + approx);
+		(double hi, double lo) = quickTwoSum(approx, correction);
+		return new DoubleDoubleSample(hi, lo);
+	}
+
+	public double ToDouble() {
+		return m_hi;
+	}
+
+	public int CompareTo(DoubleDoubleSample other) {
+		// Correct only for normalized pairs.
+		int hiComparison = m_hi.CompareTo(other.m_hi);
+		return hiComparison != 0 ? hiComparison : m_lo.CompareTo(other.m_lo);
+	}
+
+	public bool Equals(DoubleDoubleSample other) {
+		return m_hi == other.m_hi && m_lo == other.m_lo;
+	}
+
+	public override bool Equals(object? obj) {
+		return obj is DoubleDoubleSample other && Equals(other);
+	}
+
+	public override int GetHashCode() {
+		return HashCode.Combine(m_hi, m_lo);
+	}
+
+	public static DoubleDoubleSample operator +(DoubleDoubleSample left, DoubleDoubleSample right) => left.Add(right);
+	public static DoubleDoubleSample operator -(DoubleDoubleSample left, DoubleDoubleSample right) => left.Subtract(right);
+	public static DoubleDoubleSample operator *(DoubleDoubleSample left, DoubleDoubleSample right) => left.Multiply(right);
+	public static DoubleDoubleSample operator *(DoubleDoubleSample left, double right) => left.Multiply(right);
+	public static DoubleDoubleSample operator /(DoubleDoubleSample left, DoubleDoubleSample right) => left.Divide(right);
+	public static DoubleDoubleSample operator -(DoubleDoubleSample value) => value.Negate();
+
+	public static bool operator ==(DoubleDoubleSample left, DoubleDoubleSample right) => left.Equals(right);
+	public static bool operator !=(DoubleDoubleSample left, DoubleDoubleSample right) => !left.Equals(right);
+	public static bool operator <(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) < 0;
+	public static bool operator >(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) > 0;
+	public static bool operator <=(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) <= 0;
+	public static bool operator >=(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) >= 0;
+
+	// Conversion from a double is exact, so it is implicit. Conversion to a double drops the
+	// low component, so it is explicit.
+	public static implicit operator DoubleDoubleSample(double value) => FromDouble(value);
+	public static explicit operator double(DoubleDoubleSample value) => value.ToDouble();
+
+	/// <summary>Knuth's two-sum: sum + error == a + b exactly, for any two doubles.</summary>
+	private static (double Sum, double Error) twoSum(double a, double b) {
+		// The rounding error of a floating-point addition is itself representable as a double.
+		// The intermediate results below may round, but their rounding errors cancel out.
+		double sum = a + b;
+		double bKept = sum - a;        // The part of b that made it into the sum.
+		double aKept = sum - bKept;    // The part of a that made it into the sum.
+		double aLost = a - aKept;
+		double bLost = b - bKept;
+		return (sum, aLost + bLost);
+	}
+
+	/// <summary>Dekker's fast two-sum: the same guarantee as twoSum, valid only when
+	/// |a| >= |b| (or a == 0).</summary>
+	private static (double Sum, double Error) quickTwoSum(double a, double b) {
+		double sum = a + b;
+		double bKept = sum - a;
+		return (sum, b - bKept);       // a lost nothing, so only b's loss remains.
+	}
+
+	/// <summary>Exact product: product + error == a * b exactly, as long as the product is a
+	/// normal double. Below about 2e-292 the error no longer fits a double.</summary>
+	private static (double Product, double Error) twoProduct(double a, double b) {
+		// The fused multiply-add computes a*b - product with a single rounding, which is the
+		// low half of the 106-bit product. On a CPU without an FMA instruction,
+		// Math.FusedMultiplyAdd computes the same result in software, much slower.
+		double product = a * b;
+		double error = Math.FusedMultiplyAdd(a, b, -product);
+		return (product, error);
+	}
+}