// 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);
	}
}