diff --git a/0-source/Just.PreciseMath/DDMath.cs b/0-source/Just.PreciseMath/DDMath.cs
new file mode 100644
index 0000000..99b3b7b
--- /dev/null
+++ b/0-source/Just.PreciseMath/DDMath.cs
@@ -0,0 +1,57 @@
+namespace Just.PreciseMath;
+
+/// Provides mathematical functions for normalized double-double values.
+/// Results use fixed-size double-double precision, not arbitrary precision.
+public static class DDMath
+{
+ /// Returns the absolute value without discarding the low component.
+ /// Both signs of zero become positive zero; NaN is canonicalized and either infinity becomes positive infinity.
+ [Pure]
+ public static DoubleDouble Abs(DoubleDouble value)
+ {
+ return DoubleDouble.Abs(value);
+ }
+
+ /// Returns the nonnegative square root with double-double precision.
+ ///
+ /// Uses power-of-two scaling and an FMA-based Newton correction. Results are
+ /// approximate, not guaranteed correctly rounded; finite positive inputs are
+ /// tested against a relative error bound of 2^-100 across the binary64 range.
+ /// Signed zero and positive infinity are preserved. Negative nonzero values
+ /// (including negative infinity) and NaN return canonical NaN.
+ ///
+ [Pure]
+ public static DoubleDouble Sqrt(DoubleDouble value)
+ {
+ if (double.IsNaN(value.High) || value.High < 0.0)
+ {
+ return DoubleDouble.NaN;
+ }
+ if (value.High == 0.0 || double.IsPositiveInfinity(value.High))
+ {
+ return value;
+ }
+
+ // Choose an even exponent, rounding negative odd exponents down too.
+ // For positive finite normalized inputs, exponent is in [-1074, 1022]
+ // and the scaled high is in [1, 4). Scaling the high is exact; any
+ // underflow in a very sparse low is far below the relative error bound.
+ int exponent = Math.ILogB(value.High) & ~1;
+ double high = Math.ScaleB(value.High, -exponent);
+ double low = Math.ScaleB(value.Low, -exponent);
+ double estimate = Math.Sqrt(high);
+
+ // Estimate is in [1, 2]. FMA avoids rounding its square before
+ // cancellation. One Newton correction reduces O(u) error to O(u^2),
+ // where u=2^-53; both components of the input contribute to the residual.
+ double residual = Math.FusedMultiplyAdd(-estimate, estimate, high) + low;
+ double correction = residual / (2.0 * estimate);
+ DoubleDouble root = DoubleDouble.FromComponents(estimate, correction);
+
+ // The root high stays normal, even if rounding reaches the next binade.
+ // Normalize again to canonicalize a zero low, including any underflow
+ // of an exceptionally sparse correction during rescaling.
+ int rootExponent = exponent / 2;
+ return DoubleDouble.FromComponents(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent));
+ }
+}
diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs
new file mode 100644
index 0000000..a9cf5c4
--- /dev/null
+++ b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs
@@ -0,0 +1,161 @@
+using System.Numerics;
+using Shouldly;
+using Xunit;
+
+namespace Just.PreciseMath.Tests;
+
+public class PreciseMathSqrtTests
+{
+ [Theory]
+ [InlineData(1.0, 1.0)]
+ [InlineData(4.0, 2.0)]
+ [InlineData(9.0, 3.0)]
+ [InlineData(0.25, 0.5)]
+ [InlineData(2.25, 1.5)]
+ public void ExactBinarySquaresHaveExactRoots(double square, double root)
+ {
+ DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(square));
+ actual.High.ShouldBe(root);
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
+ }
+
+ [Fact]
+ public void PowersOfFourHaveExactRootsAcrossTheFiniteRange()
+ {
+ // sqrt(2^(2k)) = 2^k exactly, including the minimum subnormal input.
+ for (int exponent = -1074; exponent <= 1022; exponent += 2)
+ {
+ DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(Math.ScaleB(1.0, exponent)));
+ actual.High.ShouldBe(Math.ScaleB(1.0, exponent / 2));
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
+ }
+ }
+
+ [Theory]
+ [InlineData(2.0)]
+ [InlineData(3.0)]
+ [InlineData(5.0)]
+ [InlineData(1e-308)]
+ [InlineData(1e308)]
+ public void IrrationalRootsRetainMoreThanBinary64Precision(double value)
+ {
+ DoubleDouble input = new(value);
+ DoubleDouble actual = DDMath.Sqrt(input);
+ actual.Low.ShouldNotBe(0.0);
+ AssertSqrtBound(input, actual);
+ }
+
+ [Fact]
+ public void LowComponentsOnEitherSideOfAnExactSquareAffectTheRoot()
+ {
+ foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(1.0, low);
+ DoubleDouble actual = DDMath.Sqrt(input);
+ actual.High.ShouldBe(1.0);
+ (Math.Sign(actual.Low) == Math.Sign(low)).ShouldBeTrue();
+ AssertSqrtBound(input, actual);
+ }
+ }
+
+ [Fact]
+ public void SpecialValuesMatchBinary64AndRemainCanonical()
+ {
+ double[] values = [0.0, -0.0, double.PositiveInfinity, double.NegativeInfinity,
+ double.NaN, -1.0, -double.Epsilon, double.MinValue];
+ foreach (double value in values)
+ {
+ DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(value));
+ double expected = Math.Sqrt(value);
+ // Construction canonicalizes all NaNs to double.NaN, including
+ // domain errors; zeros must instead be compared by their sign bits.
+ BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(
+ BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected));
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
+ DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+ }
+
+ DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54));
+ DoubleDouble.IsNaN(DDMath.Sqrt(negative)).ShouldBeTrue();
+ }
+
+ [Fact]
+ public void PositiveFiniteInputsMeetExactRelativeErrorBound()
+ {
+ // Every binary64 exponent, including subnormal input exponents. Test zero,
+ // dense, and sparse lows of either sign, without relying on DD arithmetic
+ // to compute the expected square root or its square.
+ Random random = new(314159);
+ for (int exponent = -1074; exponent <= 1023; ++exponent)
+ {
+ double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
+ double low = Math.ScaleB(random.NextDouble(), exponent - 54);
+ double power = Math.ScaleB(1.0, exponent);
+ foreach (double boundary in new[] { Math.BitDecrement(power), power, Math.BitIncrement(power) })
+ {
+ if (boundary > 0.0)
+ {
+ DoubleDouble input = new(boundary);
+ AssertSqrtBound(input, DDMath.Sqrt(input));
+ }
+ }
+ foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(high, residual);
+ if (input.High > 0.0)
+ {
+ AssertSqrtBound(input, DDMath.Sqrt(input));
+ }
+ }
+ }
+
+ // Pin the tiny-input precision loss described in the optional legacy review,
+ // and guard against intermediate squaring overflow near the finite limit.
+ DoubleDouble[] boundaries = [new(1e-308), new(double.Epsilon),
+ new(Math.BitDecrement(Math.ScaleB(1.0, -1022))), new(Math.ScaleB(1.0, -1022)),
+ new(Math.BitIncrement(Math.ScaleB(1.0, -1022))), new(double.MaxValue),
+ DoubleDouble.FromComponents(double.MaxValue, Math.ScaleB(1.0, 969)),
+ DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 969)),
+ DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53)),
+ DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54)),
+ DoubleDouble.FromComponents(4.0, Math.ScaleB(1.0, -51)),
+ DoubleDouble.FromComponents(4.0, -Math.ScaleB(1.0, -52))];
+ foreach (DoubleDouble input in boundaries)
+ {
+ AssertSqrtBound(input, DDMath.Sqrt(input));
+ }
+ }
+
+ private static void AssertSqrtBound(DoubleDouble input, DoubleDouble actual)
+ {
+ DoubleDouble.IsFinite(actual).ShouldBeTrue();
+ DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+ (actual.High > 0.0).ShouldBeTrue();
+
+ // Independent exact rational oracle: x = X*2^-1074, y = Y*2^-1074.
+ // For positive x and y, |y/sqrt(x) - 1| <= t iff
+ // x*(1-t)^2 <= y^2 <= x*(1+t)^2. With t=2^-100 we
+ // cross-multiply to integers: no rounded sqrt, product, or DD conversion.
+ BigInteger x = (Units(input.High) + Units(input.Low)) << 1074;
+ BigInteger y = Units(actual.High) + Units(actual.Low);
+ BigInteger scale = BigInteger.One << 100;
+ BigInteger squared = (y * y) << 200;
+ BigInteger lower = x * (scale - 1) * (scale - 1);
+ BigInteger upper = x * (scale + 1) * (scale + 1);
+ (squared >= lower && squared <= upper).ShouldBeTrue(
+ $"Sqrt bound failed for ({input.High:R}, {input.Low:R}): ({actual.High:R}, {actual.Low:R})");
+ }
+
+ private static BigInteger Units(double value)
+ {
+ long bits = BitConverter.DoubleToInt64Bits(value);
+ int exponent = (int)((bits >> 52) & 0x7ff);
+ BigInteger significand = bits & 0xfffffffffffffL;
+ if (exponent != 0)
+ {
+ significand += BigInteger.One << 52;
+ significand <<= exponent - 1;
+ }
+ return bits < 0 ? -significand : significand;
+ }
+}
diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs
new file mode 100644
index 0000000..66ca109
--- /dev/null
+++ b/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs
@@ -0,0 +1,35 @@
+using Shouldly;
+using Xunit;
+
+namespace Just.PreciseMath.Tests;
+
+public class PreciseMathTests
+{
+ [Fact]
+ public void AbsPreservesBothComponentsAndCanonicalSpecialValues()
+ {
+ double[] scalars = [0.0, -0.0, 1.0, -1.0, double.Epsilon, -double.Epsilon,
+ double.MaxValue, double.MinValue, double.PositiveInfinity, double.NegativeInfinity, double.NaN];
+ foreach (double scalar in scalars)
+ {
+ // Binary64 is an independent reference when the input has no residual.
+ DoubleDouble actual = DDMath.Abs(new DoubleDouble(scalar));
+ AssertBits(actual, double.IsNaN(scalar) ? double.NaN : Math.Abs(scalar), 0.0);
+ DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+ }
+
+ foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54), double.Epsilon, -double.Epsilon })
+ {
+ DoubleDouble positive = DoubleDouble.FromComponents(1.0, low);
+ DoubleDouble negative = DoubleDouble.FromComponents(-1.0, -low);
+ AssertBits(DDMath.Abs(positive), 1.0, low);
+ AssertBits(DDMath.Abs(negative), 1.0, low);
+ }
+ }
+
+ private static void AssertBits(DoubleDouble actual, double high, double low)
+ {
+ BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(high));
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low));
+ }
+}
diff --git a/README.md b/README.md
index 3fcb474..7df7270 100644
--- a/README.md
+++ b/README.md
@@ -75,8 +75,32 @@ DoubleDouble quarterTurn = DoubleDouble.PiOver2;
The factors avoid deriving constants at runtime; the multiplication itself remains
approximate DD arithmetic, so conversions are not guaranteed exact round trips.
The list is mathematical and dimensionless, not a table of unit-dependent physical
-constants. Precomputed roots/logarithms do not imply general `Sqrt`/`Log` functions
-are implemented.
+constants. Precomputed logarithms do not imply a general `Log` function is implemented.
+
+## Mathematical functions
+
+The initial `DDMath` static class provides:
+
+- `Abs(DoubleDouble)`: preserves both components, maps either signed zero to
+ positive zero and either infinity to positive infinity, and returns canonical NaN.
+ It shares the existing `DoubleDouble.Abs` implementation.
+- `Sqrt(DoubleDouble)`: uses power-of-two scaling and an FMA-based Newton correction
+ to retain extended precision, including for subnormal inputs, without squaring
+ an unscaled estimate near the exponent limits. Signed zero and positive infinity
+ are preserved; negative nonzero inputs and NaN return canonical NaN.
+
+```csharp
+using Just.PreciseMath;
+
+DoubleDouble root = DDMath.Sqrt(new DoubleDouble(2.0));
+DoubleDouble magnitude = DDMath.Abs(-root);
+```
+
+Square-root tests compare the exact component sum against a `2^-100` relative
+error bound using integer inequalities. They include samples at every binary64
+exponent, boundary neighbors, both signs of the low component, and exact binary
+squares. This is a tested approximate-accuracy contract, not exhaustive coverage
+of all component pairs or a guarantee of correctly rounded results.
## Conversions and formatting
@@ -165,8 +189,8 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu
## Deferred scope
-The planned `PreciseMath` static class and its `Abs`, `Sqrt`, `Pow`, `Exp`, and
-`Log` functions are not implemented. Generic-math interfaces beyond `ISignedNumber`,
+The `DDMath.Pow`, `DDMath.Exp`, and `DDMath.Log` functions remain unimplemented.
+Generic-math interfaces beyond `ISignedNumber`,
additional text formats/general round-trip formatting, and non-arithmetic performance
benchmarks remain deferred.
Replacing allocating arithmetic boundary fallbacks is also deferred; the current