Check pow10 and log10 against the oracle, not the unrounded functions

README.md

@@ -161,8 +161,10 @@ additional judgement calls re: when precision loss is acceptable.
 
 #### log10, pow10, pow and sqrt
 
-These round their final result to nearest at 41 significant digits, `pow`
-including an integer power and `a^1`, within proven bounds of the true value:
+These round an internal approximation to nearest at 41 significant digits, `pow`
+including an integer power and `a^1`, within proven bounds of the true value.
+`sqrt` is correctly rounded; for the others, a true value within the bound of a
+rounding tie can round to either neighbour:
 
 - `pow10`: half a unit in the 41st digit plus 3.28e-8 of a unit, under
   5.0000004e-41 relative.

src/lib/implementation/LibDecimalFloatImplementation.sol

@@ -878,9 +878,11 @@ library LibDecimalFloatImplementation {
         return div(1e76, -76, signedCoefficient, exponent);
     }
 
-    /// log10(x) for a float x, rounded to nearest at 41 significant digits,
-    /// half away from zero, so within half a unit in the 41st digit plus
-    /// `LOG10_RAW_ERROR` units of 1e-50. log10(10^k) is exactly k.
+    /// log10(x) for a float x: an internal approximation within
+    /// `LOG10_RAW_ERROR` units of 1e-50, rounded to nearest at 41 significant
+    /// digits half away from zero, so within half a unit in the 41st digit
+    /// plus that error. Near a tie either neighbour can result. log10(10^k) is
+    /// exactly k.
     ///
     /// @param signedCoefficient The signed coefficient of the floating point
     /// number.
@@ -1376,8 +1378,9 @@ library LibDecimalFloatImplementation {
 
     // slither-disable-end too-many-digits
 
-    /// 10^x for a float x, rounded to nearest at 41 significant digits, half
-    /// up. 10^k is exactly 10^k for an integer k.
+    /// 10^x for a float x: an internal approximation within the bound below,
+    /// rounded to nearest at 41 significant digits half up. Near a tie either
+    /// neighbour can result. 10^k is exactly 10^k for an integer k.
     ///
     /// The fraction m of x, truncated to 1e-50, goes through `exp10Fixed`.
     /// Truncation moves 10^m by under 2.3026e-50 relative either way, and with

test/lib/LibTranscendentalOracle.sol

@@ -136,8 +136,8 @@ library LibTranscendentalOracle {
     }
 
     /// 10^(signedCoefficient 10^exponent) as power 10^powerExponent, power in
-    /// [ORACLE_ONE, 10 ORACLE_ONE], within 1e-67 relative, for an exponent of
-    /// at most 10 and an integer part that fits.
+    /// [ORACLE_ONE, 10 ORACLE_ONE], within 1e-67 relative, for an integer part
+    /// that fits int256.
     ///
     /// Relative, in units of 1e-70: the fraction f floors under a unit, 2.31
     /// of the power. y = f ln 10 / 16, at most 0.144, is within 2 below, so

test/src/lib/implementation/LibDecimalFloatImplementation.log10.t.sol

@@ -9,6 +9,7 @@ import {LibTranscendentalOracle, ORACLE_ONE, ORACLE_LN10} from "../../../lib/Lib
 import {Math} from "@openzeppelin-contracts-5.7.0/utils/math/Math.sol";
 import {LibDecimalFloat, Float} from "src/lib/LibDecimalFloat.sol";
 import {LibTestErrorBound} from "../../../lib/LibTestErrorBound.sol";
+import {LibTestExactDecimal, U512} from "../../../lib/LibTestExactDecimal.sol";
 
 contract LibDecimalFloatImplementationLog10Test is Test {
     function checkLog10(
@@ -645,24 +646,107 @@ contract LibDecimalFloatImplementationLog10Test is Test {
         assertLe(log10NearOneUlpError(p, d, negative), 500000032, "log10 error");
     }
 
-    /// log10 is log10Unrounded rounded to nearest, so within half a unit.
-    function testLog10RoundsUnrounded(int256 signedCoefficient, int256 exponent) external pure {
+    /// |a - b| for magnitudes with signs.
+    function signedDistance(bool negativeA, U512 memory a, bool negativeB, U512 memory b)
+        internal
+        pure
+        returns (U512 memory)
+    {
+        if (negativeA != negativeB) {
+            return LibTestExactDecimal.add(a, b);
+        }
+        return LibTestExactDecimal.cmp(a, b) >= 0 ? subU512(a, b) : subU512(b, a);
+    }
+
+    /// a - b for a >= b.
+    function subU512(U512 memory a, U512 memory b) internal pure returns (U512 memory) {
+        unchecked {
+            uint256 lo = a.lo - b.lo;
+            uint256 borrow = a.lo < b.lo ? 1 : 0;
+            return U512(a.hi - b.hi - borrow, lo);
+        }
+    }
+
+    /// 10^n as a U512.
+    function pow10U512(uint256 n) internal pure returns (U512 memory) {
+        return LibTestExactDecimal.mulPow10(LibTestExactDecimal.u512(1), n);
+    }
+
+    /// Whether a log10 or log10Unrounded result for a positive input is
+    /// within the documented bound of the true log, measured against the
+    /// oracle. All in units of 1e-70:
+    /// - the oracle is within 265 of the true log L;
+    /// - the README's 2e-50 absolute, 2e20;
+    /// - past a characteristic of 1e25 the sum goes through `add`, which
+    ///   loses under a unit of an int256 coefficient's last digit, at most
+    ///   10^(n - 76) for a characteristic below 10^n;
+    /// - the result is floored onto the 1e-70 grid, a unit;
+    /// - rounded, half a unit in the 41st digit of L, whose magnitude is
+    ///   below the oracle's plus 265.
+    function log10WithinOracle(
+        int256 signedCoefficient,
+        int256 exponent,
+        int256 resultCoefficient,
+        int256 resultExponent,
+        bool rounded
+    ) internal pure returns (bool) {
+        // forge-lint: disable-next-line(unsafe-typecast)
+        (int256 characteristic, uint256 fraction) = LibTranscendentalOracle.log10(uint256(signedCoefficient), exponent);
+        bool negativeExpected = characteristic < 0;
+        U512 memory characteristic70 =
+            LibTestExactDecimal.mulPow10(LibTestExactDecimal.u512(LibTestExactDecimal.abs(characteristic)), 70);
+        U512 memory expected = negativeExpected
+            ? subU512(characteristic70, LibTestExactDecimal.u512(fraction))
+            : LibTestExactDecimal.add(characteristic70, LibTestExactDecimal.u512(fraction));
+
+        U512 memory actual;
+        uint256 resultMagnitude = LibTestExactDecimal.abs(resultCoefficient);
+        if (resultExponent >= -70) {
+            // Non-negative in this branch.
+            // forge-lint: disable-next-line(unsafe-typecast)
+            uint256 shift = uint256(resultExponent + 70);
+            actual = LibTestExactDecimal.mulPow10(LibTestExactDecimal.u512(resultMagnitude), shift);
+        } else if (resultExponent >= -70 - 77) {
+            // forge-lint: disable-next-line(unsafe-typecast)
+            actual = LibTestExactDecimal.u512(resultMagnitude / 10 ** uint256(-70 - resultExponent));
+        }
+        U512 memory error = signedDistance(resultCoefficient < 0, actual, negativeExpected, expected);
+
+        U512 memory bound = LibTestExactDecimal.u512(2e20 + 265 + 1);
+        uint256 characteristicDigits =
+            uint256(LibTestExactDecimal.digits(LibTestExactDecimal.u512(LibTestExactDecimal.abs(characteristic) + 1)));
+        if (characteristicDigits > 25) {
+            bound = LibTestExactDecimal.add(bound, pow10U512(characteristicDigits - 6));
+        }
+        if (rounded) {
+            uint256 logDigits =
+                uint256(LibTestExactDecimal.digits(LibTestExactDecimal.add(expected, LibTestExactDecimal.u512(265))));
+            bound = LibTestExactDecimal.add(
+                bound,
+                logDigits >= 42
+                    ? LibTestExactDecimal.mulSmall(pow10U512(logDigits - 42), 5)
+                    : LibTestExactDecimal.u512(1)
+            );
+        }
+        return LibTestExactDecimal.cmp(error, bound) <= 0;
+    }
+
+    /// log10 is within half a unit in the 41st digit of the true log plus the
+    /// README's 2e-50, against the oracle, over every positive coefficient
+    /// and exponents past the 1e25 characteristic that leaves the fixed point
+    /// sum.
+    function testLog10WithinOracle(int256 signedCoefficient, int256 exponent) external pure {
         signedCoefficient = bound(signedCoefficient, 1, type(int256).max);
         exponent = bound(exponent, -1e40, 1e40);
-        (int256 unroundedCoefficient, int256 unroundedExponent) =
-            LibDecimalFloatImplementation.log10Unrounded(signedCoefficient, exponent);
         (int256 logCoefficient, int256 logExponent) = LibDecimalFloatImplementation.log10(signedCoefficient, exponent);
-        if (logCoefficient == 0) {
-            assertEq(unroundedCoefficient, 0, "zero");
-            return;
-        }
-        assertLe(ulpBillionths(logCoefficient, logExponent, unroundedCoefficient, unroundedExponent), 500000000, "half");
+        assertTrue(log10WithinOracle(signedCoefficient, exponent, logCoefficient, logExponent, true), "log10 oracle");
     }
 
     /// log10Unrounded(x 10^k) = log10Unrounded(x) + k exactly, as the
     /// characteristic is summed as an integer, except that a log within
     /// log10(1.001) of zero keeps digits below the 1e-50 its shift truncates
     /// to.
+    /// Both are within the bound of the oracle.
     function testLog10DecadeShift(int256 signedCoefficient, int256 exponent, int256 shift) external pure {
         signedCoefficient = bound(signedCoefficient, 1, type(int256).max);
         exponent = bound(exponent, -1e18, 1e18);
@@ -671,6 +755,11 @@ contract LibDecimalFloatImplementationLog10Test is Test {
             LibDecimalFloatImplementation.log10Unrounded(signedCoefficient, exponent);
         (int256 shiftedCoefficient, int256 shiftedExponent) =
             LibDecimalFloatImplementation.log10Unrounded(signedCoefficient, exponent + shift);
+        assertTrue(log10WithinOracle(signedCoefficient, exponent, logCoefficient, logExponent, false), "log10 oracle");
+        assertTrue(
+            log10WithinOracle(signedCoefficient, exponent + shift, shiftedCoefficient, shiftedExponent, false),
+            "shifted log10 oracle"
+        );
         (logCoefficient, logExponent) = LibDecimalFloatImplementation.add(logCoefficient, logExponent, shift, 0);
         if (logExponent < -50 || shiftedExponent < -50) {
             (int256 errorCoefficient, int256 errorExponent) =
@@ -717,13 +806,21 @@ contract LibDecimalFloatImplementationLog10Test is Test {
 
     /// x < y implies log10(x) <= log10(y) + 2E, down to adjacent 76 digit
     /// coefficients, whose logs differ by far less than a unit in the last place.
+    /// Both are within the bound of the oracle.
     function testLog10Monotone(int256 signedCoefficient, int256 gap, int256 exponent) external pure {
         signedCoefficient = bound(signedCoefficient, 1e75, 1e76 - 1e3);
         gap = bound(gap, 1, 1e3);
         exponent = bound(exponent, -1e6, 1e6);
         (int256 lowCoefficient, int256 lowExponent) = LibDecimalFloatImplementation.log10(signedCoefficient, exponent);
         (int256 highCoefficient, int256 highExponent) =
             LibDecimalFloatImplementation.log10(signedCoefficient + gap, exponent);
+        assertTrue(
+            log10WithinOracle(signedCoefficient, exponent, lowCoefficient, lowExponent, true), "low log10 oracle"
+        );
+        assertTrue(
+            log10WithinOracle(signedCoefficient + gap, exponent, highCoefficient, highExponent, true),
+            "high log10 oracle"
+        );
         Float low = LibDecimalFloat.packLossless(lowCoefficient, lowExponent);
         Float high = LibDecimalFloat.packLossless(highCoefficient, highExponent);
         assertTrue(

test/src/lib/implementation/LibDecimalFloatImplementation.pow10.t.sol

@@ -364,10 +364,74 @@ contract LibDecimalFloatImplementationPow10Test is Test {
         return (x, exponent);
     }
 
-    /// Test the current range that we can handle power10 over does not revert.
+    /// Whether a pow10 result is within the README's bound of a power P 10^e
+    /// that is itself within `referenceError` units of 10^e of the true
+    /// power T. In units of 10^e, T is below 10^n for n the digits of
+    /// P + referenceError, so a unit in T's 41st digit is at most 10^(n - 41),
+    /// and the bound is half of it plus 3.28e-8 of it.
+    function pow10WithinPower(
+        uint256 power,
+        int256 powerExponent,
+        uint256 referenceError,
+        int256 resultCoefficient,
+        int256 resultExponent
+    ) internal pure returns (bool) {
+        // P has 70 to 72 digits, a result within bound at most 41, or 1 at a
+        // power of ten.
+        if (resultCoefficient <= 0 || resultExponent < powerExponent || resultExponent > powerExponent + 71) {
+            return false;
+        }
+        // forge-lint: disable-next-line(unsafe-typecast)
+        uint256 actual = uint256(resultCoefficient) * 10 ** uint256(resultExponent - powerExponent);
+        uint256 n = 0;
+        for (uint256 rest = power + referenceError; rest > 0; rest /= 10) {
+            n++;
+        }
+        uint256 bound = 5 * 10 ** (n - 42) + 328 * 10 ** (n - 51) + referenceError;
+        return (actual > power ? actual - power : power - actual) <= bound;
+    }
+
+    /// `pow10WithinPower` against the oracle, which is within 1e-67 relative.
+    function pow10WithinOracle(
+        int256 signedCoefficient,
+        int256 exponent,
+        int256 resultCoefficient,
+        int256 resultExponent
+    ) internal pure returns (bool) {
+        (uint256 power, int256 powerExponent) = LibTranscendentalOracle.exp10(signedCoefficient, exponent);
+        return pow10WithinPower(power, powerExponent, power / 1e67 + 1, resultCoefficient, resultExponent);
+    }
+
+    /// Over the whole range pow10 handles, it does not revert and is within
+    /// its bound of the oracle.
     function testNoRevert(int224 x, int32 exponent) external pure {
         (x, exponent) = boundFloat(x, exponent);
-        LibDecimalFloatImplementation.pow10(x, exponent);
+        (int256 signedCoefficient, int256 resultExponent) = LibDecimalFloatImplementation.pow10(x, exponent);
+        assertTrue(pow10WithinOracle(x, exponent, signedCoefficient, resultExponent), "pow10 oracle");
+    }
+
+    /// The `testPow10RoundsHalfUp` inputs, whose true powers sit within 2e-9
+    /// of a unit of a tie at the 41st digit, are within the bound of their
+    /// 70 digit powers from `bc -l` at scale 220. The last is above its tie
+    /// and rounds down, which the bound allows.
+    function testPow10NearTiesWithinBound() external pure {
+        int256[2][4] memory references = [
+            [int256(10857362048), 1000000000000000000000000000000000000000250000000009641022996351405692],
+            [int256(67315644695), 1000000000000000000000000000000000000001549999999999907130159533527510],
+            [
+                int256(52839762350034614138096464404871856380432525744962),
+                3375962574923628115885758764340189836032550000000186267648444008762503
+            ],
+            [
+                int256(62876337819492131669053461917963730523385842721800),
+                4253665924693161986466925885576001212321550000000178414220893221038845
+            ]
+        ];
+        for (uint256 i = 0; i < references.length; i++) {
+            (int256 signedCoefficient, int256 exponent) = LibDecimalFloatImplementation.pow10(references[i][0], -50);
+            // forge-lint: disable-next-line(unsafe-typecast)
+            assertTrue(pow10WithinPower(uint256(references[i][1]), -69, 1, signedCoefficient, exponent), "pow10 bc");
+        }
     }
 
     /// Guard digits of exactly half a unit round up, and one below half
@@ -448,13 +512,18 @@ contract LibDecimalFloatImplementationPow10Test is Test {
     /// pow10(x + k) is pow10(x) 10^k exactly for an integer k. x is at the
     /// 1e-50 that pow10 truncates its fraction to, so the shift does not move
     /// the truncation.
+    /// Both are within the bound of the oracle.
     function testPow10DecadeShift(int256 x, int256 shift) external pure {
         x = bound(x, -1e55, 1e55);
         shift = bound(shift, -1e5, 1e5);
         (int256 signedCoefficient, int256 exponent) = LibDecimalFloatImplementation.pow10(x, -50);
         (int256 shiftedCoefficient, int256 shiftedExponent) = LibDecimalFloatImplementation.pow10(x + shift * 1e50, -50);
         assertEq(shiftedCoefficient, signedCoefficient, "coefficient");
         assertEq(shiftedExponent, exponent + shift, "exponent");
+        assertTrue(pow10WithinOracle(x, -50, signedCoefficient, exponent), "pow10 oracle");
+        assertTrue(
+            pow10WithinOracle(x + shift * 1e50, -50, shiftedCoefficient, shiftedExponent), "shifted pow10 oracle"
+        );
     }
 
     /// pow10(k) is exactly 10^k for an integer k however it is written.
@@ -467,6 +536,7 @@ contract LibDecimalFloatImplementationPow10Test is Test {
 
     /// x < y implies pow10(x) <= pow10(y) + 2E, down to adjacent inputs at
     /// 1e-70.
+    /// Both are within the bound of the oracle.
     function testPow10Monotone(int256 x, int256 gap) external pure {
         x = bound(x, -1e75, 1e75);
         gap = bound(gap, 1, 1e3);
@@ -481,5 +551,7 @@ contract LibDecimalFloatImplementationPow10Test is Test {
             ),
             "monotone"
         );
+        assertTrue(pow10WithinOracle(x, -70, lowCoefficient, lowExponent), "low pow10 oracle");
+        assertTrue(pow10WithinOracle(x + gap, -70, highCoefficient, highExponent), "high pow10 oracle");
     }
 }