Skip to content

Commit b29c945

Browse files
thedavidmeisterthedavidmeister
committed
full precision div
1 parent 040b821 commit b29c945

8 files changed

Lines changed: 195 additions & 35 deletions

‎src/error/ErrDecimalFloat.sol‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,3 +30,6 @@ error LossyConversionFromFloat(int256 signedCoefficient, int256 exponent);
3030

3131
/// @dev Thrown when attempting to exponentiate 0^b where b is negative.
3232
error ZeroNegativePower(Float b);
33+
34+
/// @dev Thrown when muldiv internal to division overflows.
35+
error MulDivOverflow(uint256 x, uint256 y, uint256 denominator);

‎src/lib/implementation/LibDecimalFloatImplementation.sol‎

Lines changed: 148 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
// SPDX-License-Identifier: CAL
22
pragma solidity ^0.8.25;
33

4-
import {ExponentOverflow, Log10Negative, Log10Zero} from "../../error/ErrDecimalFloat.sol";
4+
import {ExponentOverflow, Log10Negative, Log10Zero, MulDivOverflow} from "../../error/ErrDecimalFloat.sol";
55
import {
66
LOG_TABLES,
77
LOG_TABLES_SMALL,
@@ -10,6 +10,7 @@ import {
1010
ANTI_LOG_TABLES_SMALL
1111
} from "../../generated/LogTables.pointers.sol";
1212
import {LibDecimalFloat} from "../LibDecimalFloat.sol";
13+
import {console2} from "forge-std/console2.sol";
1314

1415
error WithTargetExponentOverflow(int256 signedCoefficient, int256 exponent, int256 targetExponent);
1516

@@ -61,6 +62,11 @@ int256 constant MAXIMIZED_ZERO_SIGNED_COEFFICIENT = NORMALIZED_ZERO_SIGNED_COEFF
6162
/// @dev The exponent of maximized zero.
6263
int256 constant MAXIMIZED_ZERO_EXPONENT = NORMALIZED_ZERO_EXPONENT;
6364

65+
/// @dev The signed coefficient of minimized zero.
66+
int256 constant MINIMIZED_ZERO_SIGNED_COEFFICIENT = 0;
67+
/// @dev The exponent of minimized zero.
68+
int256 constant MINIMIZED_ZERO_EXPONENT = 0;
69+
6470
library LibDecimalFloatImplementation {
6571
/// Negates and normalizes a float.
6672
/// Equivalent to `0 - x`.
@@ -190,14 +196,129 @@ library LibDecimalFloatImplementation {
190196
{
191197
unchecked {
192198
(signedCoefficientA, exponentA) = maximize(signedCoefficientA, exponentA);
193-
(signedCoefficientB, exponentB) = normalize(signedCoefficientB, exponentB);
199+
(signedCoefficientB, exponentB) = maximize(signedCoefficientB, exponentB);
194200

195-
int256 signedCoefficient = signedCoefficientA / signedCoefficientB;
196-
int256 exponent = exponentA - exponentB;
201+
uint256 signedCoefficientAAbs;
202+
if (signedCoefficientA < 0) {
203+
if (signedCoefficientA == type(int256).min) {
204+
signedCoefficientAAbs = uint256(type(int256).max) + 1;
205+
} else {
206+
signedCoefficientAAbs = uint256(-signedCoefficientA);
207+
}
208+
} else {
209+
signedCoefficientAAbs = uint256(signedCoefficientA);
210+
}
211+
uint256 signedCoefficientBAbs;
212+
if (signedCoefficientB < 0) {
213+
if (signedCoefficientB == type(int256).min) {
214+
signedCoefficientBAbs = uint256(type(int256).max) + 1;
215+
} else {
216+
signedCoefficientBAbs = uint256(-signedCoefficientB);
217+
}
218+
} else {
219+
signedCoefficientBAbs = uint256(signedCoefficientB);
220+
}
221+
int256 scale = 1e76;
222+
int256 adjustExponent = 76;
223+
if (signedCoefficientB / scale == 0) {
224+
scale = 1e75;
225+
adjustExponent = 75;
226+
}
227+
uint256 signedCoefficientAbs = mulDiv(signedCoefficientAAbs, uint256(scale), signedCoefficientBAbs);
228+
int256 signedCoefficient = (signedCoefficientA ^ signedCoefficientB) < 0
229+
? -int256(signedCoefficientAbs)
230+
: int256(signedCoefficientAbs);
231+
int256 exponent = exponentA - exponentB - adjustExponent;
197232
return (signedCoefficient, exponent);
198233
}
199234
}
200235

236+
/// mulDiv as seen in Open Zeppelin, PRB Math, Solady, and other libraries.
237+
/// Credit to Remco Bloemen under MIT license: https://2π.com/21/muldiv
238+
function mulDiv(uint256 x, uint256 y, uint256 denominator) internal pure returns (uint256 result) {
239+
// 512-bit multiply [prod1 prod0] = x * y. Compute the product mod 2^256 and mod 2^256 - 1, then use
240+
// use the Chinese Remainder Theorem to reconstruct the 512-bit result. The result is stored in two 256
241+
// variables such that product = prod1 * 2^256 + prod0.
242+
uint256 prod0; // Least significant 256 bits of the product
243+
uint256 prod1; // Most significant 256 bits of the product
244+
assembly ("memory-safe") {
245+
let mm := mulmod(x, y, not(0))
246+
prod0 := mul(x, y)
247+
prod1 := sub(sub(mm, prod0), lt(mm, prod0))
248+
}
249+
250+
// Handle non-overflow cases, 256 by 256 division.
251+
if (prod1 == 0) {
252+
unchecked {
253+
return prod0 / denominator;
254+
}
255+
}
256+
257+
// Make sure the result is less than 2^256. Also prevents denominator == 0.
258+
if (prod1 >= denominator) {
259+
revert MulDivOverflow(x, y, denominator);
260+
}
261+
262+
////////////////////////////////////////////////////////////////////////////
263+
// 512 by 256 division
264+
////////////////////////////////////////////////////////////////////////////
265+
266+
// Make division exact by subtracting the remainder from [prod1 prod0].
267+
uint256 remainder;
268+
assembly ("memory-safe") {
269+
// Compute remainder using the mulmod Yul instruction.
270+
remainder := mulmod(x, y, denominator)
271+
272+
// Subtract 256 bit number from 512-bit number.
273+
prod1 := sub(prod1, gt(remainder, prod0))
274+
prod0 := sub(prod0, remainder)
275+
}
276+
277+
unchecked {
278+
// Calculate the largest power of two divisor of the denominator using the unary operator ~. This operation cannot overflow
279+
// because the denominator cannot be zero at this point in the function execution. The result is always >= 1.
280+
// For more detail, see https://cs.stackexchange.com/q/138556/92363.
281+
uint256 lpotdod = denominator & (~denominator + 1);
282+
uint256 flippedLpotdod;
283+
284+
assembly ("memory-safe") {
285+
// Factor powers of two out of denominator.
286+
denominator := div(denominator, lpotdod)
287+
288+
// Divide [prod1 prod0] by lpotdod.
289+
prod0 := div(prod0, lpotdod)
290+
291+
// Get the flipped value `2^256 / lpotdod`. If the `lpotdod` is zero, the flipped value is one.
292+
// `sub(0, lpotdod)` produces the two's complement version of `lpotdod`, which is equivalent to flipping all the bits.
293+
// However, `div` interprets this value as an unsigned value: https://ethereum.stackexchange.com/q/147168/24693
294+
flippedLpotdod := add(div(sub(0, lpotdod), lpotdod), 1)
295+
}
296+
297+
// Shift in bits from prod1 into prod0.
298+
prod0 |= prod1 * flippedLpotdod;
299+
300+
// Invert denominator mod 2^256. Now that denominator is an odd number, it has an inverse modulo 2^256 such
301+
// that denominator * inv = 1 mod 2^256. Compute the inverse by starting with a seed that is correct for
302+
// four bits. That is, denominator * inv = 1 mod 2^4.
303+
uint256 inverse = (3 * denominator) ^ 2;
304+
305+
// Use the Newton-Raphson iteration to improve the precision. Thanks to Hensel's lifting lemma, this also works
306+
// in modular arithmetic, doubling the correct bits in each step.
307+
inverse *= 2 - denominator * inverse; // inverse mod 2^8
308+
inverse *= 2 - denominator * inverse; // inverse mod 2^16
309+
inverse *= 2 - denominator * inverse; // inverse mod 2^32
310+
inverse *= 2 - denominator * inverse; // inverse mod 2^64
311+
inverse *= 2 - denominator * inverse; // inverse mod 2^128
312+
inverse *= 2 - denominator * inverse; // inverse mod 2^256
313+
314+
// Because the division is now exact we can divide by multiplying with the modular inverse of denominator.
315+
// This will give us the correct result modulo 2^256. Since the preconditions guarantee that the outcome is
316+
// less than 2^256, this is the final result. We don't need to compute the high bits of the result and prod1
317+
// is no longer required.
318+
result = prod0 * inverse;
319+
}
320+
}
321+
201322
/// Add two floats together.
202323
///
203324
/// Note that because the input values can have arbitrary exponents that may
@@ -367,14 +488,9 @@ library LibDecimalFloatImplementation {
367488
return signedCoefficientA == signedCoefficientB;
368489
}
369490

370-
/// Inverts a float. Equivalent to `1 / x` with modest gas optimizations.
491+
/// Inverts a float. Equivalent to `1 / x`.
371492
function inv(int256 signedCoefficient, int256 exponent) internal pure returns (int256, int256) {
372-
(signedCoefficient, exponent) = normalize(signedCoefficient, exponent);
373-
374-
signedCoefficient = 1e76 / signedCoefficient;
375-
exponent = -exponent - 76;
376-
377-
return (signedCoefficient, exponent);
493+
return div(1e76, -76, signedCoefficient, exponent);
378494
}
379495

380496
/// log10(x) for a float x.
@@ -535,6 +651,27 @@ library LibDecimalFloatImplementation {
535651
}
536652
}
537653

654+
// function minimize(int256 signedCoefficient, int256 exponent) internal pure returns (int256, int256) {
655+
// unchecked {
656+
// if (signedCoefficient == 0) {
657+
// return (MINIMIZED_ZERO_SIGNED_COEFFICIENT, MINIMIZED_ZERO_EXPONENT);
658+
// }
659+
660+
// int256 initialExponent = exponent;
661+
662+
// while (signedCoefficient % 10 == 0) {
663+
// signedCoefficient /= 10;
664+
// exponent += 1;
665+
// }
666+
667+
// if (initialExponent > exponent) {
668+
// revert ExponentOverflow(signedCoefficient, exponent);
669+
// }
670+
671+
// return (signedCoefficient, exponent);
672+
// }
673+
// }
674+
538675
function maximize(int256 signedCoefficient, int256 exponent) internal pure returns (int256, int256) {
539676
unchecked {
540677
if (signedCoefficient == 0) {

‎test/lib/LibCommonResults.sol‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
// SPDX-License-Identifier: CAL
22
pragma solidity ^0.8.25;
33

4-
int256 constant ONES = 111111111111111111111111111111111111111;
5-
int256 constant THREES = 333333333333333333333333333333333333333;
4+
int256 constant ONES = 1111111111111111111111111111111111111111111111111111111111111111111111111111;
5+
int256 constant THREES_PACKED = 3333333333333333333333333333333333333333333333333333333333333333333;
6+
int256 constant THREES = 3333333333333333333333333333333333333333333333333333333333333333333333333333;

‎test/src/lib/LibDecimalFloat.mixed.t.sol‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -2,26 +2,26 @@
22
pragma solidity =0.8.25;
33

44
import {LibDecimalFloat, Float} from "src/lib/LibDecimalFloat.sol";
5-
import {THREES, ONES} from "../../lib/LibCommonResults.sol";
5+
import {THREES_PACKED, ONES} from "../../lib/LibCommonResults.sol";
66

77
import {Test} from "forge-std/Test.sol";
88

99
contract LibDecimalFloatMixedTest is Test {
1010
using LibDecimalFloat for Float;
1111

1212
/// (1 / 3) * 555e18
13-
function testDiv1Over3() external pure {
13+
function testDiv1Over3Mixed() external pure {
1414
Float a = LibDecimalFloat.packLossless(1, 0);
1515
Float b = LibDecimalFloat.packLossless(3, 0);
1616
Float c = a.div(b);
1717
(int256 signedCoefficientDiv, int256 exponentDiv) = LibDecimalFloat.unpack(c);
18-
assertEq(signedCoefficientDiv, THREES, "coefficient");
19-
assertEq(exponentDiv, -39, "exponent");
18+
assertEq(signedCoefficientDiv, THREES_PACKED, "coefficient");
19+
assertEq(exponentDiv, -67, "exponent");
2020

2121
Float d = c.mul(LibDecimalFloat.packLossless(555, 18));
2222
(int256 signedCoefficientMul, int256 exponentMul) = LibDecimalFloat.unpack(d);
2323

24-
assertEq(signedCoefficientMul, 184999999999999999999999999999999999999815);
25-
assertEq(exponentMul, -21);
24+
assertEq(signedCoefficientMul, 1849999999999999999999999999999999999999999999999999999999999999999);
25+
assertEq(exponentMul, -46);
2626
}
2727
}

‎test/src/lib/LibDecimalFloat.pow.t.sol‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -38,9 +38,15 @@ contract LibDecimalFloatPowTest is LogTest {
3838
}
3939

4040
function testPows() external {
41-
checkPow(5e37, -38, 3e37, -36, 9.32835820895522388059701492537313432835e38, -48);
42-
checkPow(5e37, -38, 6e37, -36, 8.71080139372822299651567944250871080139e38, -57);
43-
// // Issues found in fuzzing from here.
41+
// 0.5 ^ 30 = 9.3132257e-10
42+
checkPow(
43+
5e37, -38, 3e37, -36, 9.328358208955223880597014925373134328358208955223880597014925373134e66, -66 - 10
44+
);
45+
// 0.5 ^ 60 = 8.6736174e-19
46+
checkPow(
47+
5e37, -38, 6e37, -36, 8.710801393728222996515679442508710801393728222996515679442508710801e66, -66 - 19
48+
);
49+
// Issues found in fuzzing from here.
4450
checkPow(99999, 0, 12182, 0, 1000, 60907);
4551
checkPow(1785215562, 0, 18, 0, 3388, 163);
4652
}

‎test/src/lib/LibDecimalFloat.pow10.t.sol‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,12 @@ contract LibDecimalFloatPow10Test is LogTest {
2727
} else {
2828
Float floatPower10 = this.pow10External(float);
2929
(int256 signedCoefficientUnpacked, int256 exponentUnpacked) = floatPower10.unpack();
30+
31+
// Compensate for the implied pack and unpack.
32+
(Float resultPacked, bool lossless) = LibDecimalFloat.packLossy(signedCoefficient, exponent);
33+
(lossless);
34+
(signedCoefficient, exponent) = resultPacked.unpack();
35+
3036
assertEq(signedCoefficient, signedCoefficientUnpacked);
3137
assertEq(exponent, exponentUnpacked);
3238
}

‎test/src/lib/implementation/LibDecimalFloatImplementation.div.t.sol‎

Lines changed: 17 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -40,12 +40,12 @@ contract LibDecimalFloatImplementationDivTest is Test {
4040

4141
/// 1 / 3
4242
function testDiv1Over3() external pure {
43-
checkDiv(1, 0, 3, 0, THREES, -39);
43+
checkDiv(1, 0, 3, 0, THREES, -76);
4444
}
4545

4646
/// - 1 / 3
4747
function testDivNegative1Over3() external pure {
48-
checkDiv(-1, 0, 3, 0, -THREES, -39);
48+
checkDiv(-1, 0, 3, 0, -THREES, -76);
4949
}
5050

5151
/// 1 / 3 gas
@@ -56,41 +56,47 @@ contract LibDecimalFloatImplementationDivTest is Test {
5656

5757
/// 1e18 / 3
5858
function testDiv1e18Over3() external pure {
59-
checkDiv(1e18, 0, 3, 0, THREES, -21);
59+
checkDiv(1e18, 0, 3, 0, THREES, -58);
6060
}
6161

6262
/// 10,0 / 1e38,-37 == 1
6363
function testDivTenOverOOMs() external pure {
64-
checkDiv(10, 0, 1e38, -37, 1e39, -39);
64+
checkDiv(10, 0, 1e38, -37, 1e76, -76);
6565
}
6666

6767
/// 1e38,-37 / 2,0 == 5
6868
function testDivOOMsOverTen() external pure {
69-
checkDiv(1e38, -37, 2, 0, 5e38, -38);
69+
checkDiv(1e38, -37, 2, 0, 5e75, -75);
7070
}
7171

7272
/// 5e37,-37 / 2e37,-37 == 2.5
7373
function testDivOOMs5and2() external pure {
74-
checkDiv(5e37, -37, 2e37, -37, 25e38, -39);
74+
checkDiv(5e37, -37, 2e37, -37, 2.5e76, -76);
7575
}
7676

7777
/// (1 / 9) / (1 / 3) == 0.333..
7878
function testDiv1Over9Over1Over3() external pure {
7979
// 1 / 9
8080
(int256 signedCoefficientA, int256 exponentA) = LibDecimalFloatImplementation.div(1, 0, 9, 0);
8181
assertEq(signedCoefficientA, ONES);
82-
assertEq(exponentA, -39);
82+
assertEq(exponentA, -76);
8383

8484
// 1 / 3
8585
(int256 signedCoefficientB, int256 exponentB) = LibDecimalFloatImplementation.div(1, 0, 3, 0);
8686
assertEq(signedCoefficientB, THREES);
87-
assertEq(exponentB, -39);
87+
assertEq(exponentB, -76);
8888

8989
// (1 / 9) / (1 / 3)
9090
(int256 signedCoefficient, int256 exponent) =
9191
LibDecimalFloatImplementation.div(signedCoefficientA, exponentA, signedCoefficientB, exponentB);
92-
assertEq(signedCoefficient, 333333333333333333333333333333333333336);
93-
assertEq(exponent, -39);
92+
assertEq(signedCoefficient, THREES);
93+
assertEq(exponent, -76);
94+
95+
// (1 / 3) / (1 / 9) == 3
96+
(signedCoefficient, exponent) =
97+
LibDecimalFloatImplementation.div(signedCoefficientB, exponentB, signedCoefficientA, exponentA);
98+
assertEq(signedCoefficient, 3e76);
99+
assertEq(exponent, -76);
94100
}
95101

96102
/// forge-config: default.fuzz.runs = 100
@@ -102,7 +108,7 @@ contract LibDecimalFloatImplementationDivTest is Test {
102108
int256 di = 0;
103109
while (true) {
104110
int256 i = 1;
105-
int256 j = -39 - di;
111+
int256 j = -76 - di;
106112
while (true) {
107113
// want to see full precision on the THREES regardless of the
108114
// scale of the numerator and denominator.

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

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ contract LibDecimalFloatImplementationPow10Test is LogTest {
3232
checkPow10(1, 4, 1000, 9997);
3333
}
3434

35-
function testExactLookups() external {
35+
function testExactLookupsPow10() external {
3636
// 10^2 = 100
3737
checkPow10(2, 0, 1000, -1);
3838
// 10^3 = 1000
@@ -56,7 +56,8 @@ contract LibDecimalFloatImplementationPow10Test is LogTest {
5656
checkPow10(0.5e37, -37, 3162, -3);
5757

5858
checkPow10(0.3e37, -37, 1995, -3);
59-
checkPow10(-0.3e37, -37, 0.501253132832080200501253132832080200501e39, -39);
59+
// 10^-0.3 = 0.50118723362
60+
checkPow10(-0.3e37, -37, 0.5012531328320802005012531328320802005012531328320802005012531328320802005012e76, -76);
6061
}
6162

6263
function testInterpolatedLookupsPower() external {

0 commit comments

Comments
 (0)